From c25921e437cde34190390533ce289f0f1cbaef74 Mon Sep 17 00:00:00 2001
From: Laura Nervo <laura.nervo@esrf.fr>
Date: Thu, 24 Jan 2013 16:58:00 +0100
Subject: [PATCH] General commenting and formatting, modified inputs, unified
mis_angle calculation for cubic and hexagonal
Signed-off-by: Laura Nervo <laura.nervo@esrf.fr>
---
zUtil_Twins/gtTwinTest.m | 283 ++++++++++++++++-----------------------
1 file changed, 115 insertions(+), 168 deletions(-)
diff --git a/zUtil_Twins/gtTwinTest.m b/zUtil_Twins/gtTwinTest.m
index d093093c..fa174d71 100644
--- a/zUtil_Twins/gtTwinTest.m
+++ b/zUtil_Twins/gtTwinTest.m
@@ -1,9 +1,13 @@
-function out = gtTwinTest(phaseid, parameters, grain1, grain2, grain)
-% out = gtTwinTest(phaseid, parameters, grain1, grain2, grain)
+function out = gtTwinTest(phaseid, parameters, grain1, grain2, grain, vol_out, id)
+% out = gtTwinTest(phaseid, parameters, grain1, grain2, grain, vol_out, id)
+% grain1 = grainID for grain 1
+% grain2 = grainID for grain 2
+% grain = grain structure <cell>
+% vol_out = output from gt_find_grain_boundaries
-pxsize = (parameters.labgeo.pixelsizeu + parameters.labgeo.pixelsizev)/2;
crystal_system = parameters.cryst(phaseid).crystal_system;
spacegroup = parameters.cryst(phaseid).spacegroup;
+latticepar = parameters.cryst(phaseid).latticepar;
% find the max grainID
maxgrain = max(0, grain1);
@@ -41,185 +45,128 @@ end
% this only works for cubic systems until someone works out how to pdo te
% permutations of the g1 for the hexagonal system
-if spacegroup==225 || spacegroup==229
-
- % calculate the misorientation angle and axis between these grains
- if ~isempty(g1) && ~isempty(g2)
- %if we have both grains
- %need to search symmetry equivelents for the minimum misorientation
- rot_offset=[];
- % warning('g*symm.g maybe should be symm.g*g')
- for j=1:24
- g1equiv = g1*symm(j).g;
- netg = g1equiv\g2;
- rot_offset(j) = acos((trace(netg)-1)/2);
- end
-
- dummy = find(rot_offset ==min(rot_offset));
- g1equiv = g1*symm(dummy).g;
+
+% same calculation but using hexagonal symm operations
+% calculate the misorientation angle and axis between these grains
+if ~isempty(g1) && ~isempty(g2)
+ %if we have both grains
+ %need to search symmetry equivelents for the minimum misorientation
+ rot_offset=[];
+ % warning('g*symm.g maybe should be symm.g*g')
+ for jj=1:length(symm)
+ g1equiv = g1*symm(jj).g3;
netg = g1equiv\g2;
-
- misorientation_angle = acosd((trace(netg)-1)/2);
-
- [eigVec, eigVal]=eig(netg);
- for p=1:3
- if isreal(eigVec(:,p)) %find the real eigenvector - should be just one
- misorientation_axis=eigVec(:,p);
- end
- end
+ rot_offset(jj) = acos((trace(netg)-1)/2);
end
-
-elseif spacegroup==167 || spacegroup==663 || spacegroup==194
- % same calculation but using hexagonal symm operations
- % calculate the misorientation angle and axis between these grains
- if ~isempty(g1) && ~isempty(g2)
- %if we have both grains
- %need to search symmetry equivelents for the minimum misorientation
- rot_offset=[];
- % warning('g*symm.g maybe should be symm.g*g')
- for j=1:length(symm)
- g1equiv = g1*symm(j).g3;
- netg = g1equiv\g2;
- rot_offset(j) = acos((trace(netg)-1)/2);
- end
-
- dummy = find(rot_offset ==min(rot_offset));
- g1equiv = g1*symm(dummy).g3;
- netg = g1equiv\g2;
-
- misorientation_angle = acosd((trace(netg)-1)/2);
-
- [eigVec, eigVal]=eig(netg);
- for p=1:3
- if isreal(eigVec(:,p)) %find the real eigenvector - should be just one
- misorientation_axis=eigVec(:,p);
- eigVal=eigVal(p,p);
- end
+
+ dummy = find(rot_offset ==min(rot_offset));
+ g1equiv = g1*symm(dummy).g3;
+ netg = g1equiv\g2;
+
+ misorientation_angle = acosd((trace(netg)-1)/2);
+
+ [eigVec, eigVal]=eig(netg);
+ for p=1:3
+ if isreal(eigVec(:,p)) %find the real eigenvector - should be just one
+ misorientation_axis=eigVec(:,p);
+ eigVal=eigVal(p,p);
end
- end
-end
+ end
+end
-out.grain1 = grain1;
-out.grain2 = grain2;
+
+out.grain1 = grain1;
+out.grain2 = grain2;
out.grain1_R_vector = r_vectors(grain1,2:4);
out.grain2_R_vector = r_vectors(grain2,2:4);
-out.mis_angle = misorientation_angle;
-out.mis_axis = misorientation_axis';
-
-hkl=misorientation_axis';
+out.mis_angle = misorientation_angle;
+out.mis_axis = misorientation_axis';
-hkil=[];
-hkil(1) = hkl(1) - (1/sqrt(3))*hkl(2);
-hkil(2) = (2/sqrt(3))*hkl(2);
-hkil(3) = -(hkil(1) + hkil(2));
-hkil(4) = hkl(3);
-% account for lattice vector lengths
-hkil(1:3) = hkil(1:3) / (2/(sqrt(3)*parameters.cryst(phaseid).latticepar(1)));
-hkil(4) = hkil(4) / (1/parameters.cryst(phaseid).latticepar(3));
-out.mis_axis_hex = hkil;
-out.symm = symm;
+out.mis_axis_hex = gtCart2Hex(misorientation_axis', latticepar);
+out.symm = symm;
-out.angle = rad2deg(atan2(norm(cross(R1,R2)), dot(R1,R2)));
-
-out.center1 = grain{grain1}.center/pxsize;
-out.center2 = grain{grain2}.center/pxsize;
-out.dist_centers = norm(out.center1-out.center2);
-out.radius1 = min(grain{grain1}.stat.bbxsmean,grain{grain1}.stat.bbysmean);
-out.radius2 = min(grain{grain2}.stat.bbxsmean,grain{grain2}.stat.bbysmean);
-out.distf = parameters.index.strategy.m.end.distf;
-out.is_inside1 = (out.dist_centers < out.distf*out.radius1);
-out.is_inside2 = (out.dist_centers < out.distf*out.radius2);
+out.center1 = grain{grain1}.center;
+out.center2 = grain{grain2}.center;
+out.dist_centers = norm(out.center1-out.center2);
+out.radius1 = min(grain{grain1}.stat.bbxsmean,grain{grain1}.stat.bbysmean);
+out.radius2 = min(grain{grain2}.stat.bbxsmean,grain{grain2}.stat.bbysmean);
+out.distf = parameters.index.strategy.m.end.distf;
+out.is_inside1 = (out.dist_centers < out.distf*out.radius1);
+out.is_inside2 = (out.dist_centers < out.distf*out.radius2);
gtDBConnect();
[int, x, y] = mym(['select avintint, avbbXsize, avbbYsize from ' parameters.acq.pair_tablename ' where grainID = ',num2str(grain1)]);
-out.avgint1 = mean(int);
+out.avgint1 = mean(int);
[int, x, y] = mym(['select avintint, avbbXsize, avbbYsize from ' parameters.acq.pair_tablename ' where grainID = ',num2str(grain2)]);
-out.avgint2 = mean(int);
-
-out.ratio1 = out.dist_centers/(out.radius1*out.distf);
-out.ratio2 = out.dist_centers/(out.radius2*out.distf);
-
-out.table = [out.mis_angle out.angle out.mis_axis_hex out.grain1 out.grain2];
-out.table_header = 'mis_angle angle mis_axis_hex grain1 grain2';
-%
-% %even if we have only one grain, can deal with the boundary plane
-% test = vol_out==i;
-% [x,y,z]=ind2sub(size(test), find(test));
-% [origin, a]=lsplane([x y z]);
-%
-%
-% gb_normal=a;
-% % gb_normal(3)=-gb_normal(3);
-%
-% gb_norm = gb_normal';
-% gb_centre = origin';
-%
-% if ~isempty(g1equiv)
-% plane1=g1equiv\gb_normal;
-% elseif ~isempty(g1)
-% plane1=g1\gb_normal;
-% else
-% plane1=[];
-% end
-%
-% if ~isempty(g2)
-% plane2=g2\gb_normal;
-% else
-% plane2=[];
-% end
-%
-% if spacegroup==225 || spacegroup==229
-% % if cubic then this is the index plane
-% gb_norm_grain1=plane1';
-% gb_norm_grain2=plane2';
-%
-% elseif spacegroup==167 || spacegroup==663 || spacegroup==194
-% % if hexagonal, keep cartesian vector, but also
-% % convert the cartesian plane normal to hexagonal four index notation
-% % ditto for misorientation axis
-%
-% gb_norm_grain1=plane1';
-% gb_norm_grain2=plane2';
-%
-% if ~isempty(plane1)
-% hkl=plane1;
-% hkil=[];
-% hkil(1) = hkl(1) - (1/sqrt(3))*hkl(2);
-% hkil(2) = (2/sqrt(3))*hkl(2);
-% hkil(3) = -(hkil(1) + hkil(2));
-% hkil(4) = hkl(3);
-% % account for lattice vector lengths
-% hkil(1:3) = hkil(1:3) / (2/(sqrt(3)*parameters.cryst(phaseid).latticepar(1)));
-% hkil(4) = hkil(4) / (1/parameters.cryst(phaseid).latticepar(3));
-%
-% gb_norm_hex_grain1=hkil;
-% else
-% gb_norm_hex_grain1=[];
-% end
-%
-% if ~isempty(plane2)
-% hkl=plane2;
-%
-% hkil=[];
-% hkil(1) = hkl(1) - (1/sqrt(3))*hkl(2);
-% hkil(2) = (2/sqrt(3))*hkl(2);
-% hkil(3) = -(hkil(1) + hkil(2));
-% hkil(4) = hkl(3);
-% % account for lattice vector lengths
-% hkil(1:3) = hkil(1:3) / (2/(sqrt(3)*parameters.cryst(phaseid).latticepar(1)));
-% hkil(4) = hkil(4) / (1/parameters.cryst(phaseid).latticepar(3));
-% gb_norm_hex_grain2=hkil;
-% else
-% gb_norm_hex_grain2=[];
-% end
-%
-%
-%
-% else
-% disp('crystallography not supported!')
-% end
+out.avgint2 = mean(int);
+
+out.ratio1 = out.dist_centers/(out.radius1*out.distf);
+out.ratio2 = out.dist_centers/(out.radius2*out.distf);
+
+out.table = [out.mis_angle out.mis_axis_hex out.mis_axis out.grain1 out.grain2];
+out.table_header = 'mis_angle mis_axis_hex mis_axis grain1 grain2';
+
+if ~isempty(vol_out)
+ %even if we have only one grain, can deal with the boundary plane
+ test = vol_out==id;
+ [x,y,z]=ind2sub(size(test), find(test));
+ [origin, a]=lsplane([x y z]);
+
+
+ gb_normal=a;
+ % gb_normal(3)=-gb_normal(3);
+
+ out.gb_norm = gb_normal';
+ out.gb_centre = origin';
+
+ if ~isempty(g1equiv)
+ plane1=g1equiv\gb_normal;
+ elseif ~isempty(g1)
+ plane1=g1\gb_normal;
+ else
+ plane1=[];
+ end
+
+ if ~isempty(g2)
+ plane2=g2\gb_normal;
+ else
+ plane2=[];
+ end
+
+ if spacegroup==225 || spacegroup==229
+ % if cubic then this is the index plane
+ out.gb_norm_grain1=plane1';
+ out.gb_norm_grain2=plane2';
+
+ elseif spacegroup==167 || spacegroup==663 || spacegroup==194
+ % if hexagonal, keep cartesian vector, but also
+ % convert the cartesian plane normal to hexagonal four index notation
+ % ditto for misorientation axis
+
+ out.gb_norm_grain1=plane1';
+ out.gb_norm_grain2=plane2';
+
+ if ~isempty(plane1)
+ hkl=plane1';
+ hkil = gtCart2Hex(hkl, parameters.cryst(phaseid).latticepar);
+ out.gb_norm_hex_grain1=hkil;
+ else
+ out.gb_norm_hex_grain1=[];
+ end
+ if ~isempty(plane2)
+ hkl=plane2';
+ hkil = gtCart2Hex(hkl, parameters.cryst(phaseid).latticepar);
+ out.gb_norm_hex_grain2=hkil;
+ else
+ out.gb_norm_hex_grain2=[];
+ end
+
+ else
+ disp('crystallography not supported!')
+ end
+end
end
--
GitLab