
% volume subsampling 

load ('volume_filled');
mand_vol = reducevolume(mand_vol_fixed,[3 3 3]);
mand_iso = isosurface(mand_vol);

clear mand_vol_fixed

%-----------------------------------------------------------------
%SMALLER MESH
%mesh subsampling
% small= reducepatch(mand_iso,0.05);
% 
% tri=small.faces;
% coord=small.vertices;
% 
% n_vertex = size(coord,1);
% n_tri = size(tri,1);
% 
% figure;trimesh(tri,coord(:,1),coord(:,2),coord(:,3));colorbar    
% 
% 
% 
% 
% 
% %center mandible
% coord=coord-kron(ones(n_vertex,1),mean(coord));
% figure;trimesh(tri,coord(:,1),coord(:,2),coord(:,3));colorbar    
% 
% 
% %picking out north and south poles
% 
% maxcood=max(coord)
% mincood=min(coord)
% 
% north=find(coord(:,2)==maxcood(2)); %multiple solutions
% north=north(1);
% south=find(coord(:,2)==mincood(2));
% south=south(1);
% 
% %
% theta=ones(n_vertex,1)*pi/2;
% heatsink=theta;
% heatsink(north)=0;
% heatsink(south)=pi;
% 
% 
% tic
% theta=laplace_smooth(theta, heatsink, tri,coord,nbr,1,10000);
% toc
% 35 seconds

%figure;trimesh(tri,coord(:,1),coord(:,2),coord(:,3),theta);colorbar    

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

%BIGGER MESH indexed by b.

tri_b=mand_iso.faces;
coord_b=mand_iso.vertices;
nb_vertex=size(coord_b,1)
nb_tri=size(tri_b,1)

coord_b=coord_b-kron(ones(nb_vertex,1),mean(coord_b));
figure;trisurf(tri_b,coord_b(:,1),coord_b(:,2),coord_b(:,3));colorbar    

%degree computation

[nbr_b, degree_b] = findnbr(tri_b);


%interploate theta values from samller mesh
theta_b = zeros(nb_vertex,1);
for i=1:nb_vertex
    diff=kron(ones(n_vertex,1),coord_b(i,:))-coord;
    diff=sum(diff.^2,2);
    theta_b(i)=theta(find(diff==min(diff)));
end;


figure;trimesh(tri_b,coord_b(:,1),coord_b(:,2),coord_b(:,3),theta_b);colorbar    

    
%picking out north and south poles

maxcood=max(coord_b)
mincood=min(coord_b)

north=find(coord_b(:,2)==maxcood(2)); %multiple solutions
north=north(1)
south=find(coord_b(:,2)==mincood(2));
south=south(1)




heatsink_b=theta_b;
heatsink_b(north)=0;
heatsink_b(south)=pi;


tic
theta_b=laplace_smooth(theta_b, heatsink_b, tri_b,coord_b,nbr_b,1,2000);
toc
% 50sec for 1000 iterations
% this part is the time sink



figure;trimesh(tri_b,coord_b(:,1),coord_b(:,2),coord_b(:,3),theta_b);colorbar    
figure;patch(coord_b(:,1),coord_b(:,2),coord_b(:,3),theta_b);colorbar    


%this gives a new representation.

ring_size=500;
ring=[];
ringind=-1*ones(ring_size,1);
for i=1:ring_size
    ind=find(abs(theta_b -pi/(ring_size-1)*(i-1))<pi/ring_size);
    if length(ind)~=0
        ringind(i)=pi/199*(i-1);
        c=mean(coord_b(ind,:),1);
        %hold on
        %plot3(c(1),c(2),c(3),'.')
        ring=[ring; c];
    end;
end;

%ringind >=0 is meaningful
%hold on
%ind=find(abs(theta_b -pi/2)<pi/50);
 %c=mean(coord_b(ind,:),1);
%plot3(c(:,1),c(:,2),c(:,3),'.');colorbar

newring=ring(1,:);
past=newring;

for i=2:size(ring,1)
    if ~isequal(past, ring(i,:))
        newring=[newring; ring(i,:)];
        past=ring(i,:);
    end
end

