% imfft.m
% 
% [filtered_photo,offt,cutout]=
%                    imfft(photo,type,centerX,centerY,radius,quad,mirror)
%
% hipass-filters a given photo with a rectangular 'R' or circular 'C'
% mask type in the frequency domain 
% When mask is circular 'C', use radius and center point (X,Y).
% When mask is rectangular 'R', use center and radius as the half diagonal
%
% Use quad values to specify quadrant for the ifft cutout mask.
%     0: defined by centerX, centerY
%     1: defined by -centerX,centerY
% Use mirror to mirror quadrants
%
% Also returns fft transformed image offt
%      (can be viewed with imagesc(log(abs(fftshift(offt))+1)) )
% and the cutout image
%      (can be viewed with imagesc(log(abs(fftshift(cutout))+1)) )
%
% BASED ON fft_filter --- 1999-02-28:DAV: David Cary <d.cary@ieee.org> started.
%

%
% 2003 May 15 ©gpavlid -- if 'type' is an image it is the mask
%
function [filtered_photo,offt,cutout] = imfft( photo, type, varargin);
% centerX, centerY, radius, quad, mirror
[height,width] = size(photo);

% do the fft on original image
fprintf('fft transforming image...\n');
offt = fft2(photo);

if min(size(type))>1,
	% type is a gray-level mask image
	mask = type;
	cutout = fftshift(mask) .* offt;
	type = 'mask image';
	centerX = 0;
	centerY = 0;
	radius = 0;
	quad = 0;
	mirror = 0;
	show = 0;
else
	centerX = varargin{1};
	centerY = varargin{2};
	radius = varargin{3};
	quad = varargin{4};
	mirror = varargin{5};
	if length(varargin)>5
		show = varargin{6};
	else
		show = 0;
	end;
% compose the filter
fprintf('composing cutout mask ...\n');
cX=centerX;cY=centerY;
if upper(type)=='C',
	mask_part = circle(2*radius+1);
else
	mask_part = ones(2*radius+1);
end;
mask = zeros(size(photo));
y = -radius:radius;
x = -radius:radius;
y=y+cY;
x=x+cX;
mask(floor(height/2) + y, floor(width/2) + x) = mask_part;
mask = 1 - mask;
cutout = fftshift(mask) .* offt;
if mirror,
	fprintf('composing the mirror cutout mask...\n');
	cX=-centerX;cY=-centerY;
	if upper(type)=='C',
		mask_part = circle(2*radius+1);
	else
		mask_part = ones(2*radius+1);
	end;
	mask = zeros(size(photo));
	y = -radius:radius;
	x = -radius:radius;
	y=y+cY;
	x=x+cX;
	mask(floor(height/2) + y, floor(width/2) + x) = mask_part;
	mask = 1 - mask;
	cutout = fftshift(mask) .* cutout;
end;
if quad,
	fprintf('composing additional quadrant cutout mask...\n');
	if centerX==0 | centerY==0,
		cX=-centerY;cY=centerX;
	else
		cX=-centerX;cY=centerY;
	end;
	if upper(type)=='C',
		mask_part = circle(2*radius+1);
	else
		mask_part = ones(2*radius+1);
	end;
	mask = zeros(size(photo));
	y = -radius:radius;
	x = -radius:radius;
	y=y+cY;
	x=x+cX;
	mask(floor(height/2) + y, floor(width/2) + x) = mask_part;
	mask = 1 - mask;
	cutout = fftshift(mask) .* cutout;
if mirror,
	fprintf('composing the additional quadrant mirror cutout mask...\n');
	if centerX==0 | centerY==0,
		cX=centerY;cY=-centerX;
	else
		cX=centerX;cY=-centerY;
	end;
	if upper(type)=='C',
		mask_part = circle(2*radius+1);
	else
		mask_part = ones(2*radius+1);
	end;
	mask = zeros(size(photo));
	y = -radius:radius;
	x = -radius:radius;
	y=y+cY;
	x=x+cX;
	mask(floor(height/2) + y, floor(width/2) + x) = mask_part;
	mask = 1 - mask;
	cutout = fftshift(mask) .* cutout;
end;
end;
end;

% do the ifft with cutout filter parameters
fprintf('doing the ifft transform...\n');
filtered_photo = abs( ifft2(cutout) );
fprintf('Done.\n');

if nargout==0,
titl = sprintf( 'FFT image filtering (%s,%d,%d,%d,%d,%d)',type,centerX, centerY, radius, quad, mirror );
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);imagesc(photo);zoom on;title('original image');
%subplot(2,2,2);imagesc(log(abs(fftshift(offt))+1) );zoom on;title('fft2(original)');
%subplot(2,2,3);imagesc(log(abs(fftshift(cutout))+1) );zoom on;title('fft2(filtered)');
%subplot(2,2,4);imagesc(filtered_photo);title('filtered 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 "colormap(gray(255));mesh(log(abs(fftshift(offt))+1)*10);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 "colormap(gray(255));mesh(log(abs(fftshift(cutout))+1)*10);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');
end;

if show
	imdisp(filtered_photo,0,1,'Filtered photo');
	imdisp(log(abs(fftshift(offt))+1),0,1,'Original FFT space');
	imdisp(log(abs(fftshift(cutout))+1),0,1,'Filtered FFT space');
end;