function surf=LAPLACEcontour(stream,surf, level)
% Moo K. Chung 2009
% University of Wisconsin-Madison
% deforms amygdala surface following level set of the stream map.

d=size(stream);
temp=zeros(d(1)*d(2)*d(3),1);

temp2=stream;

temp2=reshape(temp2,d(1)*d(2)*d(3),1);

temp(find(temp2>level))=1;

% the following line is for sink to source backward propagation
%temp(find(temp2<level))=1;

temp=reshape(temp,d(1),d(2),d(3));

    
tempsurf = isosurface(temp,'noshare');
%tempsurf = isosurface(temp);

%figure_wire(tempsurf,'yellow')

n_vertex=size(surf.vertices,1); %# of vertex in amydala
m_vertex=size(tempsurf.vertices,1);   % # of vertex in contour

for i=1:n_vertex
    % for each point in amydala, we compute the distance to contour
    dist=tempsurf.vertices-kron(ones(m_vertex,1),surf.vertices(i,:));
    dist=L2norm(dist);
    
    ind=find((dist -min(dist))<=3);  % this controls the amount of smoothing
    %ind=find(dist==min(dist));
    surf.vertices(i,:)=mean(tempsurf.vertices(ind,:),1);
    
    %surf.vertices(i,:)=tempsurf.vertices(ind(1),:);

end;

%figure_wire(surf,'yellow')


function L2=L2norm(vector)
% It computes the L2-norm of n by 3 matrix and outputs n by 1 array of L2-norm.

L2=sqrt(sum(vector.^2,2));


