Home > manopt > manifolds > essential > essentialfactory.m

essentialfactory

PURPOSE ^

Manifold structure to optimize over the space of essential matrices.

SYNOPSIS ^

function M = essentialfactory(k, strSigned)

DESCRIPTION ^

 Manifold structure to optimize over the space of essential matrices.

 function M = essentialfactory(k)
 function M = essentialfactory(k, 'signed')
 function M = essentialfactory(k, 'unsigned')


 Quotient representation of the essential manifold: deals with the
 representation of the space of essential matrices M_rE. These are used in
 computer vision to represent the epipolar constraint between projected
 points in two perspective views.

 The space is represented as the quotient (SO(3)^2/SO(2)).
 See the following references for details:

   R. Tron, K. Daniilidis,
   "On the quotient representation of the essential manifold"
   IEEE Conference on Computer Vision and Pattern Recognition, 2014

 For computational purposes, each essential matrix is represented as a
 [3x6] matrix where each [3x3] block is a rotation.

 The metric used is the one induced by the submersion of M_rE in SO(3)^2.

 Tangent vectors are represented in the Lie algebra of SO(3)^2, i.e., as
 [3x6] matrices where each [3x3] block is a skew-symmetric matrix.
 Use the function tangent2ambient(X, H) to switch from the Lie algebra
 representation to the embedding space representation in R^(3x6).

 By default, k = 1, and the geometry is 'signed'.

 Optional arguments:
   "signed"    selects the signed version of the manifold
   "unsigned"  selects the unsigned version of the manifold

 See also rotationsfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = essentialfactory(k, strSigned)
