function output=laplace_smooth(input, heatsink, tri,coord,nbr,sigma,n_smooth)
%output=laplace_smooth(input,tri,coord,nbr,sigma,n_smooth)
%
% The function performs the laplace smoothing [4] on a triangle mesh.
% 
% input       : input signal to be smoothed.
% heatsink : 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.
% tri            : triangle indices
% coord       : vertex coordinates
% nbr          : first neighbor vertex list
% 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
%
% EXAMPLE: bandwidth sigma=1, number of iteration=50
% output_signal=hk_smooth(input_signal,tri,coord,nbr,1,50);
%
% (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] future publication.
%
%
%

%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)))');

K=inline('exp(-x/(4*sigma))/sum(exp(-x/(4*sigma)))');


%kernel weight computation

n_vertex=size(nbr,1);
max_degree=max(size(nbr,2));
weight=zeros(n_vertex,max_degree+1); %weight including the current vertex

for i=1:n_vertex
    current_nbr=nbr(i,nbr(i,:)~=0);
    degree=length(current_nbr);
    pipo=[coord(current_nbr,:)-kron(ones(degree,1),coord(i,:))]';
    distance=sum(pipo.*pipo);
    weight(i,1:(1+degree))=K(sigma,[0 distance]);
end;


% iterative weighted average

% increased the index by one to account for 0 indexing in nbr.
% this is necessary to weighted averaging in matrix form to speed up
% computation.
% signal(1) is the dummy entry corresponding to 0 indexing.
signal=zeros(n_vertex+1,1);
signal(2:n_vertex+1)=input;

size(input)

% find heat sinks. it is assumed input ~=heatsink is the index corresponding to heat sinks.
ind=[];
for i=1:n_vertex
    if heatsink(i)~=input(i)
        ind=[ind i];
    end;
end;
ind=find(heatsink~=input);
heatsink=heatsink(ind);

dummy=ones(1,max_degree);
for j=1:n_smooth    
    signal(ind+1)= heatsink;  %assign heat sinks
    Y=[signal signal([dummy; nbr+1])];  %1st row is a dummy entry
    signal(2:n_vertex+1)= sum(weight.*Y(2:n_vertex+1,:),2); 
end;
signal(ind+1)=heatsink;
output =signal(2:n_vertex+1);





%  OLD CODE for MNI mesh format which has the maximum degree of 6.
%  excluding the center, there are maximum total of 6 neighbors.
%
% for i=1:n_vertex
%      the first 12 vertices have only 5 neighbors
%     if (i<=12)
%         current_nbr=nbr(i,1:5);
%         pipo=coord(:,current_nbr(1:5))-kron(ones(1,5),coord(:,i));
%          distance=sum(pipo.*pipo);
%          weight(i,:)=[K(sigma,[0 distance]) 0];
%      else
%          % the remaining vertices have 6 neighbors
%          current_nbr=nbr(i,1:6);
%         pipo=coord(:,current_nbr(1:6))-kron(ones(1,6),coord(:,i));
%          distance=sum(pipo.*pipo);
%          weight(i,:)=K(sigma,[0 distance]);
%      end;
%  end;
%iterated smoothing
% signal=input_signal;
% for j=1:n_smooth
%     Y=[signal [[signal(nbr(1:12,1:5)) zeros(12,1)]; signal(nbr(13:n_vertex,:))]]; % signal + neighboring signal
%     signal=sum(weight.*Y,2);  %weighed averaging
% end;


