gtAnalyseGrainBoundarySurface.m

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

%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
%
% return pole_data - a structure containing the pole data separated into each boundary

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

for rf = 1:length(reduce_factor)

%down sample mesh using reducepatch
fv=reducepatch(fv_in, reduce_factor(rf));


%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=[];

for k = 1:length(bounds_list)
    
    poles=[];
    poles_areas=[];

    dum=find(all(faces_bounds==bounds_list(k), 2));
    if ~isempty(dum)
        
        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
            if ~isreal(areas(ii))
                areas(ii) = 0;
            end
            
            % Convert from lab to crystallographic plane/pole
            poles(end+1,:) = gtVectorLab2Cryst(normals(ii, :), g);
            poles_areas(end+1) = areas(ii);

        end


        % Apply the symmetry operators
        poles2 = [];
        for ii = 1:length(symm)
            poles2 = [poles2; poles*symm(ii).g3];
        end
        % Add this to correspond to line in gtMakePoleFigure
        poles2=[poles2; -poles2];
        poles=poles2;

        poles_areas=repmat(poles_areas', length(poles), 1);



        %%%%%%%%%% FIGURE PLOTTING %%%%%%%%%%%
        %plot a pole figure of this boundary?
        if 0
            gtMakePoleFigure(poles, 'weights', poles_areas, 'inc', 0.01, 'searchR', 0.03)
            title(sprintf('boundary %d - mesh reduced %f', bounds_list(k), reduce_factor(rf)))
        end
        %make some kind of fancy pants coloured 3d figure
        if 0
            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
        %%%%%%%%%% FIGURE PLOTTING %%%%%%%%%%%


        %keep all pole data for the grain
        poles_all=[poles_all; poles];
        poles_areas_all=[poles_areas_all; poles_areas];

    end
    
    pole_data(bounds_list(k)).poles=poles;
    pole_data(bounds_list(k)).weights=poles_areas;
    
end

end %end of function