figure
p=patch('Faces',tri_b,'Vertices',coord_b);
set(p,'FaceColor','red','EdgeColor','none');
daspect([1 1 1])
view(3); axis tight
camlight 
lighting gouraud
alpha(0.5)
hold on
plot3(newring(:,1),newring(:,2),newring(:,3),'.')



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


%find geodesic path from the north pole to the south pole
% using the method of steepest decent.


geo_ind=north;
cur_ind=north;

while cur_ind ~=south
    cur_deg=degree_b(cur_ind);
    cur_nbr=nbr_b(cur_ind,1:cur_deg);
    
    cur_theta=theta_b(cur_nbr);
    ind=find(cur_theta==max(cur_theta));
    cur_ind=cur_nbr(ind);
    geo_ind=[geo_ind cur_ind];
end
geomax_ind=geo_ind;
geosize=size(coord_gmax,1);

coord_gmax=coord_b(geomax_ind,:);
hold on;
plot3(coord_gmax(:,1), coord_gmax(:,2),coord_gmax(:,3),'-k','LineWidth',2)


% 
% geo_ind=south;
% cur_ind=south;
% while cur_ind ~=north
%     cur_deg=degree_b(cur_ind);
%     cur_ind=nbr_b(cur_ind,1:cur_deg);
%     cur_theta=theta_b(cur_ind);
%     max_ind=find(cur_theta==min(cur_theta));
%     cur_ind=cur_ind(max_ind);
%     geo_ind=[geo_ind cur_ind];
% end
% geomin_ind=geo_ind;
% coord_gmin=coord_b(geomin_ind,:);
%hold on;
%plot3(coord_gmin(:,1), coord_gmin(:,2),coord_gmin(:,3),'-k','LineWidth',2)

ring=zeros(200,3);
ringind=-1*ones(200,1);
for i=1:200
   
    if length(ind)~=0
        ringind(i)=pi/199*(i-1);
        c=mean(coord_b(ind,:),1);
        hold on
        plot3(c(1),c(2),c(3),'.')
        ring(i,:)=c;
    end;
end;


%associate angle to center line
ang=zeros(geosize,1);
for i=1:geosize
    


p=coord_gmax(1,:);

for i=2:geosize
    
    
     ind=find(abs(theta_b -pi/(n_ring-1)*(i-1))<pi/50);
     


varphi_b=pi*ones(nb_vertex,1);

heatsink_b=varphi_b;
heatsink_b(geomax_ind)=0;

varphi_b=laplace_smooth(varphi_b, heatsink_b, tri_b,coord_b,nbr_b,1,2000);









varphi_b=zeros(nb_vertex,1);
n_ring=500;
ring=zeros(n_ring,3);
ringind=-1*ones(n_ring,1);
for i=1:200
    ind=find(abs(theta_b -pi/(n_ring-1)*(i-1))<pi/50);
    
    if length(ind)>=2
        ringind(i)=pi/199*(i-1);
        ringc=coord_b(ind,:);
        c=mean(ringc,1);
        diff=coord_gmax - kron(ones(geosize,1), c);
        diff=sum(diff.^2,2);
        rp=coord_gmax(find(diff==min(diff)),:)-c;
        crr=cross(ringc,kron(ones(size(ringc,1),1),rp));

        B=sqrt(sum(rp.^2));
        A=sqrt(sum(ringc.^2,2));

        v=asin(sqrt(sum(crr.^2,2))/B./A);
        varphi_b(ind)=v;
        
        %hold on
        %plot3(c(1),c(2),c(3),'.')
        ring(i,:)=c;
    end;
end;

figure;trimesh(tri_b,coord_b(:,1),coord_b(:,2),coord_b(:,3),varphi_b);colorbar    












geo_ind=north;
cur_ind=north;

for i=1:1000
%while cur_ind ~=maxvar
    cur_deg=degree_b(cur_ind);
    cur_nbr=nbr_b(cur_ind,1:cur_deg);
    
    cur_var=varphi_b(cur_nbr);
    ind=find(cur_var==max(cur_var));
    cur_ind=cur_nbr(ind);
    geo_ind=[geo_ind cur_ind];
end
varnoth_ind=geo_ind;
coord_varnoth=coord_b(varnoth_ind,:);


