stacksize(100000000);
 
// __________________________________________________
// __________________________________________________
 
// erreur quadratique moyenne
function [err] = EQM( img_new, img_ref )
  err = mean( abs( img_new - img_ref ).^2 );
endfunction
 
 
function [result] = distance( x1,y1, x2,y2 )
	result = sqrt( (y2-y1)^2 + (x2-x1)^2 );
endfunction
 
// filtre passe-bas Gaussien
function [result] = GaussianFilter_in_FrequencySpace( img_width, img_height, cutoff_freq )
 
	center_x = img_width  / 2.0;
	center_y = img_height / 2.0;
 
	result = zeros( img_width, img_height );
 
	temp_d = 2.0 * cutoff_freq^2;
 
	for y=1:img_height
		for x=1:img_width
		result(x,y) = exp ( - (distance( x,y,center_x,center_y ))^2 / temp_d );
		end
	end
 
endfunction
 
// filtre passe-bas ideal
function [result] = removeFreq_in_FrequencySpace( img_width, img_height, seuil );
 
	center_x = img_width  / 2.0;
	center_y = img_height / 2.0;
 
	result = zeros( img_width, img_height );
 
	for y=1:img_height
		for x=1:img_width
			if ( distance( x,y,center_x,center_y ) < seuil ) then
				result(x,y) = 1;
			end
		end
	end
 
endfunction
 
// __________________________________________________
// __________________________________________________
 
 
// __________________________________________________
// chargement de l'image
img = imread ("texture.pgm");
[img_size]  = size(img);
img_width   = img_size(1);
img_height  = img_size(2);
img_nbPixel = img_width * img_height;
//imshow (img);
 
 
// __________________________________________________
// on applique la fft
img_fft = fft( img, -1 ) / sqrt( img_width*img_height );
 
// __________________________________________________
// on affiche l'image dans le dommaine des frequences
img_fft_log = log( fftshift( abs(img_fft) ) + %eps );
imwrite("imgOUT_fft_log.pgm", img_fft_log, "PGM");
 
// __________________________________________________
// on cree le filtre : soit Gaussien soit Ideal
 
img_filter = GaussianFilter_in_FrequencySpace( img_width, img_height, 50 );
 
imwrite("imgOUT_filter.pgm", img_filter*255, "PGM");
 
// __________________________________________________
// on filtre l'image dans le dommaine frequenciel
img_result_fft = fftshift( img_filter ) .* img_fft;
 
// __________________________________________________
// visualisation logaritmique de l'image filtree dans le dommaine frequentiel
img_filtered_fft_log = log( fftshift( abs(img_result_fft) ) + %eps );
imwrite("imgOUT_filtered_fft_log.pgm", img_filtered_fft_log, "PGM");
 
// __________________________________________________
// fft inverse...
img_ifft = fft( img_result_fft, 1 ) * sqrt( img_width*img_height );
// imshow( real( img_ifft ) );
imwrite("imgOUT_result.pgm", real(img_ifft), "PGM");