8
$\begingroup$

I need to visualize the points of a spherical mesh on a squared plane (this question is different than this one, I believe). For instance, let us first generate points pts on a unit sphere. (Original question used a function from MathWorld to generate these points. The updated version has the more efficient way RandomPoint[Sphere[], num]).

Then, let us use the points generated to create a mesh (Delaunay Mesh in this case):

num = 1000;(*Desired number of points*)
pts = RandomPoint[Sphere[], num];
Dmesh = DelaunayMesh[pts];
Show[Dmesh, Graphics3D[{Red, PointSize[0.01], Point[pts]}]]

Gives:

enter image description here

Is it possible to visualize a squared 'planar' version of the mesh while minimizing distortion? Any approximation will suffice since I perform computations on the sphere (the square is only useful for visualization).

I've tried creating a planar graph, but I realized that many meshes produced in this way are not planar. Then I tried making a graph with a SpringElectricalEmbedding layout (and others), but I cannot force these layouts to be 'squares':

(*Get vertex neighbors from the mesh*)
neighbors = MeshCells[Dmesh, 1][[All, 1]];
neighbors = 
  Table[neighbors[[a]][[1]] <-> neighbors[[a]][[2]], {a, 1, 
    Length[neighbors]}];

(*Create a graph layout that minimizes 'energy'*)
Graph[neighbors, GraphLayout -> "SpringElectricalEmbedding"]

enter image description here

The ideal solution would be something like:

PlanarSphereMesh[Dmesh]

enter image description here

Any hints towards this goal will be appreciated. To be clear, I don't require a 'planar graph' to be created, that was only an idea I had at first. Any method of visualization would suffice, as long as it preserves the vertex connectivity. Even more, if the 'boundaries' of the plane representation requires it, we could remove the edges at the boundaries, and only visualize the connections for 'inner' vertexes.

$\endgroup$
9
  • 4
    $\begingroup$ Just a side remark: You can generate the random points also with RandomPoint[Sphere[],num]. $\endgroup$ Mar 9 '20 at 21:42
  • 1
    $\begingroup$ "but I realized that any mesh on a sphere cannot be a planar graph" Not true, any graph that can be embedded on the sphere with no edge crossings is a planar graph. $\endgroup$
    – Szabolcs
    Mar 9 '20 at 22:00
  • 1
    $\begingroup$ "The ideal solution would be something like:" You can't do that: there must be much higher distortion. GraphLayout -> TutteEmbedding will work in principle, but it won't produce a usable result: the distortion will be extreme. Instead: $\endgroup$
    – Szabolcs
    Mar 9 '20 at 22:00
  • 4
    $\begingroup$ Borrow a simple solution from map making, for example stereographic projection. Map making is all about how to project a sphere to the plane. $\endgroup$
    – Szabolcs
    Mar 9 '20 at 22:01
  • 1
    $\begingroup$ @TumbiSapichu The boundary is the projection boundary, white vs. light blue in my answer. $\endgroup$
    – kirma
    Mar 10 '20 at 16:40
13
$\begingroup$

I wouldn't try to do this on a square (assuming you want at least some sort of systematic correspondence between edges on the sphere and their projections), at least if you don't want crossings on square sides.

As @Szabolcs stated, you can't really do this without extreme distortion. You may use map projections, though, picking your poison (for the mesh boundary):

SeedRandom[1];
With[{mapping = 
   CoordinateTransformData["Cartesian" -> "Spherical", "Mapping"]},
 GeoGraphics[
  MeshPrimitives[RegionBoundary@DelaunayMesh[RandomPoint[Sphere[], 100]], 1] /.
   {{x_Real, y_, z_} :> 
     mapping[{x, y, z}][[2 ;; 3]] / Degree - {90, 0}, Line -> GeoPath},
  GeoProjection -> "Equirectangular", GeoBackground -> LightBlue,
  GeoRange -> "World"]]

enter image description here

Here lines of the boundary of the Delaunay mesh are converted to geodesic paths between latitude-longitude points and are projected to an "empty" map, in this case choosing distortions of an equirectangular projection. Note that many graph edges cross the edge of the projection.

Instead of equirectangular projection, one can use, for instance, Hammer ("Hammer"):

enter image description here

The azimuthal equidistant ("AzimuthalEquidistant") projection can avoid crossings of the projection boundary (although lines close to it tend to have rendering problems), but with extreme distortion:

enter image description here

If you absolutely want to use PlanarGraph, you can do it on the boundary region, but it'll be laid out inside a triangle (because your boundary mesh consists of only triangles) with much worse distortion than the map projection:

SeedRandom[1];
Graph[
 MeshPrimitives[
   RegionBoundary@DelaunayMesh[RandomPoint[Sphere[], 100]], 1] /.
  Line[{x__}] -> UndirectedEdge[x], GraphLayout -> "TutteEmbedding"]

enter image description here

You might be able to apply some sort of function on vertex coordinates to spread them more evenly over map projections while maintaining graph properties, but that's too convoluted for this answer.

$\endgroup$
2
  • $\begingroup$ I will explore the Hammer projection, looks promising. As for the planar graph, maybe enforcing the position of the nodes (i.e. with VertexCoordinates) would help? Or is it the triangulation a fundamental constraint on the shape of the 'boundary' (triangle, polygon, etc)? $\endgroup$ Mar 10 '20 at 17:27
  • $\begingroup$ @TumbiSapichu Well, the Delaunay mesh boundary consists by definition of triangles, so although a graph on a genus-0 surface is planar, the "outermost" boundary of such a planar layout will be some of those triangles. Surely it's possible to adjust layout of vertices (or make edges non-straight), but the success is limited. $\endgroup$
    – kirma
    Mar 11 '20 at 5:41
7
$\begingroup$
SeedRandom[1];
dm = RegionBoundary @ DelaunayMesh[RandomPoint[Sphere[], 100]];
mc = MeshCoordinates[dm];

We can use a combination of GeoPositionXYZ and GeoPosition to get 2D projections of mc:

gp = Most /@ First @ GeoPosition @ GeoPositionXYZ[mc, 1.];
indices = List /@ MeshCells[dm, 1][[All, 1]];

Row[GeoGraphics[GeoPath[Extract[gp, indices]], GeoRange -> "World", 
      GeoProjection -> #, GeoBackground -> None, ImageSize -> 400, PlotLabel -> #] & /@ 
       {"Equirectangular", "LambertAzimuthal"}, 
  Spacer[10]]

enter image description here

We can post-process to get straight lines:

GeoGraphics[GeoPath[Extract[gp, indices]], GeoRange -> "World", 
   GeoProjection -> "Equirectangular", GeoBackground -> None, 
   ImageSize -> Large][[1]] /. 
 Line[x_] :> {Line /@ x[[All, {1, -1}]], 
    Red, PointSize[Medium], Point /@ x[[All, 1]]}

enter image description here

$\endgroup$
1
  • $\begingroup$ Ah, GeoPositionXYZ[..., 1.]! I was trying to do this stuff with it but didn't know that the last argument would help. $\endgroup$
    – kirma
    Mar 10 '20 at 10:32

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.