6
$\begingroup$

I am trying to find the solid angle taken up by a large set of triangles around a central point. I have normalized the vertices of the triangles to a unit sphere around the central point, But the triangles now are slightly inside of the sphere, because their vertices are on the sphere, but the plane between the points is in the interior of the sphere.

I would like a simple and fast way to project the triangles onto the unit sphere, so I can union the regions, then find the area of the union of all of the regions and compare that to 4*pi, the surface area of the unit sphere. If anyone can help, that would be much appreciated.

EDIT: I start with upwards of 3000 simple triangles in 3 dimensions, with vertices on the unit sphere. I want to find the solid angle taken up by all of the triangles combined. The part of this problem that I find most confusing is that the triangles may overlap, so simply summing their individual areas gives an a solid angle that is much larger than it should be.

$\endgroup$
3

2 Answers 2

14
$\begingroup$

Generate three points on a sphere:

pts = {p1, p2, p3} = Sort[Normalize /@ Table[RandomReal[{-1, 1}], {3}, {3}]];

Then plot a sphere with a region function on the 'positive' side of each of the three planes defined by the origin ({0,0,0}) and successive pairs of points:

oneSide = ContourPlot3D[x^2 + y^2 + z^2 == 1,
  {x, -1, 1}, {y, -1, 1}, {z, -1, 1},
  ContourStyle -> {Opacity[0.3], Green},
  RegionFunction -> 
   Function[{x, y, z}, 
        {x, y, z}.Cross[p1, p2] < 0 && 
        {x, y, z}.Cross[p1, p3] > 0 && 
        {x, y, z}.Cross[p2, p3] < 0],
  Mesh -> None,
  PlotPoints -> 50]

enter image description here

Then use a different color and select the negation of the region defined above:

otherSide = ContourPlot3D[x^2 + y^2 + z^2 == 1,
  {x, -1, 1}, {y, -1, 1}, {z, -1, 1},
  ContourStyle -> {Opacity[0.3], Red},
  RegionFunction -> 
   Function[{x, y, z}, 
   ! ({x, y, z}.Cross[p1, p2] < 0 && 
      {x, y, z}.Cross[p1, p3] > 0 && 
      {x, y, z}.Cross[p2, p3] < 0)],
  Mesh -> None,
  PlotPoints -> 50]

enter image description here

Putting the sides together with the points and the triangle:

Show[oneSide, otherSide, 
 Graphics3D[{PointSize[0.02], Blue, Point /@ pts}], 
 Graphics3D[{Yellow, Polygon[pts]}]]

enter image description here

Be sure to play around with the three-dimensional figure, rotating it as you like, to appreciate the answer.

Note that there is some ambiguity of the definition of your triangle, and you'll have to refine the signs to select the smallest triangle defined by your three points (which is what I presume you seek).

To calculate the region's area

Define an implicit region on the sphere and limited by the three planes:

    triangleArea = ImplicitRegion[
     x^2 + y^2 + z^2 == 1 && 
     {x, y, z}.Cross[p1, p2] < 0 && 
     {x, y, z}.Cross[p1, p3] > 0 && 
     {x, y, z}.Cross[p2, p3] < 0, 
   {x, y, z}]

(*

ImplicitRegion[ x^2 + y^2 + z^2 == 1 && -0.355097 x + 0.523151 y - 0.576394 z < 0 && -0.297567 x + 0.516951 y + 0.801686 z > 0 && 0.885233 x + 0.16822 y - 0.420479 z < 0, {x, y, z}]

*)

Then compute its area:

RegionMeasure[DiscretizeRegion[triangleArea]]

(*

2.217

*)

(This last area happens to be for a different random choice of points than the figure.)

$\endgroup$
1
  • $\begingroup$ Thank you for this answer! I've found that ContourPlot3D + RegionFunction is the best (fast and smooth) method for making Voronoi areas on a curved surfaces. $\endgroup$
    – lesobrod
    Dec 11, 2016 at 15:07
6
$\begingroup$

This is far from ideal and it has limitations with extreme triangles. The aesthetically pleasing images prompt me to post despite this. Perhaps others will improve. I modified David G. Stork function. The solid area is reasonable straight forward. The orange spheres delineate triangle:

ae[a_, b_, c_] := 
 Developer`PartitionMap[VectorAngle @@ # &, {a, b, c}, 2, 1, 1]
sa[w1_, w2_, w3_] := 
  With[{ts = Plus @@ {w1, w2, w3}/2}, 
   4 ArcTan[
     Sqrt[Times @@ (Tan /@ (0.5 (ts - #) & /@ {0, w1, w2, w3}))]]];
solid[a_, b_, c_] := sa @@ ae[a, b, c]
trg[a_, b_, c_] := 
 Show[Graphics3D[{EdgeForm[Black], Orange, 
    Sphere[#, 0.1] & /@ {a, b, c}, 
    Text[Style[solid[a, b, c], White, Bold, 
      FontFamily -> "Kartika"], {0.5, -0.1, 0}]}, Background -> Black,
    Boxed -> False], 
  ParametricPlot3D[{Sin[u] Cos[v], Sin[u] Sin[v], Cos[u]}, {u, 0, 
    Pi}, {v, 0, 2 Pi}, 
   MeshFunctions -> 
    Function[{x, y, z}, 
     Sign[{x, y, z}.Cross[a, b] {x, y, z}.Cross[c, a] {x, y, z}.Cross[
         b, c]]], Mesh -> {{0}}, MeshShading -> {Red, Blue}, 
   MeshStyle -> {Yellow, Thick}, PlotPoints -> 100], 
  ViewPoint -> 10 (a + b + c)/3]

An example:

sm[n_] := 
 Table[trg @@ (Normalize /@ RandomReal[{-1, 1}, {3, 3}]), {n}]
gr = Grid[Partition[sm[16], 4]];

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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