Skip to main content
Mathematica

Return to Answer

added 85 characters in body
Source Link
J. M.'s missing motivation ♦
J. M.'s missing motivation ♦
  • 125.5k
  • 11
  • 407
  • 581

Quoting Wikipedia (https://en.wikipedia.org/wiki/Solid_angleWikipedia), "Solid angle is the two-dimensional angle in three-dimensional space that an object subtends at a point."

Now think about composite objects. In 2D once again, suppose we have a point and a polygon whose boundary edges are a number of straight line segments. We sum the signed angles subtended at the point by all edges. If the point is within the polygon, the arcs entirely encircle the point and the angles sum to 2 pi$2\pi$ radians. If the point is outside, the angles sum to 0 radians. This is because edges on the near side of the polygon have positive angles, while edges on the far side have negative angles. They entirely cancel. These results are true even if the polygon is non-convex -- appearing to double back on itself viewed from the point.

And in 3D, suppose we have a point and a polyhedron with polygonal faces. We sum the signed solid angles subtended at the point by all faces. If the point is within the polyhedron, the entire unit-sphere is covered and the solid angles sum to 4 pi$4\pi$ (the area of a unit-radius sphere) steradians. If the point is outside, the sum is 0 steradians because solid angles of polygons facing away from the point are positive and those facing toward are negative. As in 2D, this is true even if the polyhedron is non-convex.

So the problem boils down to summing the signed solid angle of the cow's polygonal faces subtended at the point. The total is 4 pi$4\pi$ if the point is inside or 0 if outside.

MathematicaMathematica represents the cow using only triangular faces. Simple formulas for the signed solid angle subtended by a planar triangle appear in the following two papers:

Van Oosterom and Strackee, "The Solid Angle of a Plane Triangle""The Solid Angle of a Plane Triangle", IEEE Transactions on Biomedical Engineering 30(2) 125-126.

Eriksson, "On the Measure of Solid Angles""On the Measure of Solid Angles", Mathematics Magazine 63(3) 184-187.

I assume (but have not checked) that MathematicaMathematica winds faces consistently, meaning a normal vector computed from the three vertices points outward from the polyhedron. The Det[]Det[] function used below relies on this. If a face is wound incorrectly, its solid angle will have the wrong sign.

There can be numerical problems. We sum a great number (5804 faces for the cow) of small quantities. Results might not be exactly 4 pi$4\pi$ or exactly 0.

Execution time might be another problem. On a 2008 MacBook Pro running MathematicaMathematica 8.0.4.0, evaluating one point takes about 0.35 seconds. So evaluating the OP's 1000 points would take about 6 minutes.

Quoting Wikipedia (https://en.wikipedia.org/wiki/Solid_angle), "Solid angle is the two-dimensional angle in three-dimensional space that an object subtends at a point."

Now think about composite objects. In 2D once again, suppose we have a point and a polygon whose boundary edges are a number of straight line segments. We sum the signed angles subtended at the point by all edges. If the point is within the polygon, the arcs entirely encircle the point and the angles sum to 2 pi radians. If the point is outside, the angles sum to 0 radians. This is because edges on the near side of the polygon have positive angles, while edges on the far side have negative angles. They entirely cancel. These results are true even if the polygon is non-convex -- appearing to double back on itself viewed from the point.

And in 3D, suppose we have a point and a polyhedron with polygonal faces. We sum the signed solid angles subtended at the point by all faces. If the point is within the polyhedron, the entire unit-sphere is covered and the solid angles sum to 4 pi (the area of a unit-radius sphere) steradians. If the point is outside, the sum is 0 steradians because solid angles of polygons facing away from the point are positive and those facing toward are negative. As in 2D, this is true even if the polyhedron is non-convex.

So the problem boils down to summing the signed solid angle of the cow's polygonal faces subtended at the point. The total is 4 pi if the point is inside or 0 if outside.

Mathematica represents the cow using only triangular faces. Simple formulas for the signed solid angle subtended by a planar triangle appear in the following two papers:

Van Oosterom and Strackee, "The Solid Angle of a Plane Triangle", IEEE Transactions on Biomedical Engineering 30(2) 125-126.

Eriksson, "On the Measure of Solid Angles", Mathematics Magazine 63(3) 184-187.

I assume (but have not checked) that Mathematica winds faces consistently, meaning a normal vector computed from the three vertices points outward from the polyhedron. The Det[] function used below relies on this. If a face is wound incorrectly, its solid angle will have the wrong sign.

There can be numerical problems. We sum a great number (5804 faces for the cow) of small quantities. Results might not be exactly 4 pi or exactly 0.

Execution time might be another problem. On a 2008 MacBook Pro running Mathematica 8.0.4.0, evaluating one point takes about 0.35 seconds. So evaluating the OP's 1000 points would take about 6 minutes.

Quoting Wikipedia, "Solid angle is the two-dimensional angle in three-dimensional space that an object subtends at a point."

Now think about composite objects. In 2D once again, suppose we have a point and a polygon whose boundary edges are a number of straight line segments. We sum the signed angles subtended at the point by all edges. If the point is within the polygon, the arcs entirely encircle the point and the angles sum to $2\pi$ radians. If the point is outside, the angles sum to 0 radians. This is because edges on the near side of the polygon have positive angles, while edges on the far side have negative angles. They entirely cancel. These results are true even if the polygon is non-convex -- appearing to double back on itself viewed from the point.

And in 3D, suppose we have a point and a polyhedron with polygonal faces. We sum the signed solid angles subtended at the point by all faces. If the point is within the polyhedron, the entire unit-sphere is covered and the solid angles sum to $4\pi$ (the area of a unit-radius sphere) steradians. If the point is outside, the sum is 0 steradians because solid angles of polygons facing away from the point are positive and those facing toward are negative. As in 2D, this is true even if the polyhedron is non-convex.

So the problem boils down to summing the signed solid angle of the cow's polygonal faces subtended at the point. The total is $4\pi$ if the point is inside or 0 if outside.

Mathematica represents the cow using only triangular faces. Simple formulas for the signed solid angle subtended by a planar triangle appear in the following two papers:

Van Oosterom and Strackee, "The Solid Angle of a Plane Triangle", IEEE Transactions on Biomedical Engineering 30(2) 125-126.

Eriksson, "On the Measure of Solid Angles", Mathematics Magazine 63(3) 184-187.

I assume (but have not checked) that Mathematica winds faces consistently, meaning a normal vector computed from the three vertices points outward from the polyhedron. The Det[] function used below relies on this. If a face is wound incorrectly, its solid angle will have the wrong sign.

There can be numerical problems. We sum a great number (5804 faces for the cow) of small quantities. Results might not be exactly $4\pi$ or exactly 0.

Execution time might be another problem. On a 2008 MacBook Pro running Mathematica 8.0.4.0, evaluating one point takes about 0.35 seconds. So evaluating the OP's 1000 points would take about 6 minutes.

deleted 7 characters in body
Source Link
kglr
kglr
  • 400.5k
  • 18
  • 488
  • 929
threshold = 0.05; (* empirical cutoff *)
pointInsidePolygonList = Select[Take[lPoints, 100], polyhedronSolidAngle[#, polygonCoordList] > threshold &];

(I regret not knowing how to include an image with this answer.)

 
Show[g, ListPointPlot3D[pointInsidePolygonList]]

Mathematica graphics

threshold = 0.05; (* empirical cutoff *)
pointInsidePolygonList = Select[Take[lPoints, 100], polyhedronSolidAngle[#, polygonCoordList] > threshold &];

(I regret not knowing how to include an image with this answer.)

Show[g, ListPointPlot3D[pointInsidePolygonList]]

threshold = 0.05; (* empirical cutoff *)
pointInsidePolygonList = Select[Take[lPoints, 100], polyhedronSolidAngle[#, polygonCoordList] > threshold &];
 
Show[g, ListPointPlot3D[pointInsidePolygonList]]

Mathematica graphics

Source Link
Bob Werner
Bob Werner
  • 421
  • 1
  • 3
  • 8

One way is to compute the solid angle subtended by the cow viewed at the point by summing signed solid angles corresponding to the cow's polygonal faces. If the total is 4 pi, the point is inside the cow; if 0, outside.

Background

Quoting Wikipedia (https://en.wikipedia.org/wiki/Solid_angle), "Solid angle is the two-dimensional angle in three-dimensional space that an object subtends at a point."

As a 2D analogy, think of a point and a line segment in a 2D plane. Center a unit-radius circle on the point. Project the line segment onto the circle by connecting the endpoints to the center point. The arc projected on the unit-radius circle has a certain length, the angle subtended by the line segment. It is measured in radians.

Now move up to 3D. Think of a point and a polygon in 3D. Center a unit-radius sphere on the point, and use a central projection to map the polygon onto the sphere. The area of the projected spherical polygon on the unit-radius sphere is the solid angle subtended by the polygon. The measure is steradians.

These angles are signed. In 2D, if the circular arc evolves counter-clockwise as the line is traversed from beginning to end, the angle is positive; if clockwise, negative. In 3D, if the polygon's normal vector is directed away from the point, the solid angle is positive; if toward, negative. (The normal vector is orthogonal to the polygon plane and by convention points away from the polyhedron/cow.)

Now think about composite objects. In 2D once again, suppose we have a point and a polygon whose boundary edges are a number of straight line segments. We sum the signed angles subtended at the point by all edges. If the point is within the polygon, the arcs entirely encircle the point and the angles sum to 2 pi radians. If the point is outside, the angles sum to 0 radians. This is because edges on the near side of the polygon have positive angles, while edges on the far side have negative angles. They entirely cancel. These results are true even if the polygon is non-convex -- appearing to double back on itself viewed from the point.

And in 3D, suppose we have a point and a polyhedron with polygonal faces. We sum the signed solid angles subtended at the point by all faces. If the point is within the polyhedron, the entire unit-sphere is covered and the solid angles sum to 4 pi (the area of a unit-radius sphere) steradians. If the point is outside, the sum is 0 steradians because solid angles of polygons facing away from the point are positive and those facing toward are negative. As in 2D, this is true even if the polyhedron is non-convex.

One reference for the 3D case is Spiegel, "Schaum's Outline Series: Theory and Problems of Vector Analysis", McGraw-Hill, Ch 6 Sec 27.

Code

So the problem boils down to summing the signed solid angle of the cow's polygonal faces subtended at the point. The total is 4 pi if the point is inside or 0 if outside.

Mathematica represents the cow using only triangular faces. Simple formulas for the signed solid angle subtended by a planar triangle appear in the following two papers:

Van Oosterom and Strackee, "The Solid Angle of a Plane Triangle", IEEE Transactions on Biomedical Engineering 30(2) 125-126.

Eriksson, "On the Measure of Solid Angles", Mathematics Magazine 63(3) 184-187.

This is adapted from the OP's question:

lX = RandomVariate[UniformDistribution[{-0.4, 0.4}], {1000}];
lY = RandomVariate[UniformDistribution[{-0.2, 0.2}], {1000}];
lZ = RandomVariate[UniformDistribution[{-0.2, 0.2}], {1000}];
lPoints = Thread[{lX, lY, lZ}];

cow = ExampleData[{"Geometry3D", "Cow"}];

polygonCoordList = Cases[Normal[cow], Polygon[x_, ___] :> x, {0, Infinity}];

I assume (but have not checked) that Mathematica winds faces consistently, meaning a normal vector computed from the three vertices points outward from the polyhedron. The Det[] function used below relies on this. If a face is wound incorrectly, its solid angle will have the wrong sign.

The Eriksson paper's formula for the solid angle of a single triangular face:

Clear[triangleSolidAngle];

triangleSolidAngle[pnt_, {vertex1_, vertex2_, vertex3_}] := Module[
{dir1, dir2, dir3},

(* unit-length direction vectors from pnt to the three vertices *)

dir1 = Normalize[vertex1 - pnt];
dir2 = Normalize[vertex2 - pnt];
dir3 = Normalize[vertex3 - pnt];

(* solid angle *)

2*ArcTan[1 + dir1 . dir2 + dir2 . dir3 + dir3 . dir1, Det[{dir1, dir2, dir3}]]
]

The solid angle of an entire triangular-faced polyhedron:

Clear[polyhedronSolidAngle];

polyhedronSolidAngle[pnt_, triangleList_] := Total[triangleSolidAngle[pnt, #] & /@ triangleList]

There can be numerical problems. We sum a great number (5804 faces for the cow) of small quantities. Results might not be exactly 4 pi or exactly 0.

Execution time might be another problem. On a 2008 MacBook Pro running Mathematica 8.0.4.0, evaluating one point takes about 0.35 seconds. So evaluating the OP's 1000 points would take about 6 minutes.

Trial run of only the first 100 points:

threshold = 0.05; (* empirical cutoff *)
pointInsidePolygonList = Select[Take[lPoints, 100], polyhedronSolidAngle[#, polygonCoordList] > threshold &];

(I regret not knowing how to include an image with this answer.)

Show[g, ListPointPlot3D[pointInsidePolygonList]]