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

## CROSS-REFERENCE INFORMATION

This function calls:
• essential_flat Reshape a [3x6xk] matrix to a [3x3x2k] matrix
• essential_hat3 Compute the matrix representation of the cross product
• essential_sharp Reshape a [3x3x2k] matrix to a [3x6xk] matrix
• essential_closestRepresentative
• randrot Generates uniformly random rotation matrices.
• randskew Generates random skew symmetric matrices with normal entries.
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
• multiprod Multiplying 1-D or 2-D subarrays contained in two N-D arrays.
• multiskew Returns the skew-symmetric parts of the matrices in the 3D matrix X.
• multisym Returns the symmetric parts of the matrices in the 3D matrix X
• multitrace Computes the traces of the 2D slices in a 3D matrix.
• multitransp Transposing arrays of matrices.
This function is called by:
• essential_svd Sample solution of an optimization problem on the essential manifold.

## 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 %
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
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
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 Fri 08-Sep-2017 12:43:19 by m2html © 2005