% Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
% Licensed under the MIT License. See LICENSE.txt in the project root for full license information.

% Unifies faces, vertices, edges and point world coordinates into consistent structure.

% Version History
%   2017/03/28  Marc Bosch       Initial version
function [mesh]=mesh2struct(triangle_node,node_xyz,pointspacing,gtflag)
    
    uninods=unique(triangle_node);
    a=node_xyz(1,:);    
    b=node_xyz(2,:);
    c=node_xyz(3,:);
    node_xyz=[];
    node_xyz(1,:)=a;
    node_xyz(2,:)=b;
    node_xyz(3,:)=c;
    if gtflag==0
        %perform subdivision until average point spacing reached
        %take random triangle faces and estimate its pairwise distance in X and
        %Y if larger than 2*pointspacing*pointspacing+eps, then perform
        %subdivision. Repeat until reached target pointspacing
        original_node_xyz = node_xyz;
        Nfaces= max(size(triangle_node,2),size(triangle_node,1));

        Ntest=500;
        Ndivisions=50;
        dave=[];
        randpos=randperm(Nfaces);
        countvalids=0;
        for i=1:Ntest
            if max(triangle_node(:,randpos(i)))>size(original_node_xyz,2)%nodes larger than dimension are ignored
                continue;
            else                   
                countvalids=countvalids+1;
                xx=original_node_xyz(1,triangle_node(:,randpos(i)));
                yy=original_node_xyz(2,triangle_node(:,randpos(i)));
                d1=sqrt((xx(1)-xx(2))^2+(yy(1)-yy(2))^2);
                d2=sqrt((xx(1)-xx(3))^2+(yy(1)-yy(3))^2);
                d3=sqrt((xx(3)-xx(2))^2+(yy(3)-yy(2))^2);
                dave=[dave;(d1+d2+d3)/3];
            end
        end        
%         median(dave)
        if median(dave)>sqrt(pointspacing^2*2)+0.1% || ...
                %abs(Numnodesgt-size(node_xyz,2))/size(node_xyz,2)<0.5
    %         write_mesh2ply(node_xyz',triangle_node','C:\Users\boschm1\src\Algorithms\Core3D\metrics\Matlab\init1.ply');
            for j=1:Ndivisions                
                [v1 f1] = subdivide_tri( node_xyz', triangle_node' );
                %check new pointspacing
                dave=[];
                Nfaces= max(size(triangle_node,2),size(triangle_node,1));
                randpos=randperm(Nfaces);
                countvalids=0;
                for i=1:Ntest
                    if max(f1(randpos(i),:))>size(v1,1)%nodes larger than dimension are ignored
                        continue;
                    else
                        countvalids=countvalids+1;
                        xx=v1(f1(randpos(i),:),1);
                        yy=v1(f1(randpos(i),:),2);
                        d1=sqrt((xx(1)-xx(2))^2+(yy(1)-yy(2))^2);
                        d2=sqrt((xx(1)-xx(3))^2+(yy(1)-yy(3))^2);
                        d3=sqrt((xx(3)-xx(2))^2+(yy(3)-yy(2))^2);
                        dave=[dave;(d1+d2+d3)/3];
                    end
                    
                end
                 [median(dave) sqrt(pointspacing^2)+0.1]
                if median(dave)>sqrt(pointspacing^2)+0.1% || ...
                    %abs(Numnodesgt-size(v1,1))/size(v1,1)<0.5            
                    node_xyz=v1';
                    triangle_node=f1';
                else
%                     write_mesh2ply(v1,f1,'C:\Users\boschm1\src\Algorithms\Core3D\metrics\Matlab\end1.ply');
                    node_xyz=v1';
                    triangle_node=f1';
                    break;
                end
            end
        end
    end
    
    
    mesh.triangle_node=triangle_node;
    nodeuniqueIDs=unique(triangle_node);
     if length(nodeuniqueIDs)~=max(nodeuniqueIDs)
         disp('some node IDs are out of range...')
     end
  
    
    
    mesh.Nfaces = size(triangle_node,2);
    mesh.node_xyz=node_xyz;
    mesh.Nnodes = size(node_xyz,2);
    mesh.edge_node(1:mesh.Nfaces,1:2) = [triangle_node(1,:)' triangle_node(2,:)'];
    mesh.edge_node(mesh.Nfaces+1:2*mesh.Nfaces,1:2) = [triangle_node(1,:)' triangle_node(3,:)'];
    mesh.edge_node(2*mesh.Nfaces+1:3*mesh.Nfaces,1:2) = [triangle_node(2,:)' triangle_node(3,:)'];
    %remove duplicates
    %edge_node_temp=unique(mesh.edge_node,'rows');
    edge_node_temp=mesh.edge_node;
    sumval=sum(edge_node_temp,2);
    absdif=abs(edge_node_temp(:,1)-edge_node_temp(:,2));
    edgeuniquecand=[sumval absdif];
    [dummy,posuniqueedge]=unique(edgeuniquecand,'rows');
    mesh.edge_node=edge_node_temp(posuniqueedge,:)';
    mesh.Nedges = size(mesh.edge_node,2);
end


% function [mesh]=mesh2struct(triangle_node,node_xyz)
%     mesh.triangle_node=triangle_node;
%     mesh.Nfaces = size(triangle_node,2);
%     mesh.node_xyz=node_xyz;
%     mesh.Nnodes = size(node_xyz,2);
%     mesh.edge_node(1:mesh.Nfaces,1:2) = [triangle_node(1,:)' triangle_node(2,:)'];
%     mesh.edge_node(mesh.Nfaces+1:2*mesh.Nfaces,1:2) = [triangle_node(1,:)' triangle_node(3,:)'];
%     mesh.edge_node(2*mesh.Nfaces+1:3*mesh.Nfaces,1:2) = [triangle_node(2,:)' triangle_node(3,:)'];
%     %remove duplicates
%     %edge_node_temp=unique(mesh.edge_node,'rows');
%     edge_node_temp=mesh.edge_node;
%     sumval=sum(edge_node_temp,2);
%     absdif=abs(edge_node_temp(:,1)-edge_node_temp(:,2));
%     edgeuniquecand=[sumval absdif];
%     [dummy,posuniqueedge]=unique(edgeuniquecand,'rows');
%     mesh.edge_node=edge_node_temp(posuniqueedge,:)';
%     mesh.Nedges = size(mesh.edge_node,2);
% end