%This function generate a matrix describing a set of direction with pratically uniform distribution over the sphere
%This directions are obtained dividing iteratively faces of regular dodecaedron in smaller triangles 
%
%The mode parameter select the number of point returned within the following choices:
%Mode = 0 ->   12 directions: center of regular dodecaedron faces
%Mode = 1 ->   20 directions: vertex point of regular dodecaedron
%Mode = 2 ->   32 directions (Very uniform): vertex plus center of faces of regular dodecaedron
%Mode = 3 ->   60 directions: center of triangular faces obtained dividing each pentagonal dodecaedron faces to 5 triangles
%Mode = 4 ->  122 directions (Very uniform): vertex of triangular faces obtained dividing each triangle at previous step to 4 smaller triagles
%Mode = 5 ->  240 directions: center of triangular faces determined at previous step
%Mode = 6 ->  362 directions (Very uniform): vertex plus center of faces determined at previous step
%Mode = 7 ->  482 directions (Very uniform): ...
%Mode = 8 ->  960 directions: ...
%Mode = 9 -> 1442 directions (Very uniform): ...
%
%Note: This function can generate directions equivalent to Eigenmike capsules directions
%      simply calling it with "mode=2" and tilting down by 32° resulting directions
%      Anyway the order is not the same as returned by the specific script "GenerateEM32Dir.m"  
%
%> Dir = GenerateUniformTestDir(2);
%> [Dir(:,1),Dir(:,2)] = sph_rotation(Dir(:,1),Dir(:,2),0,+32/180*pi);


function TestDir = GenerateUniformTestDir(Mode)

%------------------------------------------------------------------------
%Start from reference for dodecaedron data comes from:
%http://www.csee.umbc.edu/~squire/reference/polyhedra.shtml#dodecahedron

%Dodecaedron vertex coords
%Step0_Vertex is a (Vertex x [x,y,z]) matrix -> 20 x 3
Step0_Vertex = [
     0.607,  0.000,  0.795 
     0.188,  0.577,  0.795 
    -0.491,  0.357,  0.795 
    -0.491, -0.357,  0.795 
     0.188, -0.577,  0.795 
     0.982,  0.000,  0.188 
     0.304,  0.934,  0.188 
    -0.795,  0.577,  0.188 
    -0.795, -0.577,  0.188 
     0.304, -0.934,  0.188 
     0.795,  0.577, -0.188 
    -0.304,  0.934, -0.188 
    -0.982,  0.000, -0.188 
    -0.304, -0.934, -0.188 
     0.795, -0.577, -0.188 
     0.491,  0.357, -0.795 
    -0.188,  0.577, -0.795 
    -0.607,  0.000, -0.795 
    -0.188, -0.577, -0.795 
     0.491, -0.357, -0.795 
   ];

%Dodecaedron vertex to face association
%Step0_Faces is a (Faces x Vertex) matrix  -> 12 x 5
Step0_Faces = 1+[
    0,  1,  2,  3,  4
    0,  1,  6, 10,  5
    1,  2,  7, 11,  6
    2,  3,  8, 12,  7
    3,  4,  9, 13,  8
    4,  0,  5, 14,  9
   15, 16, 11,  6, 10
   16, 17, 12,  7, 11
   17, 18, 13,  8, 12
   18, 19, 14,  9, 13
   19, 15, 10,  5, 14
   15, 16, 17, 18, 19
];



%------------------------------------------------------------------------
%Find the central point of each face dividing it in 5 triangles
%So at this step we obtain 20 + 12 = 32 vertex and 
%5 x 12 = 60 faces each with 3 associated vertex

Step1_Vertex = zeros(32,3);
Step1_Vertex(1:20,:) = Step0_Vertex; 
Step1_Faces = zeros(60,3);

for f=1:1:12
    %For each face
    
    %Find the coords of the face center    
    FaceCenter = mean( [   Step0_Vertex( Step0_Faces(f,1) ,:)
                           Step0_Vertex( Step0_Faces(f,2) ,:)
                           Step0_Vertex( Step0_Faces(f,3) ,:)
                           Step0_Vertex( Step0_Faces(f,4) ,:)
                           Step0_Vertex( Step0_Faces(f,5) ,:)
                        ]);
                    
    %Push the new vertex on the unitary sphere
    [theta,phi,r] = cart2sph(FaceCenter(1),FaceCenter(2),FaceCenter(3));
    [FaceCenter(1),FaceCenter(2),FaceCenter(3)] = sph2cart(theta,phi,1);
     
    %Add the new vertex
    Step1_Vertex(20+f,:) = FaceCenter;
    
    %Add the 5 triangular faces generated by the new vertex 
    Step1_Faces((f-1)*5+1,:) = [ 20+f, Step0_Faces(f,1), Step0_Faces(f,2) ];
    Step1_Faces((f-1)*5+2,:) = [ 20+f, Step0_Faces(f,2), Step0_Faces(f,3) ];
    Step1_Faces((f-1)*5+3,:) = [ 20+f, Step0_Faces(f,3), Step0_Faces(f,4) ];
    Step1_Faces((f-1)*5+4,:) = [ 20+f, Step0_Faces(f,4), Step0_Faces(f,5) ];
    Step1_Faces((f-1)*5+5,:) = [ 20+f, Step0_Faces(f,5), Step0_Faces(f,1) ];
