function surf=mesh_areapreserve(surf, iter)
%
%function surf=mesh_areapreserve(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] 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.
% 
%
% (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



vertex=surf.vertices';
faces=surf.faces';

% Adjacency matrix
% The incidence matrix in the old code is the adjaency matrix

tri   = surf.faces;
coord = surf.vertices;
n = size(coord,1);

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 = vertex;
vertex1 = vertex1 - repmat( mean(vertex1,2), [1 n] );
vertex1 = vertex1 ./ repmat( sqrt(sum(vertex1.^2,1)), [3 1] );


% normal to faces
% [normal,normalf] = compute_normal(vertex1,faces);
% % center of faces
% C = squeeze(mean(reshape(vertex1(:,faces),[3 3 m]), 2));
% % inner product
% I = sum(C.*normalf);
% 
% options.name = 'none';
% clf;
% options.face_vertex_color = double(I(:)>0);
% figure; plot_mesh(vertex1,faces,options);
% colormap gray(256); axis tight;
% shading faceted;

for i=1:iter
    vertex1 = vertex1*tW';
    vertex1 = vertex1 ./ repmat( sqrt(sum(vertex1.^2,1)), [3 1] );

end;


surf.vertices=vertex1';
surf.faces=faces';



%-----------------------
% FUNCTION CALLS

function A = make_sparse(i,j,x, n, m)

% make_sparse - synonymous with sparse(i,j,x)


if nargin<3
    A = sparse(i,i);
    return;
end

if nargin<4
    n = max(i);
end
if nargin<5
    m = max(j);
end

A = sparse( i, j, x, n, m );


function [normal,normalf] = compute_normal(vertex,face)

% compute_normal - compute the normal of a triangulation
%
%   [normal,normalf] = compute_normal(vertex,face);
%
%   normal(i,:) is the normal at vertex i.
%   normalf(j,:) is the normal at face j.
%
%   Copyright (c) 2004 Gabriel Peyré

% [vertex,face] = check_face_vertex(vertex,face);

nface = size(face,2);
nvert = size(vertex,2);
normal = zeros(3, nvert);

% unit normals to the faces
normalf = crossp( vertex(:,face(2,:))-vertex(:,face(1,:)), ...
                  vertex(:,face(3,:))-vertex(:,face(1,:)) );
d = sqrt( sum(normalf.^2,1) ); d(d<eps)=1;
normalf = normalf ./ repmat( d, 3,1 );

% unit normal to the vertex
normal = zeros(3,nvert);
for i=1:nface
    f = face(:,i);
    for j=1:3
        normal(:,f(j)) = normal(:,f(j)) + normalf(:,i);
    end
end
% normalize
d = sqrt( sum(normal.^2,1) ); d(d<eps)=1;
normal = normal ./ repmat( d, 3,1 );

% enforce that the normal are outward
v = vertex - repmat(mean(vertex,1), 3,1);
s = sum( v.*normal, 2 );
if sum(s>0)<sum(s<0)
    % flip
    normal = -normal;
    normalf = -normalf;
end

function z = crossp(x,y)
% x and y are (m,3) dimensional
z = x;
z(1,:) = x(2,:).*y(3,:) - x(3,:).*y(2,:);
z(2,:) = x(3,:).*y(1,:) - x(1,:).*y(3,:);
z(3,:) = x(1,:).*y(2,:) - x(2,:).*y(1,:);

function A = unique_rows(A)

% unique_row - unique rows of a matrix
%
%   A = unique(A, 'rows');
%
%   For scilab/matlab compatibility
%
%   Copyright (c) 2009 Gabriel Peyre

A = unique(A, 'rows');
