14
$\begingroup$

Define a circle $\cal{C}(p,r)$ on the surface of an ellipsoid $E$ in $\mathbb{R}^3$ to be the set of points on $E$ whose shortest geodesic distance from centerpoint $p$ is $r$. Let me assume that $r$ is small enough so that $\cal{C}(p,r)$ does not self-intersect.

Given an ellipsoid centered on the origin, with its three axes $a,b,c$ aligned along $x,y,z$, I would like to draw $\cal{C}(p,r)$, where $p$ is given either in $(x,y,z)$ coordinates, or maybe via a unit vector $u$ from the origin, $p = u t$ for a scaling $t$.

The only way I see how to accomplish this is to numerically track geodesics from $p$, shooting off geodesics at many angles around $p$, stopping each when its length reaches $r$. It is a project to implement this.

Is there some easier route, perhaps employing Regions? I am not seeking code, but just an idea for an easier approach.


Added 21May2019. Much earlier I did the opposite: Ellipses on spheres (and other surfaces) :
         EllipseOnSphere1


$\endgroup$
  • 1
    $\begingroup$ Have you not tried anything though? $\endgroup$ – MarcoB Apr 24 '18 at 1:14
  • $\begingroup$ @MarcoB: No. I added a clarification to the post to explain. Thanks. $\endgroup$ – Joseph O'Rourke Apr 24 '18 at 11:23
  • 2
    $\begingroup$ Do you know how to compute a geodesic distance on an ellipsoid? If not maybe you should pose a question on that first. $\endgroup$ – george2079 Apr 24 '18 at 15:44
  • $\begingroup$ @george2079: Good point. I cannot address this now, but will eventually. Thanks. $\endgroup$ – Joseph O'Rourke Apr 24 '18 at 19:28
  • 2
    $\begingroup$ You could use ParametricNDSolve, with the differential equation for a geodesic on an ellipsoid, and supply the starting direction of the geodesic as parameter. $\endgroup$ – Niki Estner Apr 25 '18 at 8:50
11
$\begingroup$

A quick attempt can be made on discretised surfaces, where one can use the Geodesics in Heat Algorithm to calculate the geodesic distance between a vertex i and all other vertices on a surface. This algorithm is implemented here. We first create a discretised surface of an ellipsoid:

ellipsoid = BoundaryDiscretizeRegion[ImplicitRegion[x^2/4 + y^2/9 + z^2 <= 1,{x, y, z}], MaxCellMeasure -> {"Length" -> 0.2}];

Then using the algorithm:

prep = heatDistprep[b];
sol = solveHeat[b, prep, 5, 0.1];

which when plotted:

plotHeat[plotdata_, a_, dist_] := 
ListSliceContourPlot3D[plotdata, a, ContourShading -> Automatic, 
ColorFunction -> "BrightBands", Boxed -> False, Axes -> False, 
BoxRatios -> Automatic, Contours -> {dist}]
plotHeat[sol[[1]], b, -2]

geodesic circle on ellipsoid

This requires one to choose which vertex is the centre of the ellipse (in this case vertex 5) the geodesic "radius" (in this case -2) needs to be rescaled to be converted to real distances. Of course to make this quick answer better one would need to calculate firstly which vertex corresponds to the starting point on the ellipse, and then rescale the radius to real distances. I hope to do this if I get time.

For completeness the code of the geodesics in heat algorithm is given here:

