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


sym=gtGetCubicSymOp;
R=r_vectors(grainID,2:4);
g=Rod2g(R);

%if no list supplied which boundaries are relevent?
if isempty(bounds_list)
for i=1:length(boundaries_structure)
    if boundaries_structure(i).grain1==grainID || boundaries_structure(i).grain2==grainID
        bounds_list=[bounds_list i];
    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 i=1:length(fv.vertices)
    a=round(fv.vertices(i,:));
    boundaryID=bounds_dil(a(2), a(1), a(3));
    faces_bounds(find(fv.faces==i))=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 j=1:length(dum)
            i=dum(j);
            face=fv.faces(i,:);
            %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(i,:)=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(i)=norm(v1-v2)*norm(v3-v2)*sin(angle)*0.5;

            %can get co-linear points
            if ~isreal(areas(i))
                areas(i)=0;
            end
            
            
            %convert to crystallographic plane/pole
            poles(end+1,:)=inv(g)*normals(i,:)';
            poles_areas(end+1)=areas(i);

        end


        %apply the symmetry operators
        poles2=[];
        for i=1:length(sym)
            poles2=[poles2; poles*sym(i).g];
        end
        %add this to correspond to line in gtMakePoleFigure
        poles2=[poles2; -poles2];
        poles=poles2;

        poles_areas=repmat(poles_areas', 48, 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