minor fixes / adding slow (non-parallelized version) of re-mapping into fundamental zone

Signed-off-by: Wolfgang Ludwig <wolfgang.ludwig@esrf.fr>

zUtil_Maths/gtMathsOriMat2Rod.m

@@ -26,22 +26,29 @@ elseif (size(g, 1)~=3 || size(g, 2)~=3)
 end
 
 % Set epsilon to avoid singular values
-epsilon = 1.e-9;
+epsilon = 1.e-10;
 
 % Compute orientation matrices traces
 g_trace_plus1 = g(1, 1, :) + g(2, 2, :) + g(3, 3, :) + 1;
 
 % Find regular cases and prepare multiplier (null for non regular cases)
-multiplier = double(abs(g_trace_plus1) >= epsilon)./(g_trace_plus1);
+singular = abs(g_trace_plus1) < epsilon;
 
-r = [ g(2, 3, :) - g(3, 2, :) ; ...
-      g(3, 1, :) - g(1, 3, :) ; ...
-      g(1, 2, :) - g(2, 1, :) ];
+if any(singular)
+    [angle, axis] = gtMathsOriMat2AngleAxis(g);  % slower but handling the singularity
+    r = repmat(tand(squeeze(angle) / 2)', 3, 1).* axis;
+else
+    multiplier = double(abs(g_trace_plus1) >= epsilon)./(g_trace_plus1);
 
-r = r .* multiplier([1 1 1], :, :);
- 
-% Reshape resulting matrix 
-% r = squeeze(r);
-r = reshape(r, 3, []); % Runs faster than squeeze
+    r = [ g(2, 3, :) - g(3, 2, :) ; ...
+          g(3, 1, :) - g(1, 3, :) ; ...
+          g(1, 2, :) - g(2, 1, :) ];
+
+    r = r .* multiplier([1 1 1], :, :);
+
+    % Reshape resulting matrix 
+    % r = squeeze(r);
+    r = reshape(r, 3, []); % Runs faster than squeeze
+end
 
 end % end of function

zUtil_Maths/gtMathsRod2RodInFundZone.m

 function r_vec = gtMathsRod2RodInFundZone(r_vec, symm, shift)
+    if (exist('shift', 'var'))
+        r_vec = fast_map(r_vec, symm, shift);
+    else
+        r_vec = slow_map(r_vec, symm);
+    end
+end
+
+function r_vec = fast_map(r_vec, symm, shift)
     num_symm_ops = 2 * numel(symm);
 
     % Expanding the symmetry operators and the matrix to be checked, to
@@ -43,3 +51,24 @@ function r_vec = gtMathsRod2RodInFundZone(r_vec, symm, shift)
         r_vec(:, iis) = r_vecs(:, ind_min);
     end
 end
+
+function r_vec = slow_map(r_vec, symm)
+    g = gtMathsRod2OriMat(r_vec);
+
+    num_symm_ops = numel(symm);
+
+    r_vecs = zeros(2 * num_symm_ops, 3);
+
+    for ii = 1:num_symm_ops
+        r_vecs((ii - 1) * 2 + 1, :) = gtMathsOriMat2Rod(symm(ii).g3 * g);
+        r_vecs((ii - 1) * 2 + 2, :) = gtMathsOriMat2Rod(symm(ii).g3' * g);
+    end
+
+    dists = 2 * atand(sqrt(sum(r_vecs .^ 2, 2)));
+
+    [~, ii] = min(dists);
+
+%     r_vecs
+
+    r_vec = r_vecs(ii, :)';
+end