0
$\begingroup$

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];

enter image description here

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.

$\endgroup$
  • $\begingroup$ But the region is the circle of radius equal to 1, isn't it? $\endgroup$ – José Antonio Díaz Navas Apr 24 '18 at 17:50
  • $\begingroup$ @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. $\endgroup$ – Subho Apr 24 '18 at 17:59
  • $\begingroup$ 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 ? $\endgroup$ – José Antonio Díaz Navas Apr 24 '18 at 18:19
  • $\begingroup$ @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. $\endgroup$ – Subho Apr 24 '18 at 18:32
2
$\begingroup$

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]

enter image description here

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.

$\endgroup$
  • $\begingroup$ I might be wrong. Could you also do the same in cartesian coordinates? $\endgroup$ – Subho Apr 25 '18 at 1:28
  • $\begingroup$ try working an example where you can more easily work out the exact result. $\endgroup$ – george2079 Apr 25 '18 at 1:48
  • $\begingroup$ I updated my question with a better example. Your code does produce a close enough result. $\endgroup$ – Subho Apr 25 '18 at 3:16
  • $\begingroup$ 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. $\endgroup$ – george2079 Apr 25 '18 at 14:47

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.