function surf=mesh_areapreserve2(surf, iter)
%
%function surf=mesh_areapreserve2(surf, iter)
%
% mesh_areapreseve performs area preserving mapping of a surface mesh onto a unit sphere.
%
% surf : surface mesh 
% iter  : number of iterations
%
% The weight matrix is coming from matrix "A" in the papers [1] and [2] or book [3] 
% where it is need to perform the finite element method (FEM). The code is downloaded from 
% http://brainimaging.waisman.wisc.edu/~chung/midus
%
%
%  Reference:
%
% [1] Chung, M.K. 2001. Statistical Morphometry in Neuroanatomy, 
%     PhD Thesis, McGill University.
%     http://www.stat.wisc.edu/~mchung/papers/thesis.pdf
%
% [2] Chung, M.K., Taylor, J. 2004. Diffusion Smoothing on Brain 
%     Surface via Finite Element Method,  IEEE International Symposium 
%     on Biomedical Imaging (ISBI). 432-435.
%
% [3] Chung, M.K. Computational Neuroanatomy: The Methods, World Scientific
%
% (C) 2013. Moo K. Chung
%     Department of Biostatistics and Medical Informatics
%     Waisman Laboratory for Brain Imaging
%     University of Wisconsin-Maison
%  
% Version 1.0. August 10, 2013 created
% Version 2.0. August 15, 2013. The difference from mesh_areapreserve.m is that A-matrix
% compuation is done wihin a loop for each updated mesh. In mesh_areapreserve.m, it is only done
% once for the initial mesh. Speedwise, mesh_areapreserve2.m is slower but provide a better mesh
% without overlapping regions.



  tri   = surf.faces;
  n=size(surf.vertices,1);
  
 
for i=1:iter
    vertex1=surf.vertices';
    vertex1 = vertex1 - repmat( mean(vertex1,2), [1 n] );
    vertex1 = vertex1 ./ repmat( sqrt(sum(vertex1.^2,1)), [3 1] );
    
    
    coord = vertex1';
    %-----------------------------------
    % A matrix computation. See reference [3] to see how it is computed
    Adj = sparse(tri(:,[1 2 3]),tri(:,[2 3 1]),1,n,n);
    Adj = double(Adj|Adj'); % adjacent nodes
    [i,j] = find(Adj);
    [ei,ej] = find(Adj(:,i) & Adj(:,j)); % common adjacent triangle-edge
    e1 = find([1; diff(ej)]); % 1:2:end  % ->ignors singularities
    e2 = e1+1;                % 2:2:end
    
    pi = coord(i,:);
    pj = coord(j,:);
    qi = coord(ei(e1),:); % qi = coord(ei(1:2:end),:);
    qj = coord(ei(e2),:); % qj = coord(ei(2:2:end),:);
    
    dii = pi-qi; dij = pi-qj;
    dji = pj-qi; djj = pj-qj;
    
    norm = @(x) sqrt(sum(x.^2,2)); % norm of a vector
    
    area = @(vi,vj) norm(cross(vi,vj))./2;
    A = sparse(i,j,(area(dii,dji) + area(dij,djj))./12,n,n);
    A = A + diag(sum(A,2));
    W=A;
    %----------------------------------------------
    
    
    d = full( sum(W,1) );
    D = spdiags(d(:), 0, n,n);
    iD = spdiags(d(:).^(-1), 0, n,n);
    tW = iD * W;
    
    
    vertex1 = vertex1*tW';
    vertex1 = vertex1 - repmat( mean(vertex1,2), [1 n] );
    vertex1 = vertex1 ./ repmat( sqrt(sum(vertex1.^2,1)), [3 1] );
    
    
    surf.vertices=vertex1';
    
end;
