Checking that a point is in a spherical polygon

linear algebraspherical-geometry

I am trying to find an algorithm to check if a point $p$ on the surface of a sphere lies within a spherical polygon with vertices at $v_1, v_2, \dots, v_J$. Any polygons I am considering will be "small", so the smaller of the two potential polygons on the surface of the sphere. I also expect my polygons to be spherical squares and hexagons, that is nothing too irregular. My points are arranged so that any polygon will be CCW.

Right now, I am thinking of using the following. I have not yet found a counter-example to suggest it shouldn't work, but given how hard it has been to find such an algorithm I am wondering if I am missing something as this seems very straighforward.

If $\hat{\mathbf{n}}_{j,p}$ is unit vector normal to the great circle connecting $v_j$ and $p$, then the interior angle for the points $v_j$–$p$–$v_{j+1}$ will be $$\alpha_j = \arccos(\hat{\mathbf{n}}_{j,p}\cdot\hat{\mathbf{n}}_{p,j+1})$$
and the point is in the polygon iff
$$\sum_{j=1}^J\alpha_j=2\pi.$$

Can anyone suggest a reason why this would not work? If one checked that angles were CW or CCW, is there any reason it could not be extended for more complex polygons with more of a "Pacman" shape? (Apologies for poor vernacular!)

Best Answer

Using the suggestions in Will's comments, if a $J$-sided polygon is convex and CCW, then $\hat{\mathbf{n}}_j$ points into the polygon and so for a point $\mathbf{p}$ if $\hat{\mathbf{n}}_j\cdot\mathbf{p}>0\,\forall j\in J$ then the point is in the polygon. Should the $\hat{\mathbf{n}}_j\cdot\mathbf{p}=0$ for any side $j$ then the point lies on that line.

Related Solutions

Converting between two spherical coordinate systems with an application to astronothe

Let the standard astronomical system be

$$ X(\theta,\phi, R) = R \left( \cos \theta \cos \phi, \ \cos \theta \ \sin \phi, \sin \theta\right)$$

with $\theta=0$ equatorial plane, $\theta=\pm \pi/2$ the polar directions and $\phi=0$ the Greenwich meridian.

In order to find the map to a second system with polar direction $\{\theta_0,\phi_0\}$ you have to construct the rotation matrix in this plane spanned by the two polar vectors with an angle of rotation parameter $\theta_0$.

So the single task is to find the generator of such a rotation. The usual method is to take the cross product and norm it to unity

$$ \vec {n'} = (0,0,1) \times X(\theta_0,\phi_0,1)$$ But its much easier, the rotation plane is the plane of constant $\phi_0$ and the unit vector of rotation is simply the unit vector in the equatorial plane at $\phi=\phi_0\pm\pi/2$

Then the group of rotations in this plane is

$$\text{MatrixExp}\left[t \left(\left( \begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \\ \end{array} \right) \sin \left(\phi _0\right)+\left( \begin{array}{ccc} 0 & 0 & 1 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \\ \end{array} \right) \cos \left(\phi _0\right)\right)\right] = \left( \begin{array}{ccc} \cos ^2\left(\frac{t}{2}\right)-\sin ^2\left(\frac{t}{2}\right) \cos \left(2 \phi _0\right) & \sin ^2\left(\frac{t}{2}\right) \sin \left(2 \phi _0\right) & \sin (t) \cos \left(\phi _0\right) \\ \sin ^2\left(\frac{t}{2}\right) \sin \left(2 \phi _0\right) & \sin ^2\left(\frac{t}{2}\right) \cos \left(2 \phi _0\right)+\cos ^2\left(\frac{t}{2}\right) & -\sin (t) \sin \left(\phi _0\right) \\ \sin (t) \left(-\cos \left(\phi _0\right)\right) & \sin (t) \sin \left(\phi _0\right) & \cos (t) \\ \end{array} \right)$$

Setting $t=\pi/2 - \theta_0$ for t=0, if th polar axes coincide, one finally has the compact formula

$$X' = \left( \begin{array}{c} \cos \left(\phi _0\right) \left(\sin \left(\theta _0\right) \cos (\theta ) \cos \left(\phi +\phi _0\right)+\sin (\theta ) \cos \left(\theta _0\right)\right)+\cos (\theta ) \sin \left(\phi _0\right) \sin \left(\phi +\phi _0\right) \\ \cos (\theta ) \sin \left(\phi +\phi _0\right) \cos \left(\phi _0\right)-\sin \left(\phi _0\right) \left(\sin \left(\theta _0\right) \cos (\theta ) \cos \left(\phi +\phi _0\right)+\sin (\theta ) \cos \left(\theta _0\right)\right) \\ \sin (\theta ) \sin \left(\theta _0\right)-\cos (\theta ) \cos \left(\theta _0\right) \cos \left(\phi +\phi _0\right) \\ \end{array} \right)$$

Disclaimer: This is a result produced by Mathematica GPT, so better control all of the results. For Mathematica users

  Manipulate[ Graphics3D[  col = {Red, Green, Blue}; 
   {Array[({col[[#]],   Arrow[{{0, 0, 0}, UnitVector[3, #]}]} &), 3],
    Array[({col[[#]], Arrow[{{0, 0, 0}, M[\[Theta], \[Phi]][[#]]}]} &), 3]},
     PlotRange -> {{-1.2, 1.2}, {-1.2, 1.2}, {-1.2, 1.2}}], 
    {{\[Theta], \[Pi]/2},  0, \[Pi]}, 
     {{\[Phi], 0}, -\[Pi], \[Pi]} ]

Related Question