heatDistprep[mesh0_] := Module[{a = mesh0, vertices, nvertices, edges, edgelengths, nedges, faces, faceareas, unnormfacenormals, acalc, facesnormals, facecenters, nfaces, oppedgevect, wi1, wi2, wi3, sumAr1, sumAr2, sumAr3, areaar, gradmat1, gradmat2, gradmat3, gradOp, arear2, divMat, divOp, Delta, t1, t2, t3, t4, t5, , Ac, ct, wc, deltacot, vertexcoordtrips, adjMat},
vertices = MeshCoordinates[a]; (*List of vertices*)
edges = MeshCells[a, 1] /. Line[p_] :> p; (*List of edges*)
faces = MeshCells[a, 2] /. Polygon[p_] :> p; (*List of faces*)
nvertices = Length[vertices];
nedges = Length[edges];
nfaces = Length[faces];
adjMat = SparseArray[Join[({#1, #2} -> 1) & @@@ edges, ({#2, #1} -> 1) & @@@edges]]; (*Adjacency Matrix for vertices*)
edgelengths = PropertyValue[{a, 1}, MeshCellMeasure];
faceareas = PropertyValue[{a, 2}, MeshCellMeasure];
vertexcoordtrips = Map[vertices[[#]] &, faces];
unnormfacenormals = Cross[#3 - #2, #1 - #2] & @@@ vertexcoordtrips;
acalc = (Norm /@ unnormfacenormals)/2;
facesnormals = Normalize /@ unnormfacenormals;
facecenters = Total[{#1, #2, #3}]/3 & @@@ vertexcoordtrips;
oppedgevect = (#1 - #2) & @@@ Partition[#, 2, 1, 3] & /@vertexcoordtrips;
wi1 = -Cross[oppedgevect[[#, 1]], facesnormals[[#]]] & /@Range[nfaces];
wi2 = -Cross[oppedgevect[[#, 2]], facesnormals[[#]]] & /@Range[nfaces];
wi3 = -Cross[oppedgevect[[#, 3]], facesnormals[[#]]] & /@Range[nfaces];
sumAr1 = SparseArray[Join[Map[{#, faces[[#, 1]]} -> wi1[[#, 1]] &, Range[nfaces]],Map[{#, faces[[#, 2]]} -> wi2[[#, 1]] &, Range[nfaces]],Map[{#, faces[[#, 3]]} -> wi3[[#, 1]] &, Range[nfaces]]]];
sumAr2 = SparseArray[Join[Map[{#, faces[[#, 1]]} -> wi1[[#, 2]] &, Range[nfaces]], Map[{#, faces[[#, 2]]} -> wi2[[#, 2]] &, Range[nfaces]],Map[{#, faces[[#, 3]]} -> wi3[[#, 2]] &, Range[nfaces]]]];
sumAr3 =SparseArray[Join[Map[{#, faces[[#, 1]]} -> wi1[[#, 3]] &, Range[nfaces]], Map[{#, faces[[#, 2]]} -> wi2[[#, 3]] &, Range[nfaces]], Map[{#, faces[[#, 3]]} -> wi3[[#, 3]] &, Range[nfaces]]]];
areaar = SparseArray[Table[{i, i} -> 1/(2*acalc[[i]]), {i, nfaces}]];
gradmat1 = areaar.sumAr1;
gradmat2 = areaar.sumAr2;
gradmat3 = areaar.sumAr3;
gradOp[u_] := Transpose[{gradmat1.u, gradmat2.u, gradmat3.u}];
arear2 = SparseArray[Table[{i, i} -> (2*faceareas[[i]]), {i, nfaces}]];
divMat = {Transpose[gradmat1].arear2, Transpose[gradmat2].arear2,Transpose[gradmat3].arear2};
divOp[q_] := divMat[[1]].q[[All, 1]] + divMat[[2]].q[[All, 2]] + divMat[[3]].q[[All, 3]];
Delta = divMat[[1]].gradmat1 + divMat[[2]].gradmat2 + divMat[[3]].gradmat3;
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}]; (*Required to allow addition of value assignment to Sparse Array*)
t1 = Join[faces[[All, 1]], faces[[All, 2]], faces[[All, 3]]];
t2 = Join[acalc, acalc, acalc];
Ac = SparseArray[Table[{t1[[i]], t1[[i]]} -> t2[[i]], {i, nfaces*3}]];
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];
{Ac, Delta, gradOp, divOp, nvertices, vertices, adjMat}
]

solveHeat[mesh0_, prepvals_, i0_, t0_] := Module[{nvertices, delta, t, u, Ac, Delta, g, h, phi, gradOp, divOp, vertices, plotdata},
vertices = prepvals[[6]];
nvertices = prepvals[[5]];
Ac = prepvals[[1]];
Delta = prepvals[[2]];
gradOp = prepvals[[3]];
divOp = prepvals[[4]];
delta = Table[If[i == i0, 1, 0], {i, nvertices}];
t = t0;
u = LinearSolve[(Ac + t*Delta), delta];
g = gradOp[u];
h = -Normalize /@ g;
phi = LinearSolve[Delta, divOp[h]];
plotdata = Map[Join[vertices[[#]], {phi[[#]]}] &, Range[Length[vertices]]];

{phi}
] 
$\endgroup$
  • $\begingroup$ Wonderful! Could you try to update the link to the Heat Algorithm? Seems broken... $\endgroup$ – Joseph O'Rourke Jun 18 '18 at 19:27
  • 1
    $\begingroup$ The link to the paper seems to still work (although it takes a long time to download). $\endgroup$ – Dunlop Jun 18 '18 at 21:29
  • $\begingroup$ Thank you, got it! $\endgroup$ – Joseph O'Rourke Jun 18 '18 at 22:40
  • $\begingroup$ @Dunlop That's because I updated the link. ;) $\endgroup$ – Henrik Schumacher Jun 18 '18 at 23:07
  • $\begingroup$ @HenrikSchumacher many thanks! $\endgroup$ – Dunlop Jun 19 '18 at 3:28
8
$\begingroup$

Long, long time ago, I wrote a geodesic shooter for triangle surfaces. I took the opportunity and refined it a bit. Maybe someone wants to play with it.

Clearly, this cannot compete with the heat method when it comes to speed; the heat method computes all geodesic distances from a point with only two sparse linear solves (one for the heat kernel and one for the Hodge decomposition). It was also never meant to compete; the original application was for updating curves that are constraint to a given surface. It is also supposed to be able to perform parallel transport for a set of vectors (to be specified by the "TransportedVectors" option). However, I haven't tested this feature, yet.

Note that you will need IGraphM installed for this to work.

Options[ShootGeodesic] = {
   "MaxIterations" -> 1000000,
   "TransportedVectors" -> {},
   "GeodesicData" -> Automatic
   };

ShootGeodesic[R_MeshRegion, p0_, u0_, OptionsPattern[]] :=   
 Block[{pts, faces, facenormals, p, pbag, vbag, ff, face, ν, u, P, distance, iter, bool, b, ee, edge, t, νnew, unew, ffnew, rot, maxiter, data, edgelookuptable, A12, v, transportQ},
  pts = MeshCoordinates[R];
  facenormals = Region`Mesh`MeshCellNormals[R, 2];
  faces = MeshCells[R, 2, "Multicells" -> True][[1, 1]];

  data = OptionValue["GeodesicData"];
  If[Head[data] =!= Association,
   data = GeodesyData[R];
   ];
  edgelookuptable = data["EdgeLookupTable"];
  A12 = data["EdgeFaceAdjacencyMatrix"];
  v = OptionValue["TransportedVectors"];
  transportQ = Length[v] > 0 && Dimensions[v][[2]] == 3;
  vbag = Internal`Bag[{v}];

  maxiter = OptionValue["MaxIterations"];

  ff = Region`Mesh`MeshNearestCellIndex[R, p0][[2]];
  p = RegionNearest[R, p0];
  pbag = Internal`Bag[{p}];
  face = faces[[ff]];
  ν = facenormals[[ff]];
  u = u0 - ν ν.u0;
  distance = Norm[u];
  u = u/distance;
  P = pts[[face]];
  iter = 0;
  bool = True;
  While[bool && iter < maxiter,
   iter++;
   {t, edge} = getGeodesicsols[p, u, P];
   If[t < distance,
    distance -= t;
    p = p + t u;
    Internal`StuffBag[pbag, p];
    ee = edgelookuptable[[Sequence @@ Switch[Round[edge],
        1, face[[{2, 3}]],
        2, face[[{3, 1}]],
        3, face[[{1, 2}]]
        ]]];
    ff = Complement[A12[[ee]]["AdjacencyLists"], {ff}][[1]];
    νnew = facenormals[[ff]];
    If[
     ν.νnew < 1.,
     rot = MyRotationMatrix[{ν, νnew}];
     u = rot.u;
     If[transportQ,
      v = v.Transpose[rot];
      Internal`StuffBag[vbag, v];
      ];
     ,
     u = u;
     If[transportQ,
      v = v.Transpose[rot];
      Internal`StuffBag[vbag, v];
      ];
     ];
    u = Normalize[u - νnew νnew.u];
    ν = νnew;
    face = faces[[ff]];
    P = pts[[face]];
    ,
    p = p + distance u;
    Internal`StuffBag[pbag, p];
    If[transportQ, Internal`StuffBag[vbag, v]];
    bool = False;
    ];
   ];
  If[iter == maxiter, 
   Print["Warning: MaxIterations ", maxiter, " reached!"]];
  Association[
   "Point" -> p,
   "DirectionVector" -> distance u,
   "TransportedVectors" -> Internal`BagPart[vbag, All],
   "Face" -> ff,
   "Trajectory" -> Internal`BagPart[pbag, All]
   ]
  ];

(* The working horse that handles the intersection of a geodesic with the triangle boundaries. *)
Quiet[
  Block[{YY, VV, XX, UU, PP, Y, V, X, U, P, s, t, A},
    PP = Table[Compile`GetElement[P, i, j], {i, 1, 3}, {j, 1, 3}];
    XX = Table[Compile`GetElement[X, i], {i, 1, 3}];
    UU = Table[Compile`GetElement[U, i], {i, 1, 3}];
    YY = Table[Compile`GetElement[Y, i], {i, 1, 2}];
    VV = Table[Compile`GetElement[V, i], {i, 1, 2}];
    A = Transpose[{PP[[2]] - PP[[1]], PP[[3]] - PP[[1]]}];
    With[{
      ϵ = 1. 10^(-14),
      sol1 = Inverse[Transpose[{{-1, 1}, -VV}]].(YY - {1, 0}),
      sol2 = Inverse[Transpose[{{0, -1}, -VV}]].(YY - {0, 1}),
      sol3 = Inverse[Transpose[{{1, 0}, -VV}]].YY,
      Adagger = (Inverse[A\[Transpose].A].A\[Transpose])
      },
     getGeodesicsols = Compile[{{X, _Real, 1}, {U, _Real, 1}, {P, _Real, 2}},
        Block[{V, Y, edge, Bag, sols, pos, tvals},
         Y = Adagger.(X - P[[1]]);
         V = Adagger.U;
         sols = {
           If[Abs[Compile`GetElement[V, 1] + Compile`GetElement[V, 2]] <= ϵ, {2., 0.}, sol1],
           If[Abs[Compile`GetElement[V, 1]] <= ϵ, {2., 0.}, sol2],
           If[Abs[Compile`GetElement[V, 2]] <= ϵ, {2., 0.}, sol3]
           };
         Bag = Internal`Bag[Most[{0}]];
         Do[
          If[-ϵ <= sols[[i, 1]] <= 1. + ϵ && -ϵ <= sols[[i, 2]], 
           Internal`StuffBag[Bag, i, 1]],
          {i, 1, 3}
          ];
         pos = Internal`BagPart[Bag, All];
         tvals = sols[[All, 2]];
         edge = First@pos[[Ordering[tvals[[pos]], -1]]];
         {tvals[[edge]], N[edge]}
         ]
        ];
     ];
    ];
  ];

(* Quick way to compute rotation matrices *)
Block[{angle, v, vv, u, uu, ww, e1, e2, e2prime, e3},
  uu = Table[Compile`GetElement[u, i], {i, 1, 3}];
  vv = Table[Compile`GetElement[v, i], {i, 1, 3}];
  ww = Cross[uu, vv];
  e2 = Cross[ww, uu];
  e2prime = Cross[ww, vv];
  With[{code = N[Plus[
       KroneckerProduct[vv, uu]/Sqrt[uu.uu]/Sqrt[vv.vv], 
       KroneckerProduct[e2prime, e2]/Sqrt[e2.e2]/Sqrt[e2prime.e2prime],
       KroneckerProduct[ww, ww]/ww.ww]
      ]
    },
   rotationMatrix3DVectorVector = Compile[{{u, _Real, 1}, {v, _Real, 1}}, code,
     CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
     Parallelization -> True, RuntimeOptions -> "Speed"
     ]
   ];
  ];

MyRotationMatrix[{u_, v_}] := rotationMatrix3DVectorVector[u, v];

The last two functions are copied from How to speed up RotationMatrix?

I am really fond of precomputing recycable data. So here a generator for some useful combinatorics.

Needs["IGraphM`"];
GeodesicData[R_MeshRegion] := (
   Association[
    "EdgeFaceAdjacencyMatrix" -> IGMeshCellAdjacencyMatrix[R, 1, 2],
    "EdgeLookupTable" -> 
     With[{edges = MeshCells[R, 1, "Multicells" -> True][[1, 1]]},
      SparseArray[
       Rule[
        Join[edges, Transpose@Reverse@Transpose@edges],
        Join[Range[Length[edges]], Range[Length[edges]]]
        ],
       {1, 1} Length[edges]
   ]
  ]
 ]
);

Application

Let's create a discrete ellipsoid, precompute the GeodesicData; pick a random point and a random direction; and compute a long geodesic.

R = RegionBoundary@
   BoundaryDiscretizeRegion[Ellipsoid[{0, 0, 0}, {3, 4, 2}], MaxCellMeasure -> 0.01];
data = GeodesicData[R];
SeedRandom[123];
p0 = RegionNearest[R, RandomPoint[R]];
u0 = RandomReal[{10, 1000}] RandomPoint[Sphere[]];
result = ShootGeodesic[R, p0, u0, "GeodesicData" -> data];

Show[
  R, 
  Graphics3D[{Specularity[White, 30], 
   Sphere[p0, 0.1], Gray, 
   Tube[result[["Trajectory"]], 0.01]}
  ]
 ]

enter image description here

And this is how we can draw a geodesic circle:

ff = Region`Mesh`MeshNearestCellIndex[R, p0][[2]];
ν = Region`Mesh`MeshCellNormals[R, {2, ff}];
{e1, e2} = Orthogonalize[Join[{ν}, N[IdentityMatrix[3][[Ordering[Abs[ν], 2]]]]]][[ 2 ;; 3]];
r = 3;
circle = ShootGeodesic[R, p0, r (Cos[#] e1 + Sin[#] e2), "GeodesicData" -> data
     ] & /@ Most@Subdivide[0., 2 Pi, 72];

Show[
 R, 
 Graphics3D[{Specularity[White, 30], 
  Sphere[p0, 0.1], 
  Gray, Tube[Join[#, {#[[1]]}], 0.035] &[circle[[All, "Point"]]], 
  Lighter@Lighter@Gray, Tube[{#}, 0.01] & /@ circle[[1 ;; -1 ;; 2, "Trajectory"]]}
  ]
 ]

enter image description here

Of course, we can draw geodesic circles also onto other surfaces:

R = ExampleData[{"Geometry3D", "Triceratops"}, "MeshRegion"];
data = GeodesicData[R];
SeedRandom[1234];
p0 = RegionNearest[R, RandomPoint[R]];
ff = Region`Mesh`MeshNearestCellIndex[R, p0][[2]];
ν = Region`Mesh`MeshCellNormals[R, {2, ff}];
{e1, e2} = Orthogonalize[Join[{ν}, N[IdentityMatrix[3][[Ordering[Abs[ν], 2]]]]]][[2 ;; 3]];
r = 1;
circles = Table[
   ShootGeodesic[R, p0, r (Cos[#] e1 + Sin[#] e2), "GeodesicData" -> data
      ] & /@ Most@Subdivide[0., 2 Pi, 180],
   {r, 0.2, 2, 0.2}
   ];

Show[
 R,
 Graphics3D[{
   Specularity[White, 30],
   Sphere[p0, 0.05],
   EdgeForm[], Gray,
   Tube[Join[#, {#[[1]]}], 0.02] & /@ circles[[All, All, "Point"]]
   }]
 ]

enter image description here

Remark

I do not catch boundary collisions, so this is only guaranteed to work for triangle meshes with the topology of a closed surface.

$\endgroup$
  • $\begingroup$ Beautiful! Despite the speed advantage of the heat algorithm, tracing long geodesics could prove very useful to me. $\endgroup$ – Joseph O'Rourke Jun 18 '18 at 22:43
  • $\begingroup$ @JosephO'Rourke Thanks! See also here for a more efficient implementation of the heat method (basically, the precomputation of the FEM matrices got sped up by several magnitudes). $\endgroup$ – Henrik Schumacher Jun 18 '18 at 23:09
  • $\begingroup$ This is great, I have made a link to your implementation here $\endgroup$ – Dunlop Jun 21 '18 at 6:45
3
$\begingroup$

This is not an answer but an example of what you can do. (Hopefully close enough to what you want.)

{a, b, c} = {3, 4, 2};
{r, p} = {3, {0, 1/2, 1/2}}; 

elps = Graphics3D[Ellipsoid[{0, 0, 0}, {a, b, c}]];

sph = Graphics3D[Ellipsoid[p, {r, r, r}]];

ints = Region[
 ImplicitRegion[
  x^2/a^2 + y^2/b^2 + z^2 /c^2 == 1 && 
  Norm[{x, y, z} - p, 2]^2 == r^2, {x, y, z}], Boxed -> True]

enter image description here

Show[elps, ints]

enter image description here

$\endgroup$
  • $\begingroup$ That's very nice! But I think the intersection of a sphere and an ellipse does not yield the surface circle that I seek. $\endgroup$ – Joseph O'Rourke Apr 24 '18 at 11:23
  • $\begingroup$ @JosephO'Rourke Right, you should write the right equation in ImplicitRegion. $\endgroup$ – Anton Antonov Apr 24 '18 at 14:05
  • $\begingroup$ Thanks, Anton. I will investigate ImplicitRegion. $\endgroup$ – Joseph O'Rourke Apr 24 '18 at 14:43
2
$\begingroup$

Your initial hunch was correct. You will need to numerically solve geodesics and approximate the exponential map. A region in Mathematica is essentially defined as a subset of $E^n$ with the Lebesgue measure. (Points are exceptional that RegionMeaure[pt] uses couting measure.) Thus, there is no simple representation of an object that is only defined with respect to the induced metric, which is local rather than global property.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.