Petr Sladek (sladep1) & Milan Kratochvil (kratom4) at CTU Prague
The script includes two tasks: myekv and lineq (linekv) you can see comparison of both methods the 'histrogram equalization' and the 'quasi-linear (with saturation) manual brightness correction' files: myekv.m, lineq.m
% clear everything clear all
read, if RGB, automatically grayscale.
img = imread('Vychazi_mesic_2.jpg'); % rgb or grayscale input? imgdim=size(size(img),2); if imgdim==2 orig.img = round(double(img)); fprintf('grayscale picture'); end; if imgdim==3 orig.img = round(double(rgb2gray(img))); orig.img(any(orig.img>=255))=255; fprintf('RGB picture, grayscaled'); end;
RGB picture, grayscaled
pixelsN = prod(size(orig.img)); orig.histm =hist(double(orig.img),0:255)./pixelsN; %generate vector histogram 0..255 not 256x256 matrix: %means, can be done with index matrix... % 1 + 4 + 7... = orig.hist(1) % 2 + 5 + 8... = orig.hist(2) % 3 + 6 + 9... = orig.hist(3) [sy,sx]=size(orig.histm); onesi=ones(sx,1); orig.hist=orig.histm*onesi; % sum columns into one column. see also sum()
% q_k, q_0 q_k=orig.hist(256); q_0=orig.hist(1); % index matrix 1,2,...256 index=1:1:256; % create triangular matrix (p_0..p) mult.mask % 1 0 0 0 % 1 1 0 0 % 1 1 1 0 % 1 1 1 1 ... onesmask = tril(ones(256)); %equalized histogram: % $$ sum_{i=0}^{n} x_i $$ ekv.hist(index) = (orig.hist(:)')*(onesmask(index,:)'); figure(6);clf; plot(ekv.hist); title('histogram equalization look-up-table');
%allocate mem
ekv.img=ones(sy,sx);
ekv.img=ekv.hist(orig.img+1);
%generate manual correction lut (look-up-table) linear_with_saturation: a=0; b=60; c=2; d=253; lins.lut=lineq(a,b,c,d); figure(7);clf; plot(lins.lut); title('manual brightness correction look-up-table'); % normalize to image 1=>255. (grayscale 0..1) lins.lut = lins.lut.*(1/255); % correction by manual "lut" %allocate; lins.img = orig.img; lins.img = lins.lut(orig.img+1); %histograms: ekv.histm = hist(double(ekv.img)*256,0:255)./pixelsN; ekv.hist=ekv.histm*onesi; lins.histm = hist(double(lins.img)*256,0:255)./pixelsN; lins.hist=lins.histm*onesi;
figure(1); clf; %original image subplot(2,3,1); imshow(orig.img./255); title(sprintf('imput image %d x %d pix',sx,sy)); %and it's histogram: subplot(2,3,4); plot(orig.hist); title('input image histogram'); %equalized image subplot(2,3,2); imshow(ekv.img); title('histogram equalization'); %and it's histg subplot(2,3,5); plot(ekv.hist); title('image histogram'); % _ %linear _/ subplot(2,3,3); imshow(lins.img); title(sprintf('linear+saturation a=%d,b=%d,c=%d,d=%d',a,b,c,d)); %and it's histg subplot(2,3,6); plot(lins.hist); title('image histogram'); %detail view: figure(2);clf; imshow(orig.img./255); title(sprintf('Detail view: imput image %d x %d pix',sx,sy)); figure(3);clf; imshow(ekv.img); title('Detail view: histogram equalization'); figure(4);clf; imshow(lins.img); title(sprintf('Detail view: linear+saturation a=%d,b=%d,c=%d,d=%d',a,b,c,d));
% input parametrs: n = 1; Do = 140; %filter parameters a=0; b = 1; %size of image [im_height,im_width] = size(ekv.img); %euclidean distance matrix [U,V] = rc2uv(im_height,im_width); D = sqrt(U.^2+V.^2); %generate butterworth filter from euclidean distance matrix H = 1./(1+(D./Do).^(2*n)); %spectra of image IM = fftshift(fft2(ekv.img)); %filter image with HFE (matrix dot product) IMF=IM.*H; %image reconstruction from spectra IMF (inverse fft) imf = ifft2(fftshift(IMF)); %print output: figure(5); clf; imshow(imf); title(sprintf('soften image'));