0002 % Manifold structure to optimize over the space of essential matrices.
0003 %
0004 % function M = essentialfactory(k)
0005 % function M = essentialfactory(k, 'signed')
0006 % function M = essentialfactory(k, 'unsigned')
0007 %
0008 %
0009 % Quotient representation of the essential manifold: deals with the
0010 % representation of the space of essential matrices M_rE. These are used in
0011 % computer vision to represent the epipolar constraint between projected
0012 % points in two perspective views.
0013 %
0014 % The space is represented as the quotient (SO(3)^2/SO(2)).
0015 % See the following references for details:
0016 %
0017 %   R. Tron, K. Daniilidis,
0018 %   "On the quotient representation of the essential manifold"
0019 %   IEEE Conference on Computer Vision and Pattern Recognition, 2014
0020 %
0021 % For computational purposes, each essential matrix is represented as a
0022 % [3x6] matrix where each [3x3] block is a rotation.
0023 %
0024 % The metric used is the one induced by the submersion of M_rE in SO(3)^2.
0025 %
0026 % Tangent vectors are represented in the Lie algebra of SO(3)^2, i.e., as
0027 % [3x6] matrices where each [3x3] block is a skew-symmetric matrix.
0028 % Use the function tangent2ambient(X, H) to switch from the Lie algebra
0029 % representation to the embedding space representation in R^(3x6).
0030 %
0031 % By default, k = 1, and the geometry is 'signed'.
0032 %
0033 % Optional arguments:
0034 %   "signed"    selects the signed version of the manifold
0035 %   "unsigned"  selects the unsigned version of the manifold
0036 %
0037 % See also rotationsfactory
0038 
0039 % Please cite the Manopt paper as well as the research paper:
0040 %     @InProceedings{tron2014essential,
0041 %       Title        = {On the quotient representation of the essential manifold},
0042 %       Author       = {Tron, R. and Daniilidis, K.},
0043 %       Booktitle    = {IEEE Conference on Computer Vision and Pattern Recognition},
0044 %       Year         = {2014},
0045 %       Organization = {{IEEE CVPR}}
0046 %     }
0047 
0048 
0049 % This file is part of Manopt: www.manopt.org.
0050 % Original author: Roberto Tron, Aug. 8, 2014
0051 % Contributors: Bamdev Mishra, May 15, 2015.
0052 %
0053 %
0054 % RT: General implementation note: to streamline component-wise
0055 % computations, in tangentProjection and exponential,
0056 % we flatten out the arguments into [3 x 3 x 2K] arrays, compute the
0057 % components all together, and then sharp the result again into [3 x 6 x K]
0058 % arrays.
0059 
0060 
0061     % Optional parameters to switch between the signed and unsigned
0062     % versions of the manifold.
0063     if ~exist('strSigned', 'var') || isempty(strSigned)
0064         strSigned = 'signed';
0065     end
0066     switch(strSigned)
0067         case 'signed'
0068             flagSigned = true;
0069         case 'unsigned'
0070             flagSigned = false;
0071         otherwise
0072             error('Second argument can be either empty, ''signed'', or ''unsigned''.');
0073     end
0074 
0075     
0076     if ~exist('k', 'var') || isempty(k)
0077         k = 1;
0078     end
0079     
0080     if k == 1
0081         M.name = @() sprintf('Quotient representation of the essential manifold, %s', strSigned);
0082     elseif k > 1 && k == round(k)
0083         M.name = @() sprintf('Product of %d quotient representations of the essential manifold, %s', k, strSigned);
0084     else
0085         error('k must be an integer no less than 1.');
0086     end
0087     
0088     M.dim = @() k*5;
0089     
0090     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0091     
0092     M.norm = @(x, d) norm(d(:));
0093     
0094     M.typicaldist = @() pi*sqrt(2*k);
0095     
0096     M.proj = @tangentProjection;
0097     function HProjHoriz=tangentProjection(X,H)
0098         % Project H on the tangent space of SO(3)^2
0099         HProj = essential_sharp(multiskew(multiprod(multitransp(essential_flat(X)), essential_flat(H))));
0100         
0101         % Compute projection on vertical component
0102         p = vertproj(X, HProj);
0103         
0104         HProjHoriz = HProj - multiprod(p/2,[essential_hat3(permute(X(3,1:3,:),[2 3 1])) essential_hat3(permute(X(3,4:6,:),[2 3 1]))]);% BM: okay
0105     end
0106     
0107     
0108     M.tangent = @(X, H) essential_sharp(multiskew(essential_flat(H)));
0109     
0110     M.egrad2rgrad=@egrad2rgrad;
0111     function rgrad = egrad2rgrad(X, egrad)
0112         rgrad = M.proj(X, egrad);
0113     end
0114     
0115     M.ehess2rhess = @ehess2rhess;
0116     function rhess = ehess2rhess(X, egrad, ehess, S)
0117         % Reminder: S contains skew-symmeric matrices. The actual
0118         % direction that the point X is moved along is X*S.
0119         RA = p1(X);
0120         RB = p2(X);
0121         SA = p1(S);
0122         SB = p2(S);
0123         
0124         G = egrad; 
0125         GA = p1(G);
0126         GB = p2(G);
0127         
0128         H = ehess; 
0129         
0130         % RT: We now compute the connection, i.e. the part of the derivative
0131         % given by the curvature of the space (as opposed to a simple
0132         % Euclidean derivative).
0133         
0134         % The following is the vectorized version of connection=-[multisym(GA'*RA)*SA multisym(GB'*RB)*SB];
0135         connection = tangent2ambient(X,-cat(2,...
0136             multiprod(multisym(multiprod(multitransp(GA), RA)), SA),...
0137             multiprod(multisym(multiprod(multitransp(GB), RB)), SB)));
0138         rhess = M.proj(X,H + connection);
0139     end
0140     
0141     
0142     
0143     M.exp = @exponential;
0144     function Y = exponential(X, U, t)
0145         if nargin == 3
0146             U = t*U;
0147         end
0148         
0149         UFlat = essential_flat(U);
0150         exptUFlat = rot3_exp(UFlat);
0151         Y = essential_sharp(multiprod(essential_flat(X), exptUFlat));
0152     end
0153     
0154     M.retr = @exponential;
0155     
0156     M.log = @logarithm; 
0157     function U = logarithm(X, Y)
0158         
0159         QX = [X(:,1:3,:);X(:,4:6,:)];
0160         QY = [Y(:,1:3,:);Y(:,4:6,:)];
0161         QYr = essential_closestRepresentative(QX,QY,'flagSigned',flagSigned);
0162         Yr = [QYr(1:3,:,:) QYr(4:6,:,:)];
0163         U = zeros(size(X));
0164         U(:,1:3,:) = rot3_log(multiprod(multitransp(X(:,1:3,:)),Yr(:,1:3,:)));
0165         U(:,4:6,:) = rot3_log(multiprod(multitransp(X(:,4:6,:)),Yr(:,4:6,:)));
0166     end
0167     
0168     M.hash = @(X) ['z' hashmd5(X(:))];
0169     
0170     M.rand = @() randessential(k);
0171     function Q = randessential(N)
0172         % Generates random essential matrices.
0173         %
0174         % function Q = randessential(N)
0175         %
0176         % Q is a [3x6] matrix where each [3x3] block is a uniformly distributed
0177         % matrix.
0178         
0179         % This file is part of Manopt: www.manopt.org.
0180         % Original author: Roberto Tron, Aug. 8, 2014
0181         % Contributors:
0182         % Change log:
0183         
0184         if nargin < 1
0185             N = 1;
0186         end
0187         
0188         Q = [randrot(3,N) randrot(3,N)];
0189     end
0190     
0191     M.randvec = @randomvec;
0192     function U = randomvec(X)
0193         U = tangentProjection(X, essential_sharp(randskew(3, 2*k)));
0194         U = U / sqrt(M.inner([],U,U));
0195     end
0196     
0197     M.lincomb = @matrixlincomb;
0198     
0199     M.zerovec = @(x) zeros(3, 6, k);
0200     
0201     M.transp = @transport;
0202     function S2 = transport(X1, X2, S1)
0203         % Transport a vector from the tangent space at X1 to the tangent
0204         % space at X2. This transport uses the left translation of the
0205         % ambient group and preserves the norm of S1. The left translation
0206         % aligns the vertical spaces at the two elements.
0207         
0208         % Group operation in the ambient group, X12=X2'*X1
0209         X12 = essential_sharp(multiprod(multitransp(essential_flat(X2)),essential_flat(X1)));
0210         X12Flat = essential_flat(X12);
0211         
0212         % Left translation, S2=X12*S*X12'
0213         S2 = essential_sharp(multiprod(X12Flat,multiprod(essential_flat(S1),multitransp(X12Flat))));
0214     end
0215     
0216     M.pairmean = @pairmean;
0217     function Y = pairmean(X1, X2)
0218         V = M.log(X1, X2);
0219         Y = M.exp(X1, .5*V);
0220     end
0221     
0222     M.dist = @(x, y) M.norm(x, M.log(x, y)); 
0223     
0224     M.vec = @(x, u_mat) u_mat(:);
0225     M.mat = @(x, u_vec) reshape(u_vec, [3, 6, k]);
0226     M.vecmatareisometries = @() true;
0227     
0228     
0229     
0230     p1 = @(X) X(:,1:3,:);
0231     p2 = @(X) X(:,4:6,:);
0232     
0233     
0234     vertproj = @(X,H) multiprod(X(3,1:3,:),permute(vee3(H(:,1:3,:)),[1 3 2]))+multiprod(X(3,4:6,:),permute(vee3(H(:,4:6,:)),[1 3 2]));
0235     
0236     tangent2ambient = @(X, H) essential_sharp(multiprod(essential_flat(X), essential_flat(H)));
0237     
0238     
0239 end
0240 
0241 
0242 %% Some functions used by the essential factory
0243 
0244 function v = vee3(V)
0245     v = squeeze([V(3,2,:)-V(2,3,:); V(1,3,:)-V(3,1,:); V(2,1,:)-V(1,2,:)])/2;
0246 end
0247 
0248 
0249 % Compute the exponential map in SO(3) using Rodrigues' formula
0250 %  function R = rot3_exp(V)
0251 % V must be a [3x3xN] array of [3x3] skew-symmetric matrices.
0252 function R = rot3_exp(V)
0253     v = vee3(V);
0254     nv = cnorm(v);
0255     idxZero = nv < 1e-15;
0256     nvMod = nv;
0257     nvMod(idxZero) = 1;
0258     
0259     vNorm = v./([1;1;1]*nvMod);
0260     
0261     % Matrix exponential using Rodrigues' formula
0262     nv = shiftdim(nv,-1);
0263     c = cos(nv);
0264     s = sin(nv);
0265     [VNorm,vNormShift] = essential_hat3(vNorm);
0266     vNormvNormT = multiprod(vNormShift,multitransp(vNormShift));
0267     R=multiprod(eye(3),c)+multiprod(VNorm,s)+multiprod(vNormvNormT,1-c);
0268 end
0269 
0270 
0271 
0272 % Compute the logarithm map in SO(3)
0273 %  function V = rot3_log(R)
0274 % V is a [3x3xN] array of [3x3] skew-symmetric matrices
0275 function V = rot3_log(R)
0276     skewR = multiskew(R);
0277     ctheta = (multitrace(R)'-1)/2;
0278     stheta = cnorm(vee3(skewR));
0279     theta = atan2(stheta,ctheta);
0280     
0281     V=skewR;
0282     for ik=1:size(R,3)
0283         V(:,:,ik)=V(:,:,ik)/sincN(theta(ik));
0284     end
0285 end
0286 
0287 
0288 function sx = sincN(x)
0289     sx = sin(x)./x;
0290     sx(x==0) = 1;
0291 end
0292 
0293 function nv = cnorm(v)
0294     nv = sqrt(sum(v.^2));
0295 end
0296 
0297 
0298

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005