end



%------------------------------------------------------------------------
%Step2: Divide each triangular face to 4 new triangular faces
%described by connection of the middle point of each side
%So at this step we obtain 32 + (60*3/2) = 122 vertex because each 
%of the side middle point is shared between 2 faces and 60 x 4 = 240 faces

Step2_Vertex = zeros(122,3);
Step2_Vertex(1:32,:) = Step1_Vertex;
VertexNum = 32;
Step2_Faces = zeros(240,3);
VertexPtr = zeros(3,1);

for f=1:1:60
    %For each face
    
    for s = 1:1:3
        %For each of the 3 sides
        
        %Find the coords of side middle point    
        MiddlePoint = mean([    Step1_Vertex( Step1_Faces(f,s) ,:)
                                Step1_Vertex( Step1_Faces(f,mod(s,3)+1) ,:)
                           ]);
        [theta,phi,r] = cart2sph(MiddlePoint(1),MiddlePoint(2),MiddlePoint(3));
        [MiddlePoint(1),MiddlePoint(2),MiddlePoint(3)] = sph2cart(theta,phi,1);

        %WARNING: because each side middle point is shared between 2 faces 
        %we should check if it was just added to the vertex matrix

        %The following line compute the distance between the new vertex and each point of the Step2_Vertex matrix
        %then apply a small treshold to verify presence of coincident point and find the index of first occurrence 
        res = find( sqrt( (Step2_Vertex(:,1)-ones(122,1)*MiddlePoint(1)).^2 + (Step2_Vertex(:,2)-ones(122,1)*MiddlePoint(2)).^2 + (Step2_Vertex(:,3)-ones(122,1)*MiddlePoint(3)).^2 ) < 0.001 , 1);        
        if isempty(res) 
            %The new vertex is not present
            
            %Append to vertex matrix and set accordingly the pointer
            VertexNum = VertexNum + 1;
            Step2_Vertex(VertexNum,:) = MiddlePoint;
            VertexPtr(s) = VertexNum; 
        else
            %Point to vertex just present in matrix
            VertexPtr(s) = res;
        end
               
    end
    
    %Add the 4 triangular faces generated by the new Side vertex 
    Step2_Faces((f-1)*4+1,:) = [ Step1_Faces(f,1), VertexPtr(1), VertexPtr(3) ];
    Step2_Faces((f-1)*4+2,:) = [ Step1_Faces(f,2), VertexPtr(2), VertexPtr(1) ];
    Step2_Faces((f-1)*4+3,:) = [ Step1_Faces(f,3), VertexPtr(3), VertexPtr(2) ];
    Step2_Faces((f-1)*4+4,:) = [ VertexPtr(1),     VertexPtr(2), VertexPtr(3) ];
       
end


  
  

%------------------------------------------------------------------------
%Step3: Find the central point of each triangular face
%So at this step we obtain 122 + 240 = 362 vertex  

Step3_Vertex = zeros(362,3);
Step3_Vertex(1:122,:) = Step2_Vertex;

for f=1:1:240
    %For each face
    
    %Find the coords of the face center    
    FaceCenter = mean( [   Step2_Vertex( Step2_Faces(f,1) ,:)
                           Step2_Vertex( Step2_Faces(f,2) ,:)
                           Step2_Vertex( Step2_Faces(f,3) ,:)
                       ]);
    [theta,phi,r] = cart2sph(FaceCenter(1),FaceCenter(2),FaceCenter(3));
    [FaceCenter(1),FaceCenter(2),FaceCenter(3)] = sph2cart(theta,phi,1);
     
    %Add the new vertex to the matrix
    Step3_Vertex(122+f,:) = FaceCenter;
      
end



%------------------------------------------------------------------------
%Step4: Divide again each triangular face to 4 new triangular faces
%described by connection of the middle point of each side
%So at this step we obtain 122 + (240*3/2) = 482 vertex because each 
%of the side middle point is shared between 2 faces and 240 x 4 = 960 faces

Step4_Vertex = zeros(482,3);
Step4_Vertex(1:122,:) = Step2_Vertex;
VertexNum = 122;
Step4_Faces = zeros(960,3);
VertexPtr = zeros(3,1);

