$ \DeclareMathOperator{\arccosh}{arccosh} \DeclareMathOperator*{\argmin}{arg\,min} \DeclareMathOperator{\Exp}{Exp} \newcommand{\geo}[2]{\gamma_{\overset{\frown}{#1,#2}}} \newcommand{\geoS}{\gamma} \newcommand{\geoD}[2]{\gamma_} \newcommand{\geoL}[2]{\gamma(#2; #1)} \newcommand{\gradM}{\nabla_{\M}} \newcommand{\gradMComp}[1]{\nabla_{\M,#1}} \newcommand{\Grid}{\mathcal G} \DeclareMathOperator{\Log}{Log} \newcommand{\M}{\mathcal M} \newcommand{\N}{\mathcal N} \newcommand{\mat}[1]{\mathbf{#1}} \DeclareMathOperator{\prox}{prox} \newcommand{\PT}[3]{\mathrm{PT}_{#1\to#2}#3} \newcommand{\R}{\mathbb R} \newcommand{\SPD}[1]{\mathcal{P}(#1)} \DeclareMathOperator{\Tr}{Tr} \newcommand{\tT}{\mathrm{T}} \newcommand{\vect}[1]{\mathbf{#1}} $

Illustrating the gradient of a second order difference

This example illustrates how to compute a gradient of a second order difference and how to vizualize the result using an export to Asymptote. If we are given three points on the manifold , the 2-sphere, the second order difference can be interpreted as (twice) the distance of the midpoint to , since in Euclidean space . The first figure illustrates this situation on the sphere for points on the equator and near the north pole, to be precise

1
2
3
4
5
M = Sn(2);
x = [1;0;0];
z = [0;1;0];
y = M.geopoint([0;0;1], M.geopoint(x,z,0.5),0.05);
c = M.midPoint(x,z);
Illustration of the mid point model (length of green curve) on the sphere \(\mathbb S^2\).
Illustration of the mid point model (length of green curve) on the sphere \(\mathbb S^2\).

The gradient with respect to of this second order absolute difference is derived in [1] and implemented within MVIRT in gradSecondOrderMidpointModel or can be computed using DxGeo and DyGeo of the manifold class as

1
2
3
4
5
% we take the inner differential (for x and z)
inner = -M.log(c,y)./M.dist(c,y);
gradX = M.AdjDxGeo(x,z,0.5,inner);
gradY = - M.log(y,c)./M.dist(c,y);
gradZ = M.AdjDyGeo(x,z,0.5,inner);

To illustrate the vectors and the points in Matlab, we can use quiver3 and plotS2 together as

1
2
3
4
5
6
7
8
9
10
11
12
pts = cat(2,x,y,z);
vecs = -cat(2,gradX,gradY,gradZ);
quiver3(pts(1,:),pts(2,:),pts(3,:),vecs(1,:),vecs(2,:),vecs(3,:));
vecs1 = permute(cat(3,pts,vecs),[1,3,2]);
hold on
plotS2(cat(2,x,y,z,c),'o','MarkerSize',3);
pts1 = cat(2,x,y,z,c);
plotS2(M.geodesic(x,z,'t',0:1/100:1),'-');
geo1 = M.geodesic(x,z,'t',0:1/100:1);
plotS2(M.geodesic(y,c,'t',0:1/100:1),'--b');
geo2 = M.geodesic(y,c,'t',0:1/100:1);
hold off

plotting the decent direction tangent vectors , all 4 points, and two geodesics. The corresponding image in Asymptote is shown in the second Figure

Illustration of the gradients at \(x,y,z\) of \(d_2\).
Illustration of the gradients at \(x,y,z\) of \(d_2\).

We perform a step towards the negative gradient direction with step size , i.,e.

1
2
3
  xN = M.exp(x,-gradX);
  yN = M.exp(y,-gradY);
  zN = M.exp(z,-gradZ);

yields finally the third figure

A step towards \(-\gradM d_2\).
A step towards \(-\gradM d_2\).

and a reduction of the second order absolute difference from about (old dark green segment) to (light green segment).

Creating the additional arrays of points

1
2
3
pts2 = cat(2,xN,yN,zN,M.midPoint(xN,zN));
geo3 = M.geodesic(xN,zN,'t',0:1/100:1);
geo4 = M.geodesic(yN,M.midPoint(xN,zN),'t',0:1/100:1);

we can use the exportSpherePCT2Asy export to export points, curves (connected points) and tangent vectors to Asymptote with

1
2
3
4
5
6
7
exportSpherePCT2Asy({pts1,pts2},{geo1,geo2,geo3,geo4},{vecs1},
{ [0;0;.66],[.33;0;.66],...
  [.6;.6;.6],[0;0.5;0],[0.5;.6;.5],[0;1;0],...
  [0;0.5;1]},
  'File','gradTV2onS2/grad.asy',...
  'DotSize',[5,5,1,1,1,1.2,2],...
  'Camera',[1,.8,.6],'ArrowHead',10);

where the second and third line specify the colors of all 7 sets.

See also

References

  1. Bačák, M, Bergmann, R, Steidl, G and Weinmann, A (2016). A second order non-smooth variational model for restoring manifold-valued images. SIAM Journal on Scientific Computing. 38 A567–A597