# Finding the solid angle subtended by a closed surface of data points in 3D space

Suppose, I have a closed 2D region of data points.

 𝓡 =
DiscretizeRegion[
ImplicitRegion[(x - 2)^2 + (y - 3)^2 <= 1 && z == 5, {x, y, z}],
MaxCellMeasure -> 0.0001];
points = RandomPoint[𝓡, 500];


Now, I want to compute the solid angle subtended by the region so covered at the origin. The formula is: $$\int\int_{\mathcal{R}}\sin \theta \;d\theta\;d \phi$$

I tried three ways:

pointsInPolar = ToPolarCoordinates[#] & /@ points;
Integrate[Sin[θ],
Element[{r, θ, ϕ},
MeshRegion[pointsInPolar, Point[Range[500]]]]]
(*472.178*)
Integrate[Sqrt[x^2 + y^2]/
Sqrt[x^2 + y^2 + z^2], {x, y, z} ∈ MeshRegion[points, Point[Range[500]]]]
(*291.617*)
Integrate[Sqrt[x^2 + y^2]/
Sqrt[x^2 + y^2 + z^2], {x, y, z} ∈ ConvexHullMesh[points]]


The last one flags an error.

How to go about doing it?

Edit 1: Here is a little better region, in the sense that I can easily determine its solid angle:

𝓡 =
DiscretizeRegion[
ImplicitRegion[x^2 + y^2 <= 1/2 && z == 1/Sqrt[2], {x, y, z}],
MaxCellMeasure -> 0.0001];
points = RandomPoint[𝓡, 500000];


The exact answer is $2 \pi (1-\cos \theta)$ which in this case evaluates to 1.8403.

• But the region is the circle of radius equal to 1, isn't it? Apr 24, 2018 at 17:50
• @JoséAntonioDíazNavas No, its a disc in hindsight. I constructed this simple example to understand how to go about finding a solid angle for more arbitrary 2D regions. Apr 24, 2018 at 17:59
• yes I know the circle is an ellipse due to the perspective. You must obtain the angle between the normal to the region and the vector pointing from the origin to the geometrical center in order to project the region to the plane perpendicular to this one. Or am I missing something ? Apr 24, 2018 at 18:19
• @JoséAntonioDíazNavas I didn't mean ellipse. By disc I meant circle together with its interior points. And yes the projection idea is correct. The problem seems to be to convert the random data set to a properRegion object in order to use Integrate. Apr 24, 2018 at 18:32

It seems you want to compute the convex hull in spherical coordinates and integrate over that:

ph = ConvexHullMesh[(ToSphericalCoordinates[#] & /@ points)[[All,2 ;; 3]]]
Integrate[ Sin[t], {t, p} \[Element] ph]


0.0630096

also

range = MinMax /@
Transpose[(ToSphericalCoordinates[#] & /@ points)[[All, 2 ;; 3]]];
NIntegrate[ Sin[t] Boole[RegionMember[ph, {t, p}]],
Evaluate@{t, Sequence @@ range[[1]]},
Evaluate@{p, Sequence @@ range[[2]]}]


0.0629917

edit: a bit of warning, this breaks if the region breaks the phi=+/-Pi boundary:

𝓡 = DiscretizeRegion[
ImplicitRegion[(y)^2 + (z)^2 < 1/2 && x == -1/Sqrt[2], {x, y, z}],
MaxCellMeasure -> 0.0001];
points = RandomPoint[𝓡, 500];
spherical = (ToSphericalCoordinates[#] & /@ points)[[All, 2 ;; 3]];
Show[{ph=ConvexHullMesh[spherical], Graphics[Point[spherical]]},
Axes -> True, AspectRatio -> 1]


Integrate[ Sin[t], {t, p} [Element] ph]

8.73094 (*wrong*)

This is readily fixed by adding 2 Pi to the negative phi values before generating the convex hull, but you need to devise some heuristic to detect that condition.

• I might be wrong. Could you also do the same in cartesian coordinates? Apr 25, 2018 at 1:28
• try working an example where you can more easily work out the exact result. Apr 25, 2018 at 1:48
• I updated my question with a better example. Your code does produce a close enough result. Apr 25, 2018 at 3:16
• good, it actually converges nice if you up the # points quite a lot. Note there is an issue if the area includes the jump at phi=+/-pi. Apr 25, 2018 at 14:47