# Draw circle on ellipsoid

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) :

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

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_] :=
ColorFunction -> "BrightBands", Boxed -> False, Axes -> False,
BoxRatios -> Automatic, Contours -> {dist}]
plotHeat[sol[[1]], b, -2]


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}]];
arear2 = SparseArray[Table[{i, i} -> (2*faceareas[[i]]), {i, nfaces}]];
divOp[q_] := divMat[[1]].q[[All, 1]] + divMat[[2]].q[[All, 2]] + divMat[[3]].q[[All, 3]];
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}];
]

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]];
divOp = prepvals[[4]];
delta = Table[If[i == i0, 1, 0], {i, nvertices}];
t = t0;
u = LinearSolve[(Ac + t*Delta), delta];
h = -Normalize /@ g;
phi = LinearSolve[Delta, divOp[h]];
plotdata = Map[Join[vertices[[#]], {phi[[#]]}] &, Range[Length[vertices]]];

{phi}
]

• Wonderful! Could you try to update the link to the Heat Algorithm? Seems broken... Jun 18, 2018 at 19:27
• The link to the paper seems to still work (although it takes a long time to download). Jun 18, 2018 at 21:29
• Thank you, got it! Jun 18, 2018 at 22:40
• @Dunlop That's because I updated the link. ;) Jun 18, 2018 at 23:07
• @HenrikSchumacher many thanks! Jun 19, 2018 at 3:28

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 = RegionMeshMeshCellNormals[R, 2];
faces = MeshCells[R, 2, "Multicells" -> True][[1, 1]];

data = OptionValue["GeodesicData"];
data = GeodesyData[R];
];
edgelookuptable = data["EdgeLookupTable"];
v = OptionValue["TransportedVectors"];
transportQ = Length[v] > 0 && Dimensions[v][[2]] == 3;
vbag = InternalBag[{v}];

maxiter = OptionValue["MaxIterations"];

ff = RegionMeshMeshNearestCellIndex[R, p0][[2]];
p = RegionNearest[R, p0];
pbag = InternalBag[{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;
InternalStuffBag[pbag, p];
ee = edgelookuptable[[Sequence @@ Switch[Round[edge],
1, face[[{2, 3}]],
2, face[[{3, 1}]],
3, face[[{1, 2}]]
]]];
νnew = facenormals[[ff]];
If[
ν.νnew < 1.,
rot = MyRotationMatrix[{ν, νnew}];
u = rot.u;
If[transportQ,
v = v.Transpose[rot];
InternalStuffBag[vbag, v];
];
,
u = u;
If[transportQ,
v = v.Transpose[rot];
InternalStuffBag[vbag, v];
];
];
u = Normalize[u - νnew νnew.u];
ν = νnew;
face = faces[[ff]];
P = pts[[face]];
,
p = p + distance u;
InternalStuffBag[pbag, p];
If[transportQ, InternalStuffBag[vbag, v]];
bool = False;
];
];
If[iter == maxiter,
Print["Warning: MaxIterations ", maxiter, " reached!"]];
Association[
"Point" -> p,
"DirectionVector" -> distance u,
"TransportedVectors" -> InternalBagPart[vbag, All],
"Face" -> ff,
"Trajectory" -> InternalBagPart[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[CompileGetElement[P, i, j], {i, 1, 3}, {j, 1, 3}];
XX = Table[CompileGetElement[X, i], {i, 1, 3}];
UU = Table[CompileGetElement[U, i], {i, 1, 3}];
YY = Table[CompileGetElement[Y, i], {i, 1, 2}];
VV = Table[CompileGetElement[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,
},
getGeodesicsols = Compile[{{X, _Real, 1}, {U, _Real, 1}, {P, _Real, 2}},
Block[{V, Y, edge, Bag, sols, pos, tvals},
sols = {
If[Abs[CompileGetElement[V, 1] + CompileGetElement[V, 2]] <= ϵ, {2., 0.}, sol1],
If[Abs[CompileGetElement[V, 1]] <= ϵ, {2., 0.}, sol2],
If[Abs[CompileGetElement[V, 2]] <= ϵ, {2., 0.}, sol3]
};
Bag = InternalBag[Most[{0}]];
Do[
If[-ϵ <= sols[[i, 1]] <= 1. + ϵ && -ϵ <= sols[[i, 2]],
InternalStuffBag[Bag, i, 1]],
{i, 1, 3}
];
pos = InternalBagPart[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[CompileGetElement[u, i], {i, 1, 3}];
vv = Table[CompileGetElement[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[
"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]}
]
]


And this is how we can draw a geodesic circle:

ff = RegionMeshMeshNearestCellIndex[R, p0][[2]];
ν = RegionMeshMeshCellNormals[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"]]}
]
]


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 = RegionMeshMeshNearestCellIndex[R, p0][[2]];
ν = RegionMeshMeshCellNormals[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"]]
}]
]


# Remark

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

• Beautiful! Despite the speed advantage of the heat algorithm, tracing long geodesics could prove very useful to me. Jun 18, 2018 at 22:43
• @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). Jun 18, 2018 at 23:09
• This is great, I have made a link to your implementation here Jun 21, 2018 at 6:45

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]


Show[elps, ints]


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

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.