function surf=REGULARIZEarea(surf, alpha);
% Regularize area of triangles using only larger traingles with given ratio alpha
% (c) 2012 Moo K. Chung
% mkchung@wisc.edu


% This code is based on 
%
% 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. 

%surf=amygsurf;
n_vertex=size(surf.vertices,1);
n_faces=size(surf.faces,1);
[nbr, degree] = findnbr(surf.faces);

area=GETarea(surf);

%figure;plot(area,'.');
areac=mean(area);

%c=mean(surf.vertices);
%r=surf.vertices-kron(ones(n_vertex,1),c);
%radius = mean(L2norm(r));

% determine index above alpha-percent quantile.
areaq=quantile(area,alpha);
ind=find(area > areaq);  % larger area

%plot(area,'.')

for j=1:length(ind)
    p=surf.vertices(surf.faces(ind(j),:),:);

    p1=p(1,:);
    p2=p(2,:);
    p3=p(3,:);
    pc=mean(p);    % center of the triangle

    ratio=sqrt(areac/area(ind(j))); % shrinkage ratio. if ratio=1, do not move. 
    
    p1=  ratio*p1 + (1-ratio)*pc;
    p2=  ratio*p2 + (1-ratio)*pc;
    p3=  ratio*p3 + (1-ratio)*pc;
    
    surf.vertices(surf.faces(ind(j),1),:)=p1;
    surf.vertices(surf.faces(ind(j),2),:)=p2;
    surf.vertices(surf.faces(ind(j),3),:)=p3;
    
end;

%------------------------------------------
function [nbr, degree] = findnbr(tri)
% From tri=surf.faces, we determine the connectivity info and the
% corresponding degree of each vertex.

n_tri=size(tri,1);
n_vertex = (n_tri+4)/2;
degree=zeros(n_vertex,1);

for j=1:n_tri
    degree(tri(j,:))=degree(tri(j,:))+1;
end

max_degree=max(degree);

% find out the neighboring vertices

nbr=zeros(n_vertex,max_degree);

for i_tri=1:n_tri
    for j=1:3
        cur_point = tri(i_tri,j);
        for k=1:3
            if (j ~= k)
                nbr_point= tri(i_tri,k);
                if find(nbr(cur_point,:)==nbr_point)
                    ;
                else
                    n_nbr = min(find(nbr(cur_point,:) == 0));
                    nbr(cur_point,n_nbr) = nbr_point;
                end;
            end;
        end;
    end;
end;

%--------------------------------------------------------

function area=GETarea(surf);
% from surf, computes the area of each triangle element.
coord=surf.vertices;
tri=surf.faces;
n_faces=size(surf.faces,1);
area=zeros(n_faces,1);

for i=1:n_faces
        
   p = coord(tri(i,:),:);
   p1 = p(1,:)-p(3,:);   
   p2 = p(2,:)-p(3,:);
   q = cross(p1,p2);
   area(i)= sqrt(q*q')/2;
end;
