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.
• cat CAT concatenation of two TT/MPS tensors.
• norm NORM Norm of a TT/MPS tensor.
• round ROUND Approximate TTeMPS tensor within a prescribed tolerance.
• norm NORM Norm of a TT/MPS block-mu tensor.
• round ROUND Approximate TTeMPS tensor within a prescribed tolerance.
• round ROUND Approximate TTeMPS operator within a prescribed tolerance.
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
• multiprod Matrix multiply 2-D slices of 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 a 3D array
• multitrace Computes the traces of the 2D slices in a 3D matrix.
• multitransp Transpose the matrix slices of an N-D array (no complex conjugate)
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 % Change log:
0053 %   Jan. 4, 2021 (NB):
0054 %       Compatibility with Octave 6.1.0: moved quite a few nested functions
0055 %       to end-of-file functions, and made tangent2ambient accessible as
0056 %       M.tangent2ambient (which should have been the case already).
0057
0058 % RT: General implementation note: to streamline component-wise
0059 % computations, in tangentProjection and exponential,
0060 % we flatten out the arguments into [3 x 3 x 2K] arrays, compute the
0061 % components all together, and then sharp the result again into [3 x 6 x K]
0062 % arrays.
0063
0064
0065     % Optional parameters to switch between the signed and unsigned
0066     % versions of the manifold.
0067     if ~exist('strSigned', 'var') || isempty(strSigned)
0068         strSigned = 'signed';
0069     end
0070     switch(strSigned)
0071         case 'signed'
0072             flagSigned = true;
0073         case 'unsigned'
0074             flagSigned = false;
0075         otherwise
0076             error('Second argument can be either empty, ''signed'', or ''unsigned''.');
0077     end
0078
0079
0080     if ~exist('k', 'var') || isempty(k)
0081         k = 1;
0082     end
0083
0084     if k == 1
0085         M.name = @() sprintf('Quotient representation of the essential manifold, %s', strSigned);
0086     elseif k > 1 && k == round(k)
0087         M.name = @() sprintf('Product of %d quotient representations of the essential manifold, %s', k, strSigned);
0088     else
0089         error('k must be an integer no less than 1.');
0090     end
0091
0092     M.dim = @() k*5;
0093
0094     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0095
0096     M.norm = @(x, d) norm(d(:));
0097
0098     M.typicaldist = @() pi*sqrt(2*k);
0099
0100     M.proj = @tangentProjection;
0101     function HProjHoriz=tangentProjection(X,H)
0102         % Project H on the tangent space of SO(3)^2
0103         HProj = essential_sharp(multiskew(multiprod(multitransp(essential_flat(X)), essential_flat(H))));
0104
0105         % Compute projection on vertical component
0106         p = vertproj(X, HProj);
0107
0108         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
0109     end
0110
0111
0112     M.tangent = @(X, H) essential_sharp(multiskew(essential_flat(H)));
0113
0117     end
0118
0119     M.ehess2rhess = @ehess2rhess;
0120     function rhess = ehess2rhess(X, egrad, ehess, S)
0121         % Reminder: S contains skew-symmetric matrices. The actual
0122         % direction that the point X is moved along is X*S.
0123         RA = p1(X);
0124         RB = p2(X);
0125         SA = p1(S);
0126         SB = p2(S);
0127
0129         GA = p1(G);
0130         GB = p2(G);
0131
0132         H = ehess;
0133
0134         % RT: We now compute the connection, i.e. the part of the derivative
0135         % given by the curvature of the space (as opposed to a simple
0136         % Euclidean derivative).
0137
0138         % The following is the vectorized version of connection=-[multisym(GA'*RA)*SA multisym(GB'*RB)*SB];
0139         connection = tangent2ambient(X,-cat(2,...
0140             multiprod(multisym(multiprod(multitransp(GA), RA)), SA),...
0141             multiprod(multisym(multiprod(multitransp(GB), RB)), SB)));
0142         rhess = M.proj(X,H + connection);
0143     end
0144
0145
0146
0147     M.exp = @exponential;
0148     function Y = exponential(X, U, t)
0149         if nargin == 3
0150             U = t*U;
0151         end
0152
0153         UFlat = essential_flat(U);
0154         exptUFlat = rot3_exp(UFlat);
0155         Y = essential_sharp(multiprod(essential_flat(X), exptUFlat));
0156     end
0157
0158     M.retr = @exponential;
0159
0160     M.log = @logarithm;
0161     function U = logarithm(X, Y)
0162
0163         QX = [X(:,1:3,:);X(:,4:6,:)];
0164         QY = [Y(:,1:3,:);Y(:,4:6,:)];
0165         QYr = essential_closestRepresentative(QX,QY,'flagSigned',flagSigned);
0166         Yr = [QYr(1:3,:,:) QYr(4:6,:,:)];
0167         U = zeros(size(X));
0168         U(:,1:3,:) = rot3_log(multiprod(multitransp(X(:,1:3,:)),Yr(:,1:3,:)));
0169         U(:,4:6,:) = rot3_log(multiprod(multitransp(X(:,4:6,:)),Yr(:,4:6,:)));
0170     end
0171
0172     M.hash = @(X) ['z' hashmd5(X(:))];
0173
0174     M.rand = @() randessential(k);
0175     function Q = randessential(N)
0176         % Generates random essential matrices.
0177         %
0178         % function Q = randessential(N)
0179         %
0180         % Q is a [3x6] matrix where each [3x3] block is a uniformly distributed
0181         % matrix.
0182
0183         if nargin < 1
0184             N = 1;
0185         end
0186
0187         Q = [randrot(3,N) randrot(3,N)];
0188     end
0189
0190     M.randvec = @randomvec;
0191     function U = randomvec(X)
0192         U = tangentProjection(X, essential_sharp(randskew(3, 2*k)));
0193         U = U / sqrt(M.inner([],U,U));
0194     end
0195
0196     M.lincomb = @matrixlincomb;
0197
0198     M.zerovec = @(x) zeros(3, 6, k);
0199
0200     M.transp = @transport;
0201     function S2 = transport(X1, X2, S1)
0202         % Transport a vector from the tangent space at X1 to the tangent
0203         % space at X2. This transport uses the left translation of the
0204         % ambient group and preserves the norm of S1. The left translation
0205         % aligns the vertical spaces at the two elements.
0206
0207         % Group operation in the ambient group, X12=X2'*X1
0208         X12 = essential_sharp(multiprod(multitransp(essential_flat(X2)),essential_flat(X1)));
0209         X12Flat = essential_flat(X12);
0210
0211         % Left translation, S2=X12*S*X12'
0212         S2 = essential_sharp(multiprod(X12Flat,multiprod(essential_flat(S1),multitransp(X12Flat))));
0213     end
0214
0215     M.pairmean = @pairmean;
0216     function Y = pairmean(X1, X2)
0217         V = M.log(X1, X2);
0218         Y = M.exp(X1, .5*V);
0219     end
0220
0221     M.dist = @(x, y) M.norm(x, M.log(x, y));
0222
0223     M.vec = @(x, u_mat) u_mat(:);
0224     M.mat = @(x, u_vec) reshape(u_vec, [3, 6, k]);
0225     M.vecmatareisometries = @() true;
0226
0227     M.tangent2ambient = @tangent2ambient;
0228
0229 end
0230
0231
0232 %% Some functions used by the essential factory
0233
0234 function v = vee3(V)
0235     v = squeeze([V(3,2,:)-V(2,3,:); V(1,3,:)-V(3,1,:); V(2,1,:)-V(1,2,:)])/2;
0236 end
0237
0238
0239 % Compute the exponential map in SO(3) using Rodrigues' formula
0240 %  function R = rot3_exp(V)
0241 % V must be a [3x3xN] array of [3x3] skew-symmetric matrices.
0242 function R = rot3_exp(V)
0243     v = vee3(V);
0244     nv = cnorm(v);
0245     idxZero = nv < 1e-15;
0246     nvMod = nv;
0247     nvMod(idxZero) = 1;
0248
0249     vNorm = v./([1;1;1]*nvMod);
0250
0251     % Matrix exponential using Rodrigues' formula
0252     nv = shiftdim(nv,-1);
0253     c = cos(nv);
0254     s = sin(nv);
0255     [VNorm,vNormShift] = essential_hat3(vNorm);
0256     vNormvNormT = multiprod(vNormShift,multitransp(vNormShift));
0257     R=multiprod(eye(3),c)+multiprod(VNorm,s)+multiprod(vNormvNormT,1-c);
0258 end
0259
0260
0261
0262 % Compute the logarithm map in SO(3)
0263 %  function V = rot3_log(R)
0264 % V is a [3x3xN] array of [3x3] skew-symmetric matrices
0265 function V = rot3_log(R)
0266     skewR = multiskew(R);
0267     ctheta = (multitrace(R)'-1)/2;
0268     stheta = cnorm(vee3(skewR));
0269     theta = atan2(stheta,ctheta);
0270
0271     V=skewR;
0272     for ik=1:size(R,3)
0273         V(:,:,ik)=V(:,:,ik)/sincN(theta(ik));
0274     end
0275 end
0276
0277
0278 function sx = sincN(x)
0279     sx = sin(x)./x;
0280     sx(x==0) = 1;
0281 end
0282
0283 function nv = cnorm(v)
0284     nv = sqrt(sum(v.^2));
0285 end
0286
0287 function M = vertproj(X, H)
0288     M = multiprod(X(3, 1:3, :), permute(vee3(H(:, 1:3, :)), [1 3 2])) + ...
0289         multiprod(X(3, 4:6, :), permute(vee3(H(:, 4:6, :)), [1 3 2]));
0290 end
0291
0292 function M = p1(X)
0293     M = X(:, 1:3, :);
0294 end
0295
0296 function M = p2(X)
0297     M = X(:, 4:6, :);
0298 end
0299
0300 function U = tangent2ambient(X, H)
0301     U = essential_sharp(multiprod(essential_flat(X), essential_flat(H)));
0302 end```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005