% Statistical and Computational Methods in Brain Image Analysis
%Chapter 4 Gaussian Kernel Smoothing
%
% (C) 2013 Moo K. Chung 
% University of Wisconsin-Madison
% mkchung@wisc.edu

%-----------------------------------------------------
% GAUSSIAN KERNEL SMOOTHING IN 1D

% truncated Gaussian Kernel  in 1D
K=inline('exp(-(x.^2)/2/sigma^2)');
dx= -5:5;

sum(K(10,dx))
weight=K(10,dx)/sum(K(10,dx))
sum(weight)

t=1:100;
mu=(t-50).^2/500;
noise= normrnd(0, 2, 1,100);
Y=mu + noise
figure; plot(Y, '.');

smooth=conv(Y,weight,'same');
hold on; plot(smooth, 'r'); 


figure_bg('w')
print_pdf('w2')


%-----------------------------------
% GAUSSIAN KERNEL SMOOTHING IN 2D

K=inline('exp(-(x.^2+y.^2)/2/sigma^2)'); 
[dx,dy]=meshgrid([-10:10]);  

sigma=100;
weight=K(sigma,dx,dy)/sum(sum(K(sigma,dx,dy)));

% toy example
signal= imread('toy-key.tif');
signal=(double(signal)-219)/36; % makes the image into real number
figure; subplot(2,2,1); imagesc(signal); colormap('bone'); colorbar %596 x 368
title('Signal')

% Gaussian white noise
noise= normrnd(0, 5, 596, 368);
f = signal + noise;
subplot(2,2,2); imagesc(f); colormap('bone'); colorbar
title('Signal + N(0,5)')

weight=K(1,dx,dy)/sum(sum(K(1,dx,dy)));
smooth=conv2(f,weight,'same');
subplot(2,2,3); imagesc(smooth); colormap('bone'); colorbar
title('Kernel smoothing with \sigma=1')

weight=K(10,dx,dy)/sum(sum(K(10,dx,dy)));
smooth=conv2(f,weight,'same');
subplot(2,2,4); imagesc(smooth); colormap('bone'); colorbar
title('Kernel smoothing with \sigma=10')

%--------------------------------
% Increased Gaussianness
for i=1:50
    noise= normrnd(0, 5, 596, 368);
    weight=K(100,dx,dy)/sum(sum(K(100,dx,dy)));
    f=signal+noise;
    smooth=conv2(f,weight,'same');
    pixelvalue(i)=f(314, 150);
    pixelvalues(i)=smooth(314, 150);
end;

figure; subplot(1,2,1); qqplot(pixelvalue)
subplot(1,2,2); qqplot(pixelvalues)



