Processing in MATLAB
Experiment 1 - See experiment page for description
% EXPERIMENT 1: SIGNALS AND FILTERING clear; % Create a sine wave modulated in pulses t = 0:0.001:1.023; % Time index s = sin(2*pi*30*t); % Sine wave at 30 Hz p = (1+square(2*pi*7*t))/2; % Square wave at 7 Hz (0/1 amplitudes only) n = s.*p; % Create train of sinusoidal pulses x = n + 0.1 * 2*(rand(1,1024)-.5); % Add some noise xerror=sum((x-n).^2); % Compute error as sum of the squares %plot(t,s+5,t,p+2,t,x); % Compute averaged signal (low pass filter) a = conv([1 2 4 2 1],x); % Average using 12421 filter a = a(3:1026)/10; % Extract central elements and normalize aerror=sum((a-n).^2); % Compute error m = medfilt1(x,5); merror=sum((m-n).^2); % Compute bilateral filter (low pass with edge preservation) xl1 = [0 x(1:1023)]; % Get values 1 and 2 units delay xl2 = [0 0 x(1:1022)]; xr1 = [x(2:1024) 0]; % Get values 1 and 2 units advance xr2 = [x(3:1024) 0 0]; mp = (.5) ^ 2; % Mean edge value squared ms = (1.5) ^2; % Mean spatial distance squared wl1 = exp(-(x-xl1).^2/mp) * exp(-1^2/ms); % Compute gaussian weights wl2 = exp(-(x-xl2).^2/mp) * exp(-2^2/ms); wr1 = exp(-(x-xr1).^2/mp) * exp(-1^2/ms); wr2 = exp(-(x-xr2).^2/mp) * exp(-2^2/ms); y = (xl2.*wl2 + xl1.*wl1 + x + xr1.*wr1 + xr2.*wr2); % Apply weights y = y./( wl2 + wl1 + 1 + wr1 + wr2 ); % Normalize yerror = sum((y-n).^2); % Compute error %plot(t,9+x,t,6+a,t,3+m,t,y); % Use Hamming window and compute magnitudes f=t(1:512); w=hamming(1024)'; xf=abs(fft(x.*w)); xf=log(xf(1:512)); nf=abs(fft(n.*w)); nf=log(nf(1:512)); af=abs(fft(a.*w)); af=log(af(1:512)); mf=abs(fft(m.*w)); mf=log(mf(1:512)); yf=abs(fft(y.*w)); yf=log(yf(1:512)); %plot(f,20+xf,f,15+af,f,10+mf,f,5+yf,f,nf); %Cepstrum calculations xc=abs(cceps(x)); nc=abs(cceps(n)); ac=abs(cceps(a)); mc=abs(cceps(m)); yc=abs(cceps(y)); plot(t,4+xc,t,3+ac,t,2+mc,t,1+yc,t,nc,1,-1);
Experiment 2 - See experiment page for description
% EXPERIMENT 2: SMOOTHING FILTERS - CREATE TEST FILES FROM PHOTO
% Initializations
clear, close all
N=256; % Define processing unit size
h=[1 2 4 2 1 ;
2 4 8 4 2 ;
4 8 16 8 4 ;
2 4 8 4 2 ;
1 2 4 2 1 ]; % Low-pass filter
h=h/sum(sum(h)); % Normalize
% Read main image
A=imread('cells.jpg'); % Read image file
X=(sum(A,3,'double')/3)/256; % Convert to grayscale
S=min([size(X,1) size(X,2)]); % Compute size of image
S=floor(S/N); % Compute scale factor
for iteration=1:floor(S/(size(h,1)-1))
X=filter2(h,X); % Low-pass filter before decimation
end
X=X(1:S:end,1:S:end); % Decimate in 2-D
S1=floor((size(X,1)-N)/2); % Compute crop area
S2=floor((size(X,2)-N)/2);
X=X((1+S1):(S1+N),(1+S2):(S2+N)); % Crop image
imwrite(X,'sfimg1.bmp'); % Write as low contrast image
% Amplify contrast
C1=min(min(X)); % Get lowest value
C2=max(max(X)); % Get highest value
X=(X-C1)/(C2-C1); % Rescale linearly
imwrite(X,'sfimg2.bmp'); % Write amplified image
image(cat(3,X,X,X)); % Show image
% EXPERIMENT 2: SMOOTHING FILTERS - PROCESS THE TEST FILES
% Initializations
clear, close all
dB_Scale=11; % Define magnitude log display scale
h=[1 2 4 2 1 ;
2 4 8 4 2 ;
4 8 16 8 4 ;
2 4 8 4 2 ;
1 2 4 2 1 ]; % Low-pass filter
h=h/sum(sum(h)); % Normalize
% Read main image
X=imread('sfimg2.bmp'); % Read image file
X=double(X)/256; % Convert to floating point
N=min(size(X)); % Get size of the array (must be even)
% Compute Fourier Transform
F=fft2(X); % Perform 2-D FFT
M=log(eps+abs(F))/dB_Scale; % Put magnitude on log scale
S=circshift(M,[N/2 N/2]); % Shift so DFT is centered
imwrite(S,'sfimg2ft.bmp'); % Save DFT image
S=0.5+angle(F)/(2*pi); % Get phase angles
S=circshift(S,[N/2 N/2]); % Shift so DFT is centered
imwrite(S,'sfimg2ph.bmp'); % Save DFT image
% Filter the vertical signal
H=ones(N,N); % Initialize freq. domain filter
H(1:N,1)=0; % Clear the vertical line
H(1,1)=1; % Preserve DC value
F=H.*F; % Apply filter
M=log(eps+abs(F))/dB_Scale;
S=circshift(M,[N/2 N/2]); % Shift so DFT is centered
imwrite(S,'sfimg3ft.bmp'); % Write the magnitude file
Y=ifft2(F); % Compute inverse 2-D FFT
imwrite(Y,'sfimg3.bmp'); % Write the real space file
% Amplify contrast
C1=min(min(Y)); % Get lowest value
C2=max(max(Y)); % Get highest value
Y=(Y-C1)/(C2-C1); % Rescale linearly
imwrite(Y,'sfimg4.bmp'); % Write image
% Lowpass filters
LPY=filter2([1 2 1; 2 4 2; 1 2 1]/16,Y);
imwrite(LPY,'sfimg5.bmp'); % Write first low-pass filtered image
LPY=filter2(h,Y);
imwrite(LPY,'sfimg6.bmp'); % Write second image
% Fourier transform
F=fft2(LPY);
M=log(eps+abs(F))/dB_Scale;
S=circshift(M,[N/2 N/2]); % Shift so DFT is centered
S=max(cat(3,zeros(N,N),S),[],3);
imwrite(S,'sfimg6ft.bmp');
% Bilateral filter
MS = 2.5; % Spatial distance parameter for Gaussian distribution
RS = ceil(MS*3); % Use bounding radius where Gaussian is small (ie. exp(-9))
MR = mean(mean(abs(Y-mean(mean(Y))))); % Get mean radiosity deviation
MR = 1.1*MR; % Use radiosity parameter of 1.1 mean deviation
BLY = zeros(N,N); % Bilateral sum (accumulated)
BLW = zeros(N,N); % Bilateral weights (accumulated)
YPad = [ zeros(RS,N+2*RS);
[zeros(N,RS) Y zeros(N,RS)];
zeros(RS,N+2*RS)]; % Build zero padded matrix
for i1 = -RS:RS % Loop through usable shift vectors
for i2= -RS:RS
WS = i1^2+i2^2; % Compute Euclidean distance
if WS < RS^2 % Use rotationally-symmetric Gaussian
WS = exp(-WS/(MS^2)); % Compute spatial weighting scalar
S = YPad( (1+RS+i1):(N+RS+i1), (1+RS+i2):(N+RS+i2) ); % Shift the Y matrix
WR = exp(-((Y-S).^2)/(MR^2)); % Compute radiosity weighting matrix
WR = WS*WR; % Combine weights
BLY = BLY + (WR.*S) ; % Accumulate bilateral summation
BLW = BLW + WR ; % Accumulate bilateral weighting
end
end
end
BLY=BLY./BLW; % Normalize by weighting applied
BLY=max(cat(3,zeros(N,N),BLY),[],3); % Limit to values >=0
imwrite(BLY,'sfimg7.bmp'); % Write image file
%Fourier transform
F=fft2(BLY); % Compute 2-D FFT
M=log(eps+abs(F))/dB_Scale; % Get magnitude spectrum
S=circshift(M,[N/2 N/2]); % Shift so DFT is centered
S=max(cat(3,zeros(N,N),S),[],3); % Limit to values >=0
S=min(cat(3,ones(N,N),S),[],3); % Limit to values <=1
imwrite(S,'sfimg7ft.bmp'); % Write image file
image(cat(3,S,S,S)); % Display final image
Experiment 3 - See experiment page for description
% EXPERIMENT 3: UNBLURRING IMAGE - CREATE TEST IMAGE 11/23/04
clear, close all
dB_Scale=11;
N=512; % Define processing unit size
h=[1 2 4 2 1 ;
2 4 8 4 2 ;
4 8 16 8 4 ;
2 4 8 4 2 ;
1 2 4 2 1 ]; % Low-pass filter
h=h/sum(sum(h));
A=imread('windgen.jpg');
A=A(:,(1+end/5):end,:);
S=min([size(A,1) size(A,2)]); % Compute size of image
S=floor(S/N); % Compute scale factor
X=(sum(A,3,'double')/3)/256; % Compute grayscale
S2=floor(S/((size(h,1)-1)/2)); % Filter iterations needed
F=[1];
for iteration=1:S2 % Filter as many times as indicated
F=conv2(F,h); % Low-pass filter to reduce aliasing
end
X=filter2(F,X);
X=X(1:S:end,1:S:end); % Decimate in 2-D
S1=floor((size(X,1)-N)/2); % Compute crop area
S2=floor((size(X,2)-N)/2);
X=X((1+S1):(S1+N),(1+S2):(S2+N),:); % Crop gray image
C1=min(min(X));
C2=max(max(X));
X=0.5 + 0.5* (X-C1)/(C2-C1); % Balance image
imwrite(X,'ubimg1.bmp'); % Write image original
Y=imread('ubimg1.bmp'); % Read quantized original
b=zeros(39,39); % Build blurring filter
b(20-3,20-4)=.7;
b(20+3,20+4)=.3;
Y=filter2(b,Y); % Apply filter
C1=min(min(Y((1+5):(end-5),(1+5):(end-5))));
C2=max(max(Y((1+5):(end-5),(1+5):(end-5))));
Y=0.5 + 0.5* (Y-C1)/(C2-C1); % Linear rebalance
Y=max(cat(3,zeros(N,N),Y),[],3); % Clip values < 0
Y=min(cat(3,ones(N,N),Y),[],3); % Clip values > 1
imwrite(Y,'ubimg2.bmp'); % Write blurred image
image(cat(3,Y,Y,Y)); % Show test input image
% EXPERIMENT 3: UNBLURRING IMAGE - PROCESS TEST IMAGE 11/23/04
% Initialization
clear, close all
h=[1 2 4 2 1 ;
2 4 8 4 2 ;
4 8 16 8 4 ;
2 4 8 4 2 ;
1 2 4 2 1 ]; % Low-pass filter
h=h/sum(sum(h)); % Normalize the filter
X=double(imread('ubimg2.bmp'))/256; % Read double-precision image
N=min(size(X)); % Get base image size (ought to be a square image)
% Create and apply a Gaussian window
W=zeros(N,N); % Reset window
D=[1:N]; % Create ascending sequence
D=(D-(N/2+1)).^2; % Compute squared distance from center
W(N/2+1,:)=D; % Create central row of square distances
for iteration=1:(N/2) % Loop through remaining rows
W(iteration,:)=D + (iteration-(N/2+1))^2; % Use Euclidean metric (sum of squares)
if iteration>1
W(N-(iteration-2),:)=W(iteration,:); % Matrix has vertical symmetry
end
end
W=exp(-W/((N/2)^2)); % Create Gaussian window using squared distances
imwrite(W,'ubimg3.bmp'); % Write window image
Y=X.*W; % Apply window in spatial domain
imwrite(Y,'ubimg4.bmp'); % Write windowed image
FX=fft2(Y); % Compute Fourier transform
S=log(abs(FX)); % Get the log of the magnitude
S=(S-min(min(S)))/(max(max(S))-min(min(S))); % Linear normalization to [0,1]
S=circshift(S,[N/2 N/2]); % Shift so DC value is centered
imwrite(S,'ubimg4ft.bmp'); % Write the magnitude spectrum image
% Compute 2-D real cepstrum
CX=real(ifft2(log(abs(FX)))); % Compute real cepstrum
CX=(CX-min(min(CX)))/(max(max(CX))-min(min(CX))); % Rescale to [0,1]
S=circshift(CX,[N/2 N/2]); % Shift so spectral-DC value is centered
imwrite(S,'ubimg4ct.bmp'); % Write the cepstrum image
% Attempt to locate echo and build filter
EX=CX.*W; % Apply Gaussian window to cepstral data
EX(1,1)=max(max(EX)); % Get maximum cepstrum value
VF=filter2([0 1 0; 0 2 0; 0 1 0]/4,EX);
HF=filter2([0 0 0; 1 2 1; 0 0 0]/4,EX);
EX=EX-max(cat(3,VF,HF),[],3); % Subtract from directional mean values
EX=max(cat(3,EX,zeros(N,N)),[],3); % Extract positive values only
EX=(EX-min(min(EX)))/(max(max(EX))-min(min(EX))); % Rescale
Cutoff=30*mean(mean(EX)); % Define cutoff
S=zeros(N,N);
for S1=1:30
for S2=1:30
if EX(S1,S2)>Cutoff
S(S1,S2)=EX(S1,S2); % Build data points where ceptrum exceeds cutoff
end
end
end
EX=S.^0.4; % Amplitude adjustment (fudge factor)
EX(2:end,2:end)=fliplr(flipud(EX(2:end,2:end))); % Flip the directional sense
S=circshift(EX,[N/2 N/2]); % Shift to center of image
imwrite(S,'ubimg5.bmp'); % Write the computed filter points
EX=EX/sum(sum(EX));
FX=fft2(EX); % Compute filter response
S=(FX-min(min(FX)))/(max(max(FX))-min(min(FX)));
S=circshift(FX,[N/2 N/2]); % Shift to center of image
imwrite(S,'ubimg5ft.bmp'); % Write the filter's magnitude response
% Reverse filter using deconvolution in FFT domain
FX=fft2(X); % Compute image FFT
S=fft2(EX); % Compute estimated filter's FFT
%(in this case, the denominator stays away from zero)
FY=FX./S; % Complex division (deconvolution in spatial domain)
Y=ifft2(FY); % Inverse Fourier transform
imwrite(Y,'ubimg6.bmp'); % Write the final image file
Experiment 4 - See experiment page for description
% EXPERIMENT 4: DECONVOLUTION MICROSCOPY - CREATE TEST IMAGES
clear, close all
N=300; % Define processing unit size
db_Scale=10; % Arbitrary log scale
Images=7; % Number of images
ArrayXYZ=zeros(N,N,Images); % Reset 3D array
CrossXZ=zeros(Images,N); % Cross section view
CrossZY=zeros(N,Images); % Cross section view
h=[1 2 4 2 1 ;
2 4 8 4 2 ;
4 8 16 8 4 ;
2 4 8 4 2 ;
1 2 4 2 1 ]; % Low-pass filter
h=h/sum(sum(h)); % Normalize
% Gather images
for img=1:Images % Loop through numbered images
A=imread(['onion',int2str(img),'.jpg']); % Read image file
S=min([size(A,1) size(A,2)]); % Compute size of image
S=floor(S/N); % Compute scale factor
X=double(sum(A,3)/3)/256; % Compute grayscale
S2=floor(S/((size(h,1)-1)/2)); % Filter iterations needed
F=[1];
for iteration=1:S2 % Filter as many times as indicated
F=conv2(F,h); % Low-pass filter to reduce aliasing
end
X=filter2(F,X);
X=X(1:S:end,1:S:end); % Decimate in 2-D
S1=floor((size(X,1)-N)/2); % Compute crop area
S2=floor((size(X,2)-N)/2);
X=X((1+S1):(S1+N),(1+S2):(S2+N),:); % Crop gray image
C1=mean(mean(X)); % Get mean value
C2=mean(mean(abs(X-C1))); % Get variance
X=0.5+(0.5/(4*C2))*(X-C1); % Rescale linearly
X=max(cat(3,zeros(N,N),X),[],3); % Clip values < 0
X=min(cat(3,ones(N,N),X),[],3); % Clip values > 1
ArrayXYZ(:,:,img)=X; % Store in appropriate array layer
CrossXZ(img,:)=X(N/2,:); % Build cross section XZ
CrossZY(:,img)=X(:,N/2); % Build cross section ZY
imwrite(X,['dcimg',int2str(img),'.bmp']); % Write gray images
end
imwrite(CrossXZ,'dcxz.bmp'); % Write cross section image
imwrite(CrossZY,'dczy.bmp'); % Write cross section image
X=mean(ArrayXYZ,3);
imwrite(X,'dcmean.bmp'); % Write image of mean values
%image(cat(3,X,X,X));
% Compute Fourier transform
F=fft2(X);
M=log(eps+abs(F))/db_Scale;
S1=floor(N/2); % Recalculate bounds for FFT graph
S2=S1+1;
S=[M(S2:N,S2:N) M(S2:N,1:S1); M(1:S1,S2:N) M(1:S1,1:S1)];
S=max(cat(3,zeros(N,N),S),[],3);
S=min(cat(3,ones(N,N),S),[],3);
imwrite(S,'dcmeanft.bmp');
% EXPERIMENT 4: DECONVOLUTION MICROSCOPY - PROCESS IMAGES
clear, close all
dB_Scale=11;
N=300; % Define processing unit size
Images=7; % Number of images
ArrayXYZ=zeros(N,N,Images); % Reset 3D array
psf=zeros(N,N,Images); % Reset point-spread function
X=zeros(N,N);
X(1,:)=(min([0:(N-1); N:-1:1],[],1)).^2;
for iteration=2:N
X(iteration,:)=X(1,:)+X(1,iteration); % create Euclidean squared metric
end
for img=1:ceil(Images/2) % Loop through psf layers
S=exp(-X/((0.5+4*(img-1))^2)); % Compute Gaussian distribution
S=S/sum(sum(S)); % Normalize on this layer
S=S*exp(-(img-1)^2/(1.0^2)); % Fade with distance from central layer
psf(:,:,img)=S;
if img>1
psf(:,:,Images-(img-2))=S;
end
end
psf=psf/sum(sum(sum(psf))); % Normalize 3-D psf array
X=zeros(N,N); % Reset psf visualization array
for img=1:Images % Create psf graphic, superimposing each layer
S=psf(:,:,img);
C1=floor(Images/2);
C1=(mod(img-1+C1,Images)-C1)*floor(N/Images);
C2=sum(sum(S))/40;
X=X+circshift(S/C2,[N/2+C1 N/2]);
end
X=min(cat(3,ones(N,N),X),[],3);
C1=floor(N/Images);
imwrite(X(:,(end/2-C1):(end/2+C1)),'dcpsf.bmp');
% Gather images
for img=1:Images % Loop through numbered images
X=imread(['dcimg',int2str(img),'.bmp']); % Write gray images
X=double(X)/256;
ArrayXYZ(:,:,img)=X;
end
% Show frequency response of psf
X=fft2(psf(:,:,2));
M=log(eps+abs(X(:,:,1)))/dB_Scale;
M=(M-min(min(M)))/(max(max(M))-min(min(M)));
S=circshift(M,[N/2 N/2]);
imwrite(S,'dcpsfft.bmp');
% Perform Deconvolution in Frequency Domain
X=fftn(psf);
X=max(cat(4,0.0001+zeros(N,N,Images),X),[],4);
F=fftn(ArrayXYZ);
F=F./X;
ArrayXYZ=ifftn(F);
% Write the deconvolved images
for img=1:Images % Loop through numbered images
X=real(ArrayXYZ(:,:,img));
C1=mean(mean(X)); % Get mean value
C2=mean(mean(abs(X-C1))); % Get variance
X=0.5+(0.5/(4*C2))*(X-C1); % Rescale linearly
X=max(cat(3,zeros(N,N),X),[],3);
X=min(cat(3,ones(N,N),X),[],3);
imwrite(X,['dcres',int2str(img),'.bmp']); % Write gray images
end
Experiment 5 - Extension of smoothing filter to color photos
% EXPERIMENT 5: NOISE FILTER - CREATE TEST FILES
% Initializations
clear, close all
N=320; % Define processing unit size
% Read main image
A=imread('land.jpg'); % Read image file
S=min([size(A,1) size(A,2)]); % Compute size of image
S=floor(S/N); % Compute scale factor
A=A(1:S:end,1:S:end,:); % Decimate in 2-D
S1=floor((size(A,1)-N)/2); % Compute crop area
S2=floor((size(A,2)-N)/2);
A=A((1+S1):(S1+N),(1+S2):(S2+N),:); % Crop image
imwrite(A,'nfimg1.bmp'); % Write image
%S=(1-sign(rand(N,N)-0.1))/2;
%S=2*S.*(rand(N,N)-0.5);
%B=cat(3,S,S,S);
B= 0.2 * 2 * (rand(N,N,3)-.5); % Build noise array
A= uint8(double(A)+256*B);
imwrite(A,'nfimg2.bmp'); % Write noisy image
% EXPERIMENT 5: NOISE FILTER - PROCESS COLOR BILATERAL FILTER
% Read main image
clear, close all
NN=double(imread('nfimg1.bmp'))/256; % Read input image
A=double(imread('nfimg2.bmp'))/256; % Read input image
N=size(A,1); % Get image size
% Filter by linear low-pass filter
h=[ 1 2 4 2 1;
2 4 8 4 2;
4 8 16 8 4;
2 4 8 4 2;
1 2 4 2 1 ];
h=h/sum(sum(h));
B=zeros(N,N,3);
for color=1:3
B(:,:,color)=filter2(h,A(:,:,color));
end
imwrite(B,'nfimg3.bmp');
% Filter using bilateral filter
MS=3.9; % Mean spatial variation of filter
MR=0.37; % Mean radiosity variation of filter
MS=3+rand(1);
MR=0.3+0.1*rand(1);
R=ceil(3*MS); % Radius of bilateral action
WT=zeros(N,N); % Reset weights
B=zeros(N,N,3); % Reset output image
for C1=-R:R
for C2=-R:R
DS=C1^2+C2^2; % Compute distance squared
if DS < R^2 % If within significant proximity...
S=circshift(A,[C1 C2 0]); % Compute the shifted matrix and weights
W = exp(-DS/(MS^2)) * exp(-sum((S-A).^2,3)/(MR^2));
WT = WT + W; % Total all the weights
B = B + S.*cat(3,W,W,W); % Compute the sum
end
end
end
for color=1:3
B(:,:,color)=B(:,:,color)./WT;
end
imwrite(B,'nfimg4.bmp');
