# Fast spherical polygon

Cross posted on Wolfram Community

### Question

Given a list of points on a sphere and the sphere's radius, I'd like to plot a spherical polygon with those points as vertices.

And this needs to be fast, fast enough for the user to not "feel" generation time.

One should be able to style them too. Most importantly the surface, but an edge style would be nice as well.

### Motivation

I think it will be useful in many applications.

I don't have time for this but I thought it would be a nice feature to have to improve code I was playing with lately, mostly based on another J.M. answer - Voronoi grid on a sphere

arc[center_?VectorQ, {start_?VectorQ, end_?VectorQ}] :=  Module[{ang, co, r},
ang = VectorAngle[start - center, end - center]; co = Cos[ang/2]; r = EuclideanDistance[center, start]; {{start, center + r/co Normalize[(start + end)/2 - center], end}, co}
]
points = {2 π #1, ArcCos[2 #2 - 1]} & @@@ RandomReal[1, {10, 2}];
sp = Append[Sin[#2] Through[{Cos, Sin}[#1]], Cos[#2]] & @@@ points;
proc[] := (
ch = ConvexHullMesh[sp]; verts = MeshCoordinates[ch]; polys = First /@ MeshCells[ch, 2]; voro = Normalize[ Cross[verts[[#2]] - verts[[#1]],  verts[[#3]] - verts[[#1]]]] & @@@ polys; edges =  arc[{0, 0, 0}, voro[[##]]] & /@ Select[Subsets[Range[Length[polys]], {2}], Length[Intersection @@ polys[[#]]] >= 2 &];
);

proc[];

DynamicModule[{run = True},  Graphics3D[{ {Opacity[.75],
DynamicWrapper[EventHandler[Sphere[],
"MouseMoved" :>  Module[{pos = MousePosition["Graphics3DBoxIntercepts", True], pt},  If[ Not @ TrueQ @ pos , pt = RegionIntersection[Sphere[], Line @ pos]; If[pt =!= EmptyRegion,  sp[[-1]] = First@Nearest[pt[], pos[]]; proc[]]  ]]]
, TrackedSymbols :> {run}
]
}
, {AbsoluteThickness, Dynamic[BSplineCurve[#, SplineDegree -> 2,  SplineKnots -> {0, 0, 0, 1, 1, 1}, SplineWeights -> {1, #2, 1}] & @@@ edges]}, {Red, Sphere[Most @ sp, .02], Dynamic @ Sphere[Last @ sp, .02]}
}
, PlotRange -> 1.1    , SphericalRegion -> True   , ImageSize -> 500]
] To look more like this: • I had to upvote, of course. Yes, Voronoi, and Lloyd as well, are the reasons why I'm trying to research on quickly rendering spherical polygons. Apr 21, 2017 at 17:02
• Is it really possible that the graphics hardware supports fewer clipping planes than you might possibly want? Apr 22, 2017 at 1:47
• @IgorRivin Documentation rarely is superfluous so I think there is a good reason why this note is there. Don't know so much about technical aspects but one can always ask the support.
– Kuba
Apr 22, 2017 at 19:39
• @Kuba I cannot imagine any modern graphics hardware unable to deal with fewer than hundreds of thousands clipping planes. As for why it is in the documentation, that is a fair point. Apr 22, 2017 at 20:47

This might or might not be helpful. I want to show it nevertheless. If you have a convex polygon or a non-convex polygon that is already triangulated, you could subdivide the polygon by introducing new points that you push onto the sphere's surface.

For convex polygons (that don't span more than a hemisphere) a simple initial triangulation can be obtained by creating a new point that is in the centre of all vertices. Starting from this, you can subdivide each triangle until you have a mesh that is narrow enough to give a good visual result.

Dividing a triangle can be done in several ways. To get regular mesh, you could look at all three edges and introduce a new point onto the longest edge. Then you use this new point and replace the old triangle by 2 new ones. This is how the replacement looks. Note that the newly introduced bottom point is already projected onto the sphere. We can use a maxLen (or delta, appreciate my consistent naming) to stop the iteration at a certain level of detail. A compiled function subdivideTriC that does the subdivision together with a highlevel subdivideTriangles that iterates everything is given here

subdivideTriC = Compile[{{pts, _Real, 2}, {maxLen, _Real, 0}},
Module[{maxEdge = {1, 2}, vertex = 3, max = 0., tmp = 0.,
newPoint = {0., 0.}},
Do[
tmp = Norm[pts[[Mod[i, 3] + 1]] - pts[[i]]];
If[tmp > max,
max = tmp;
maxEdge = {i, Mod[i, 3] + 1};
vertex = Mod[i + 1, 3] + 1;
], {i, 3}
];
If[max < maxLen,
{pts},
newPoint = pts[[maxEdge[]]] + pts[[maxEdge[]]];
newPoint = Norm[pts[]]*newPoint/Norm[newPoint];
{{pts[[vertex]], pts[[maxEdge[]]], newPoint}, {pts[[vertex]],
newPoint, pts[[maxEdge[]]]}}
]
], Parallelization -> True, RuntimeAttributes -> {Listable}
];

subdivideTriangles[polys_, delta_] :=
Partition[
Partition[
Flatten@NestWhile[subdivideTriC[#, delta] &, polys,
Flatten[#1] != Flatten[#2] &, 2], 3], 3];


A subdivision of 56 triangles from your example down to a delta of 0.2 takes about 0.07 seconds on my machine. To show all those polygons from your example, I use your verts and polys and I can run

Graphics3D[
With[{pts = subdivideTriangles[verts[[#]], .2]}, {EdgeForm[],
RandomColor[], Polygon[pts, VertexNormals -> pts]}] & /@ polys] This method has one drawback: Since all triangles are subdivided in parallel, I cannot reuse existing vertices which creates a lot of additional points.

A better implementation would use a structure like GraphicsComplex uses, where each vertex is contained only one time in the surface and the polygons only refer to the id of the vertex. Unfortunately, at this time it is not clear to me how I can parallelize such a subdivision.

Another drawback is that I have to handle two datasets: the list of points and the list of polygon ids. Without tricks, I cannot return such a data-structure from a compiled function.

Nevertheless, let me share the approach. Here, I select the largest edge of all triangles and refine this. In the process, the two triangles that share the edge are subdivided:

subdivideTriangles2[{p_List, polys : {{_Integer ..} ..}}] := Module[{
r = Norm[p[]],
len = Length[p],
largest = {1, 2},
maxDist = 0.0,
p1 = 1,
p2 = 1,
newPoint, b = InternalBag[Most[{0}]]
},
Do[
p1 = polys[[poly, Mod[i, 3] + 1]];
p2 = polys[[poly, Mod[i + 1, 3] + 1]];
If[Norm[p[[p2]] - p[[p1]]] > maxDist,
maxDist = Norm[p[[p2]] - p[[p1]]];
largest = {p1, p2};
],
{i, 0, 3}, {poly, Length[polys]}
];
newPoint = p[[largest[]]] + p[[largest[]]];
newPoint = r*(newPoint)/Norm[newPoint];
Do[
If[
Length[Intersection[poly, largest]] === 2,
p1 = First[Complement[poly, largest]];
InternalStuffBag[
b, {{p1, largest[], len + 1}, {p1, len + 1, largest[]}}, 2],
InternalStuffBag[b, poly, 1]
]
, {poly, polys}
];
{Append[p, newPoint], Partition[InternalBagPart[b, All], 3]}
]

pts = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
Manipulate[
Graphics3D[{Opacity[.4], Sphere[{0, 0, 0}, .95], Red,
GraphicsComplex[#1, Polygon[#2]] & @@
Nest[subdivideTriangles2, {N[pts], {{1, 2, 3}}}, i]
}], {i, 0, 150, 1}] Finally, I was looking through the literature and one publication caught my eye: A marching method for the triangulation of surfaces. There it is described how you can triangulate along a surface iteratively. Since a sphere can be regared as implicit surface as well the approach would fit. It iteratively expands a wave front of polygons which would work in our case to fill a (even non-convex) polygon on the sphere. Unfortunately, I'm very pessimistic that this would even half way meet the speed expectations since you have to track vertex angles etc. and the method creates only one new polygon at a time as well.

## Automatic ClipPlanes approach

It is obviously limited to convex polygons, but should be useful anyway.

pts ought to be in order, like in a Polygon.

### Definition

sphericalPolygon // ClearAll

sphericalPolygon[
pts_List, center_List: {0, 0, 0}, radius : _?Positive : 1
] :=  Style[
, ClipPlanes -> Table[With[
{ ord = If[
VectorAngle[Cross @@ segment[[;; 2]], segment[]] < Pi/2.,
Identity, Reverse]
}
, InfinitePlane[
# + center & /@ ord @ Prepend[segment[[;; 2]], {0, 0, 0}]
]
]
, {segment, Partition[# - center & /@ pts, 3, 1, {1, 1}]}]
]


### Example

n = 5; h = .9;
center = 2 {1, 1, 1};
poly = # + center & /@ ({##, h} & @@@
Reverse @ CirclePoints[Sqrt[1 - h^2], n]);

Graphics3D[{
{Red, sphericalPolygon[poly, center]},
{Blue, sphericalPolygon[Normalize /@ RandomReal[{-1, 1}, {3, 3}]]},
{Green,
sphericalPolygon[
1 + Normalize /@ RandomReal[{-1, 1}, {3, 3}], {1, 1, 1}, .5]},
{Opacity @ .5, Sphere[center, .99], Sphere[{0, 0, 0}, .99],
Sphere[1. + {0, 0, 0}, .49]}
}] Here is a slightly fast method for rendering spherical polygons, based on combining Kuba's idea in this answer with spherical linear interpolation ("slerping") of the polygon's vertices.

The code for the function SphericalPolygon[] is at the end of this post; for now, here are a few examples of its use:

dodFacs = First /@ First[Normal[MapAt[Map[Normalize, #] &,
N[PolyhedronData["Dodecahedron", "GraphicsComplex"]], 1]]];

Graphics3D[MapIndexed[{ColorData[101, #2[]], SphericalPolygon[#1]} &,
dodFacs], Boxed -> False, Lighting -> "Neutral"] BlockRandom[SeedRandom[128, Method -> "MersenneTwister"];
tris = First /@ MeshPrimitives[ConvexHullMesh[RandomPoint[Sphere[], 10]], 2];
Graphics3D[MapIndexed[{ColorData[108, #2[]], SphericalPolygon[#1]} &, tris],
Boxed -> False, Lighting -> "Neutral"]] Now, I present the associated routines. Some auxiliaries first:

(* https://mathematica.stackexchange.com/a/97854 *)
vecang[v1_?VectorQ, v2_?VectorQ] := Module[{n1 = Norm[v1], n2 = Norm[v2]},
2 ArcTan[Norm[v1 n2 + n1 v2], Norm[v1 n2 - n1 v2]]]

(* slightly faster than the equivalent RotationMatrix[{vv1, vv2}];
from http://dx.doi.org/10.1080/10867651.1999.10487509 *)
vectorRotate[vv1_?VectorQ, vv2_?VectorQ] :=
Module[{v1 = Normalize[vv1], v2 = Normalize[vv2], c, d, d1, d2, t1, t2},
d = v1.v2;
If[TrueQ[Chop[1 + d] == 0],
c = UnitVector[3, First[Ordering[Abs[v1], 1]]];
t1 = c - v1; t2 = c - v2; d1 = t1.t1; d2 = t2.t2;
IdentityMatrix - 2 (Outer[Times, t2, t2]/d2 -
2 t2.t1 Outer[Times, t2, t1]/(d2 d1) + Outer[Times, t1, t1]/d1),

c = Cross[v1, v2];
d IdentityMatrix + Outer[Times, c, c]/(1 + d) - LeviCivitaTensor.c]]

(* exponential map for sphere *)
sphereExp[q_?VectorQ, p_?VectorQ] /; Length[q] == Length[p] + 1 :=
With[{n = Norm[p]}, vectorRotate[{0, 0, 1}, q].Append[p Sinc[n], Cos[n]]]

(* inverse of exponential map for sphere *)
sphereLog[q_?VectorQ, p_?VectorQ] /; Length[q] == Length[p] :=
Most[vectorRotate[q, {0, 0, 1}].p]/Sinc[vecang[p, q]]

(* spherical average by Buss and Fillmore, http://dx.doi.org/10.1145/502122.502124 *)
SphericalPolygonCentroid[pts_?MatrixQ] := Module[{k = 0, n = Length[pts], cp, h, pr},
cp = Normalize[Total[pts]];
pr = InternalEffectivePrecision[pts];
While[cp = sphereExp[cp, h = Sum[sphereLog[cp, p], {p, pts}]/n];
k++; Norm[h] > 10^(-pr/2) && k <= 30];
cp]

(* slerp for two points *)
slerp = Compile[{{q1, _Real, 1}, {q2, _Real, 1}, {f, _Real}},
Module[{n1 = Norm[q1], n2 = Norm[q2], omega, so},
omega = 2 ArcTan[Norm[q1 n2 + n1 q2], Norm[q1 n2 - n1 q2]];
If[Chop[so = Sin[omega]] == 0, q1, Sin[{1 - f, f} omega].{q1, q2}/so]]];

(* stripped down version of functions in
https://mathematica.stackexchange.com/a/10385 *)
sphericalLinearInterpolation[pts_?MatrixQ] := Module[{times},
times = Accumulate[Prepend[vecang @@@ Partition[pts, 2, 1, 1], 0]];
{Last[times], sphericalInterpolatingFunction[times, pts]}]

sphericalInterpolatingFunction[times_, vecs_][t_?NumericQ] :=
With[{k = GeometricFunctionsBinarySearch[times, t]},
slerp[vecs[[k]], vecs[[k + 1]], Rescale[t, times[[{k, k + 1}]], {0, 1}]]]


Finally, here is the routine SphericalPolygon[]. Note that I only use adaptive sampling for the spherical polygon's edge via ParametricPlot3D[], and then use "Chebyshev sampling" so that evaluation points are clustered near the edge and center. Also of note is that since the object lives on a sphere, the points can also be used as the normals in VertexNormals.

SphericalPolygon[pts_?MatrixQ, p_: 8] := SphericalPolygon[pts, {p, p}]

SphericalPolygon[pts_?MatrixQ, {p_, q_}] := Module[{ch, cp, en, pt, ql, sp, spl},
cp = SphericalPolygonCentroid[pts];
{en, sp} = sphericalLinearInterpolation[ArrayPad[pts, {{0, 1}}, "Periodic"]];

ch = Sin[π Range[p]/(2 p)]^2; (* rescaled Chebyshev points *)
spl = Most[First[Cases[ParametricPlot3D[sp[t], {t, 0, en}, PlotPoints -> q],
Line[l_] :> l, ∞]]];
ql = Length[spl];

pt = Prepend[Apply[Join, Outer[Normalize[#1.{cp, #2}] &,
Transpose[{Append[Reverse[Most[ch]], 0], ch}],
spl, 1]], cp];

GraphicsComplex[pt, {EdgeForm[],
Polygon[PadLeft[Partition[Range[ql] + 1, 2, 1], {Automatic, 3}, 1]
~Join~
Flatten[Apply[Join[Reverse[#1], #2] &,
Partition[Partition[Range[p ql] + 1, ql],
{2, 2}, {1, 1},
{{1, 1}, {-1, 1}}],
{2}], 1]]}, VertexNormals -> pt]]

• Here is a nice application. Aug 31, 2017 at 13:54