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.

Improve this question
$\endgroup$
4
  • $\begingroup$ But the region is the circle of radius equal to 1, isn't it? $\endgroup$
    – José Antonio Díaz Navas
    Commented Apr 24, 2018 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
    Commented Apr 24, 2018 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
    Commented Apr 24, 2018 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
    Commented Apr 24, 2018 at 18:32

1 Answer 1

Reset to default
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.

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

Your Answer

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

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