gtAnalyseGrainBoundarySurface2.m

function [poles_all,poles_areas_all]=gtAnalyseGrainBoundarySurface2(fv_in, reduce_factor, grainID, bounds_list, boundaries_structure, r_vectors, bounds_dil)

%version 2 - rather than return the poles/pole figure for the boundaries,
%colour code the boudaries according to the index plane ( see gtSSTCmap )


%nb - this was part of gtAnalyseBoundaryIndexPlaneMesh.m
%this has been split into this file (analyse mesh) and
% gtGrainBoundarySurfaceMesh (make mesh)
%
%grain_surface is fv , the facevertex data for patch
%
% bounds_list - which boundaries to analyse
% reduce_factor (should be less than 1) passed to reduce patch

symm = gtCrystGetSymmetryOperators('cubic');
R = r_vectors(grainID, 2:4);
g = gtMathsRod2RotMat(R);

%if no list supplied which boundaries are relevent?
if isempty(bounds_list)
    for ii = 1:length(boundaries_structure)
        if boundaries_structure(ii).grain1==grainID || boundaries_structure(ii).grain2==grainID
            bounds_list=[bounds_list ii];
        end
    end
end


%reduce mesh - single value only
fv=reducepatch(fv_in, reduce_factor);


%loop through vertices, assigning each to a boundary of the original grain
%may not pick up some very small boundaries, but I guess this doesn't
%matter- for grain 1 these are of size 1-39 voxels volume ie max diameter 4
%voxels
faces_bounds=zeros(size(fv.faces));
for ii = 1:length(fv.vertices)
    a = round(fv.vertices(ii, :));
    boundaryID=bounds_dil(a(2), a(1), a(3));
    faces_bounds(find(fv.faces==ii))=boundaryID;
end

%record down sampled grain surface / boundaries data?
grain_surface.fv=fv;
grain_surface.faces_bounds=faces_bounds;

%record all normals and all triangle area for weighting
normals=zeros(length(fv.faces), 3);
areas=zeros(length(fv.faces), 1);


poles_all=[];
poles_areas_all=[];

f1=figure;

%llop through the boundaries being considered, analysing the planes of the
%facets
for k = 1:length(bounds_list)

    poles=[];
    poles_areas=[];
    bad_faces=[];

    dum=find(all(faces_bounds==bounds_list(k), 2));
    if ~isempty(dum)
        
        %loop through the triangles of this boundary
        for jj = 1:length(dum)
            ii = dum(jj);
            face = fv.faces(ii, :);
            %get the three vertices -
            v1=fv.vertices(face(1), [2 1 3]);
            v2=fv.vertices(face(2),  [2 1 3]);
            v3=fv.vertices(face(3),  [2 1 3]);
            
            %normal:
            n=cross(v1-v2, v3-v2);
            %fix this to follow lsplane output - change sign of third component
            n(3)=-n(3);
            normals(ii, :) = n/norm(n);%may not need to save this info

            %area
            angle=acos(dot((v1-v2)/norm(v1-v2), (v3-v2)/norm(v3-v2)));
            areas(ii) = norm(v1-v2)*norm(v3-v2)*sin(angle)*0.5;

            %can get co-linear points - these leave NaNs in poles
            %can also get colinear points that don't causecomplex areas -
            %low angle test
            if ~isreal(areas(ii)) || abs(angle)<0.00001
                areas(ii)=0;
                bad_faces=[bad_faces; jj];
            end
            
            %convert to crystallographic plane/pole
            poles(end+1,:) = gtVectorLab2Cryst(g, normals(ii, :));
            poles_areas(end+1) = areas(ii); 
        end

        % Remove NaN faces from dum
        dum(bad_faces)=[];
        poles(bad_faces, :)=[];
        poles_areas(bad_faces)=[];
        
        
        % Apply the symmetry operators
        poles2=[];
        poles_for_pole_fig=[]
        
        for ii = 1:length(symm)
            % Deal with this in a 3d array
            poles2 = cat(3, poles2, poles*symm(ii).g);
            % Keep a 2d array for gtMakePoleFigure
            poles_for_pole_fig = [poles_for_pole_fig; poles*symm(ii).g];
        end
        % Add this to correspond to line in gtMakePoleFigure
        poles2=cat(3,poles2,-poles2);
        poles_for_pole_fig=[poles_for_pole_fig; -poles_for_pole_fig];
        poles_areas_for_pole_fig=repmat(poles_areas', 2*length(symm), 1);
     
        % Keep all pole data for the grain
        poles_all=[poles_all; poles_for_pole_fig];
        poles_areas_all=[poles_areas_all; poles_areas_for_pole_fig];
        
        % Get those poles which lie in the SST
        a=find(poles2(:,2,:)>=0 & poles2(:,2,:)<=poles2(:,1,:) & poles2(:,3,:)>=poles2(:,1,:));
        % a is effectively an index in an 2d array, size(poles2, 1) by size(poles2, 3) (48)
        [a1,a2]=ind2sub([size(poles2, 1), size(poles2, 3)], a);
        % But - this scrambles the order (i think)  should sort by a1 to keep poles3 in
        % Orignal order
        [a1, ii] = sort(a1);
        a2 = a2(ii);

        % Get the relevent poles back out of the 3d thing
        % temp use poles3
        poles3=zeros(size(poles));
        for p=1:size(poles2, 1)
            poles3(p,:)=poles2(a1(p),:,a2(p));
        end
  
        % Get the colour
        x=poles3(:,1);
        y=poles3(:,2);
        z=poles3(:,3);

        alpha=atan(x./z); % in sst alpha 0->pi/4
        phi=atan2(y,x);

        rgb=[1-(alpha/(pi/4)) (alpha/(pi/4)).*(1-(phi/(pi/4))) (alpha/(pi/4)).*(phi/(pi/4))];
        for ii = 1:size(rgb, 1)
            rgb(ii,:) = rgb(ii, :)/max(rgb(ii, :));
        end

        figure(f1);
        fv_fig.vertices=fv.vertices;
        fv_fig.faces=fv.faces(dum,:);
        p=patch(fv_fig);
        set(p, 'facevertexcdata', rgb);
        shading faceted;


        %poles=poles2;
        % Not really needed if 3d thing works
        %poles_areas=repmat(poles_areas', 48, 1);

        % Plot a pole figure of this boundary?
        if 0
            gtMakePoleFigure(poles_for_pole_fig, 'weights', poles_areas_for_pole_fig, 'inc', 0.01, 'searchR', 0.03)
            title(sprintf('boundary %d - mesh reduced %f', bounds_list(k), reduce_factor))
        end
        
%         %make some kind of fancy pants coloured 3d figure
%         %grain figure with the individual boundaries coloured
%         figure(f1);
%         fv_fig.vertices=fv.vertices;
%         fv_fig.faces=fv.faces(dum,:);
%         p=patch(fv_fig)
%         if bounds_list(k)==3076
%             set(p, 'facecolor', 'r')
%         else
%             set(p, 'facecolor', rand(3,1))
%         end
    end
end

%optional - add the facets that do not belong uniquely to one boundary,
%coloured black.
if 1
    dum=find(faces_bounds(:,1)~=faces_bounds(:,2) || faces_bounds(:,1)~=faces_bounds(:,3));
    figure(f1);
    fv_fig.vertices=fv.vertices;
    fv_fig.faces=fv.faces(dum,:);
    p=patch(fv_fig);
set(p, 'facevertexcdata', zeros(length(dum), 3))
shading faceted
end

axis equal;
axis vis3d;

end % end of function