-
Signed-off-by:
Laura Nervo <laura.nervo@esrf.fr> git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@577 4c865b51-4357-4376-afb4-474e03ccb993
Signed-off-by:
Laura Nervo <laura.nervo@esrf.fr> git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@577 4c865b51-4357-4376-afb4-474e03ccb993
gtAnalyseGrainBoundarySurface.m 3.98 KiB
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