function output=LAPLACE3Dsmooth(input, source, sink, n_smooth)
%output=LAPLACE3Dsmooth(input, source, sink, n_smooth)
%
% The function performs the laplace smoothing in 3D mesh.
% 
% input       : input signal to be smoothed.
% heat source + sink : heat sink for heat equation. This is equivalent to the boundary 
%                    condition for the laplace equation. 
%                WARNING: it is assumed input ~=heatsink is the index
%                corresponding to heat sinks.
%
% sigma      : bandwidth. The bandwidth corresponds to time in a heat equation [3]. 
%                   The bandwidth given in old publications [1] [2] are related to 
%                    new publication [3] as our new_sigma = old_sigma^2/2.
% n_smooth  : number of iterations. 5 iterantions are sufficient for
%                     amygdale. In [4], 30 iterations were used but you
%                     don't need that many iterations. 
%
% EXAMPLE: bandwidth sigma=1, number of iteration=5
% output_signal=hk_smooth(input_signal,tri,coord,nbr,1,5);
%
% (C) Moo K. Chung, mkchung@wisc.edu
% 
% Update history Feb 5, 2004; July 20, 2004; Oct 28, 2004; Sept 18, 2007; Nov 13, 2007.
% http://www.stat.wisc.edu/softwares/hk/hk.html
%
% If you use this code, please reference one of the following papers. The details on
% the mathematical basis of of the algorithm can be found in these papers.
%
% [1] Chung, M.K., Robbins,S., Dalton, K.M., Davidson, Alexander, A.L., R.J., Evans, A.C. 2005. 
%     Cortical thickness analysis in autism via heat kernel smoothing. NeuroImage 25:1256-1265
%     http://www.stat.wisc.edu/~mchung/papers/ni_heatkernel.pdf
%
% [2] Chung, M.K., Robbins, S., Evans, A.C. 2005. Unified statistical approach to cortical 
%     thickness analysis. Information Processing in Medical Imaging (IPMI). Lecture Notes in 
%     Computer Science (LNCS) 3565:627-638. Springer-Verlag
%     http://www.stat.wisc.edu/~mchung/papers/IPMI/hk.IPMI.2005.pdf
%
% [3] Chung, M.K. Hartley, R., Dalton, K.M., Davidson, R.J. 2007. Encoding cortical surface 
%     by spherical harmonics. Special Issue of Statistica Sinica on Statistical Challenges 
%     and Advances in Brain Science, in press
%     http://www.stat.wisc.edu/~mchung/papers/BMI.198.pdf
%
% [4] Chung, M.K., Worsley, K.J., Nacewicz, B.M., Dalton, K.M., Davidson, R.J.
%       2010. General multivariate linear modeling of surface shapes using SurfStat. 
%      NeuroImage. 53:491-505. 
%
%
%

%discrete heat kernel
% For publications [1] and [2], change this line to  
% K=inline('exp(-x/(2*sigma^2))/sum(exp(-x/(2*sigma^2)))');

% debugging line
%input=amygcube; source=amyg;sink=-cube;

d=size(input);
source=double(source);
sink=double(sink);
input=double(input);

step=5;sigma=10;
kernel=zeros(2*step+1,2*step+1,2*step+1);
for i=-step:step
    for j=-step:step
        for k=-step:step
            kernel(step+1+i,step+1+j,step+1+k)= exp(-(i^2+j^2+k^2)/sigma);
        end;
    end;
end;
kernel=kernel/sum(sum(sum(kernel)));

for i=1:n_smooth  
    
    output = convn(input,kernel, 'same');

    %figure;imagesc(output(:,:,rout));colorbar

    %figure;imagesc(input(:,:,rout));colorbar
    %figure;imagesc(sink(:,:,rout));colorbar

   
     A=output.*(sink);
     %figure;imagesc(A(:,:,rout));colorbar
     
    temp=output+A;
    %figure;imagesc(temp(:,:,rout));colorbar
    
    C=temp.*source;
    %figure;imagesc(C(:,:,rout));colorbar
    
    temp=temp-temp.*source + source +sink;
    
    %figure;imagesc(temp(:,:,rout));colorbar


    input=temp;

end;

temp=output+output.*sink;
temp=temp-temp.*source + source +sink;

output=temp;



  