hold on;
plot3(coord_varnoth(:,1), coord_varnoth(:,2),coord_varnoth(:,3),'-k','LineWidth',2)


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




% surface normal weighted by triangle area
n=zeros(n_vertex,3);  % normal
%an=zeros(n_vertex,1); %total angle
for i=1:n_vertex
    nei= nbr(i,:);
    nei= nei(find(nei~=0));
    degree=length(nei);
    
    p=coord(nei,:)-kron(ones(degree,1),coord(i,:));
    p1=p;
    p2=[p(2:degree,:); p(1,:)];
    q=cross(p1,p2,2);
    
    area=sqrt(diag(q*q'))/2;    
    %an(i) = sum(asin(2*area./(sqrt(sum(p1.^2,2)).*sqrt(sum(p2.^2,2))))); % total angle
    area=area/sum(area);
    n(i,:)=area'*q;
end;


n=(n+0.00001)./kron(ones(1,3),sqrt(sum((n+0.00001).^2,2)));
% checking correctness
figure;trimesh(tri,coord(:,1),coord(:,2),coord(:,3));colorbar
hold on; quiver3(coord(1:10:n_vertex,1),coord(1:10:n_vertex,2),coord(1:10:n_vertex,3),...
    n(1:10:n_vertex,1),n(1:10:n_vertex,2),n(1:10:n_vertex,3))
%the above procedure is not quite correct since surface normal can have two
%directions outward and inward. So we need to figure out which way is outside.

%generating another smaller surface

coord_s=coord;
for j=1:2
    for i=1:n_vertex
        nei= nbr(i,:);
        nei= nei(find(nei~=0));
        coord_s(i,:)=mean(coord_s(nei,:));
    end;
end;

%figure;trimesh(tri,coord_s(:,1),coord_s(:,2),coord_s(:,3));colorbar    

% this gives the external direction to expand
nn=coord-coord_s;
nn=nn./kron(ones(1,3),sqrt(sum(nn.^2,2)));

%figure;trimesh(tri,coord(:,1),coord(:,2),coord(:,3));colorbar
%hold on; quiver3(coord(1:10:n_vertex,1),coord(1:10:n_vertex,2),coord(1:10:n_vertex,3),...
%    nn(1:10:n_vertex,1),nn(1:10:n_vertex,2),nn(1:10:n_vertex,3))

temp=sum(n.*nn,2);  %we find if surface normal and external direction to expand matche.
ind=find(temp<0);
n(ind,:)=-n(ind,:); 

figure;trimesh(tri,coord(:,1),coord(:,2),coord(:,3));colorbar
hold on; quiver3(coord(1:10:n_vertex,1),coord(1:10:n_vertex,2),coord(1:10:n_vertex,3),...
   n(1:10:n_vertex,1),n(1:10:n_vertex,2),n(1:10:n_vertex,3))
% this is correctly done now

figure;trimesh(tri,coord_p(:,1),coord_p(:,2),coord_p(:,3));colorbar


% area regularization using only max and min of areas.      
for i=1:n_vertex
    
    nei= nbr(i,:);
    nei= nei(find(nei~=0));
    degree=length(nei);
    
    p=coord_p(nei,:)-kron(ones(degree,1),coord_p(i,:));
    p1=p;
    p2=[p(2:degree,:); p(1,:)];
    p3=[p(degree,:); p(1:(degree-1),:)];
    q1=cross(p1,p2,2);
    q2=cross(p3,p1,2);
    
    area=(sqrt(diag(q1*q1')) +sqrt(diag(q2*q2')))/2;    
    
    maxa=max(area);
    mina=min(area);
    alpha=(maxa-mina)/maxa;
    pp=coord_p(nei(find(area==max(area))),:);
    
    coord_p(i,:)= alpha*coord_p(i,:) + (1-alpha)*pp;
    
end;

%figure;trimesh(tri,coord_p(:,1),coord_p(:,2),coord_p(:,3),dist);colorbar

for j=1:4
    for i=1:n_vertex
        nei= nbr(i,:);
        nei= nei(find(nei~=0));
        coord_p(i,:)=mean(coord_p(nei,:));
    end;
end;

figure;trimesh(tri,coord(:,1),coord(:,2),coord(:,3));colorbar





































