%
% 3D DFT demonstration
% removal of periodical noise from corrupted images
%

img = double(imread( 'lena2.jpg'));											% load an image
[height,width] = size( img);								% get size of image
offt = fft2( img);													% get the 2D-DFT of the image
offt( 100, 400) = max( offt(:));						% insert a high value in magnitude of DFT
photo = abs( ifft2( offt) );								% get the noisy image back with 2D i-DFT

noiseX=144-floor(size(img,1)/2);						% fix the position of the noise frequency
noiseY=356-floor(size(img,2)/2);
noiseR=3;																		% fix the radius of the noise frequency

% consruct the circular mask
cX = noiseX;
cY = noiseY;
mask_part = circle( 2*noiseR+1);
mask = zeros( size( photo));
y = -noiseR:noiseR;
x = -noiseR:noiseR;
y = y+cY;
x = x+cX;
mask( floor( height/2) + y, floor( width/2) + x) = mask_part;
mask = 1 - mask;

cutout = fftshift(mask) .* offt;						% apply the mask over the 2D DFT
filtered_photo = abs( ifft2( cutout) );			% reconstruct the image from masked DFT

% display the results
titl = sprintf( 'FFT image filtering (%s,%d,%d,%d)','C', noiseX, noiseY, noiseR );
fig=figure( ...
 			'Name',titl, ...
			'NumberTitle','off', 'HandleVisibility', 'on', ...
			'tag', 'Fig2', ...
			'Visible','on', 'Resize', 'on',...
			'BusyAction','Queue','Interruptible','on',...
			'Menubar','none');
colormap( gray( 255)); axis image; zoom on;
subplot(2,2,1); subimage( photo/max(photo(:))); title('original image');
orig_fft_img = log( abs( fftshift( offt))+1);
subplot(2,2,2); subimage( orig_fft_img/max( orig_fft_img(:))); title('fft2(original)');
fprintf( 'Use "figure;colormap(gray(255));mesh(log(abs(fftshift(offt))+1)*10);axis tight;rotate3d on;"\n to get a 3D representation of the original FFT\n' );
filt_fft_img = log( abs( fftshift( cutout))+1);
subplot(2,2,3); subimage( filt_fft_img/max( filt_fft_img(:))); title('fft2(filtered)');
fprintf( 'Use "figure;colormap(gray(255));mesh(log(abs(fftshift(cutout))+1)*10);axis tight;rotate3d on;"\n to get a 3D representation of the cutout FFT\n' );
subplot(2,2,4); subimage( filtered_photo/max( filtered_photo(:))); title('filtered image');