for f=1:1:240
    %For each face
    
    for s = 1:1:3
        %For each of the 3 sides
        
        %Find the coords of side middle point    
        MiddlePoint = mean([    Step2_Vertex( Step2_Faces(f,s) ,:)
                                Step2_Vertex( Step2_Faces(f,mod(s,3)+1) ,:)
                           ]);
        [theta,phi,r] = cart2sph(MiddlePoint(1),MiddlePoint(2),MiddlePoint(3));
        [MiddlePoint(1),MiddlePoint(2),MiddlePoint(3)] = sph2cart(theta,phi,1);

        %WARNING: because each side middle point is shared between 2 faces 
        %we should check if it was just added to the vertex matrix

        %The following line compute the distance between the new vertex and each point of the Step4_Vertex matrix
        %then apply a small treshold to verify presence of coincident point and find the index of first occurrence 
        res = find( sqrt( (Step4_Vertex(:,1)-ones(482,1)*MiddlePoint(1)).^2 + (Step4_Vertex(:,2)-ones(482,1)*MiddlePoint(2)).^2 + (Step4_Vertex(:,3)-ones(482,1)*MiddlePoint(3)).^2 ) < 0.001 , 1);        
        if isempty(res) 
            %The new vertex is not present
            
            %Append to vertex matrix and set accordingly the pointer
            VertexNum = VertexNum + 1;
            Step4_Vertex(VertexNum,:) = MiddlePoint;
            VertexPtr(s) = VertexNum; 
        else
            %Point to vertex just present in matrix
            VertexPtr(s) = res;
        end
               
    end
    
    %Add the 4 triangular faces generated by the new Side vertex 
    Step4_Faces((f-1)*4+1,:) = [ Step2_Faces(f,1), VertexPtr(1), VertexPtr(3) ];
    Step4_Faces((f-1)*4+2,:) = [ Step2_Faces(f,2), VertexPtr(2), VertexPtr(1) ];
    Step4_Faces((f-1)*4+3,:) = [ Step2_Faces(f,3), VertexPtr(3), VertexPtr(2) ];
    Step4_Faces((f-1)*4+4,:) = [ VertexPtr(1),     VertexPtr(2), VertexPtr(3) ];
       
end


  




switch Mode
    
    case 0
        %12 central point of dodecaedron faces
        TestCoords = Step1_Vertex(21:32,:);
    case 1
        %20 vertex point of dodecaedron
        TestCoords = Step0_Vertex;
    case 2
        %32 point coming from Step 0 + Step 1
        TestCoords = Step1_Vertex;
    case 3
        %60 center of triangular faces obtained dividing each pentagonal faces to 5 triangles
        TestCoords = zeros(60,3);
        for f=1:1:60
            TestCoords(f,:) = mean( [   Step1_Vertex( Step1_Faces(f,1) ,:)
                                        Step1_Vertex( Step1_Faces(f,2) ,:)
                                        Step1_Vertex( Step1_Faces(f,3) ,:)
                                    ]);
        end
    case 4
        %122 vertex of triangular faces obtained dividing each triangle at previous step to 4 triagle 
        TestCoords = Step2_Vertex;
    case 5
        %240 center of faces determined at previous step
        TestCoords = Step3_Vertex(122:362,:);
    case 6
        %362 point coming from Step 4 and Step 5
        TestCoords = Step3_Vertex;
    case 7 
        %482 point coming to vertex
        TestCoords = Step4_Vertex;
    case 8
        %960 center of triangular faces
        TestCoords = zeros(960,3);
        for f=1:1:960
            TestCoords(f,:) = mean( [   Step4_Vertex( Step4_Faces(f,1) ,:)
                                        Step4_Vertex( Step4_Faces(f,2) ,:)
                                        Step4_Vertex( Step4_Faces(f,3) ,:)
                                    ]);
        end
    case 9
        %1442 (960 + 482)
        TestCoords = zeros(1442,3);
        for f=1:1:960
            TestCoords(f,:) = mean( [   Step4_Vertex( Step4_Faces(f,1) ,:)
                                        Step4_Vertex( Step4_Faces(f,2) ,:)
                                        Step4_Vertex( Step4_Faces(f,3) ,:)
                                    ]);
        end
        TestCoords(961:1442,:) = Step4_Vertex;
    otherwise
end

%------------------------------------------------------------------------
%Transform from cartesian to spherical coords matrix TestDir (Dir x 2) (Dir x [Theta Phi])

TestDir = zeros(size(TestCoords,1),3);
[TestDir(:,1),TestDir(:,2),TestDir(:,3)] = cart2sph(TestCoords(:,1),TestCoords(:,2),TestCoords(:,3));
TestDir = TestDir(:,1:2);


end



