sugr_bR_from_E_uv

% b and R from E using set of pairs of directions

 Alg. 20 

 [b, R, code] = sugr_bR_from_E_uv(E,u,v)

 E     = 3x3 E-matrix
 u,v   = Nx3 matrices of directions

 b     = 3x1 normalized base vector
 R     = 3x3 rotation matrix
 code  = code for signes of b and R
           0 S(b)= UZU', R = VWU
           1 S(b)= UZ'U', R = VWU
           2 S(b)= UZU', R = VW'U
           3 S(b)= UZ'U', R = VW'U
            
 Wolfgang Förstner 8/2013
 wfoerstn@uni-bonn.de

 See also sugr_E_Matrix

0001 %% b and R from E using set of pairs of directions
0002 %
0003 % Alg. 20
0004 %
0005 % [b, R, code] = sugr_bR_from_E_uv(E,u,v)
0006 %
0007 % E     = 3x3 E-matrix
0008 % u,v   = Nx3 matrices of directions
0009 %
0010 % b     = 3x1 normalized base vector
0011 % R     = 3x3 rotation matrix
0012 % code  = code for signes of b and R
0013 %           0 S(b)= UZU', R = VWU
0014 %           1 S(b)= UZ'U', R = VWU
0015 %           2 S(b)= UZU', R = VW'U
0016 %           3 S(b)= UZ'U', R = VW'U
0017 %
0018 % Wolfgang Förstner 8/2013
0019 % wfoerstn@uni-bonn.de
0020 %
0021 % See also sugr_E_Matrix
0022 
0023 function [b, R, code] = sugr_bR_from_E_uv(E, u, v)
0024 
0025 % set basic matrices
0026 W = [0 1 0; - 1, 0, 0; 0, 0, 1];
0027 Z = [0 1 0; - 1, 0, 0; 0, 0, 0];
0028 I = size(u, 1);
0029 
0030 % svd of E
0031 [U, S, V] = svd(E);
0032 % enforce U and V to be proper rotations
0033 U = U * det(U);
0034 V = V * det(V);
0035 % check sign of b and R
0036 code = - 1;
0037 count = 0;
0038 for c = 0:3 % for all four cases set b and R
0039  
0040     count = count + 1;
0041     switch c
0042         case 0
0043             S = U * Z'*U'; R = V * W'*U';
0044         case 1
0045             S = U * Z * U';      R=V*W' * U';
0046         case 2
0047             S = U * Z'*U'; R = V * W * U';
0048         case 3
0049             S = U * Z * U';      R=V*W*U';
0050     end
0051     b = [S(3, 2); S(1, 3); S(2, 1)];
0052     b = b / norm(b);
0053     % check whether all 3D points are in direction of u and v
0054     sign_s = zeros(I, 1);
0055     sign_r = zeros(I, 1);
0056     for i = 1:I
0057         ui = u(i, :)';
0058         vi = v(i, :)';
0059         wi = R'*vi;
0060         m = cross(cross(b, ui), b);
0061         sign_s(i) = sign(det([b, m, cross(ui, wi)]));
0062         sign_r(i) = sign_s(i) * sign(m'*wi);
0063     end
0064     % check: the majority of points need to be in direction of u and v
0065     %     signs = [sign_s,sign_r];
0066     %     correct_sign = [mean(sign_s),mean(sign_r)];
0067     if mean(sign_s) > 0 && mean(sign_r) > 0
0068         code = c;
0069         return
0070     end
0071 end
0072 
0073 
0074         
0075