% questo script non implementa nessuna contromisura per le
% periodicità della FFT o altri fenomeni indesiderati
%
function image_filt1
	
	% leggi da un file l'immagine
	a = double(imread('imm1.png'));

	% se l'immegine non è in scala di grigio, calcola la componente
	% di luminosità
	if size(a,3) > 1
		a = 0.299*a(:,:,1)+0.114*a(:,:,2)+0.587*a(:,:,3);
	end

	% plot dell'immagine
	figure(1)
	subplot(2,1,1)
	imagesc(a,[0,255])
	axis image
	colormap gray

	% vettori delle frequenze (per mettere le etichette corrette
	% agli assi)
	[nr,nc]=size(a);
	u = -nr/2:nr/2-1;
	v = -nc/2:nc/2-1;

	% calcola la DFT
	A=fft2(a);
	A(1,1) = (2 * (A(2,1) + A(1,2) + A(end,1) + A(1,end)) ...
		-(A(2,2) + A(end,2) + A(2,end) + A(end,end)))/4;

	% plot della DFT (con compressione logaritmica)
	subplot(2,1,2)
	imagesc(u,v,log(1+abs(fftshift(A))))
	axis image
	colormap gray

	% deviazione standard del rumore aggiunto
	noise_sigma = 30;

	% immagine rumorosa
	a_noisy = a + noise_sigma*randn(size(a));
	% saturazione
	a_noisy = min(max(a_noisy,0),255);

	% plot dell'immagine rumorosa
	figure(2)
	subplot(2,1,1)
	imagesc(a_noisy,[0,255])
	axis image
	colormap gray

	% DFT dell'immagine rumorosa
	A_noisy=fft2(a_noisy);
	A_noisy(1,1) = (2 * (A_noisy(2,1) + A_noisy(1,2) + A_noisy(end,1) + A_noisy(1,end)) ...
		-(A_noisy(2,2) + A_noisy(end,2) + A_noisy(2,end) + A_noisy(end,end)))/4;

	% plot della DFT rumorosa
	subplot(2,1,2)
	imagesc(u,v,log(1+abs(fftshift(A_noisy))))
	axis image
	colormap gray

	% loop su alcune scelte possibili per la varianza della
	% risposta in frequenza del filtro Gaussiano
	for sigma2 = 10.^(5:-1/3:1)

		A_denoise = A_noisy;

		% filtro passa basso "naif", perché non funziona?
		% A_denoise(1+(19:nc-19),:) = 0;
		% A_denoise(:,1+(4:nr-4)) = 0;

		% filtro Gaussiano
		[uH,vH] = meshgrid(fftshift(u),fftshift(v));
		H = exp(-(uH.^2+vH.^2)/sigma2);
		A_denoise = A_denoise .* H;

		% visualizza la DFT dell'immagine filtrata
		figure(3)
		subplot(2,1,2)
		imagesc(u,v,log(1+abs(fftshift(A_denoise))))
		axis image
		colormap gray

		% calcola la IDFT e rinormalizzala
		a_denoise = real(ifft2(A_denoise));
		a_denoise = (a_denoise - min(a_denoise(:))) / ...
			(max(a_denoise(:))-min(a_denoise(:)))*255;
		
		% visualizza l'immagine filtrata
		subplot(2,1,1)
		imagesc(a_denoise)
		axis image
		colormap gray

		% calcola l'errore (visto che abbiamo l'immagine
		% originale si può fare)
		MSE = sum(abs(a(:) - a_denoise(:)).^2)/numel(a);
		fprintf('sigma2: %6.0f MSE: %7.2f\n',sigma2,MSE);
		
		drawnow
	end
end
