General commenting and formatting, modified inputs, unified mis_angle...

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

-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