Find three closest points which triangle contains a given point on a sphere

Question

I have a 3D sphere with points on its surface. These points are represented in spherical coordinates, so azimuth, elevation and r.

For example my dataset is a matrix with all the available points on the given sphere:

  azimuth       elevation     r
[[  0.          90.           1.47      ]
 [  0.          85.2073787    1.47      ]
 [  0.          78.16966379   1.47      ]
 [  0.          70.30452954   1.47      ]
 [  0.          62.0367492    1.47      ]
 [  0.          53.56289304   1.47      ]
 [  0.          45.           1.47      ]
 [  0.          36.43710696   1.47      ]
 [  0.          27.9632508    1.47      ]
 [  0.          19.69547046   1.47      ]
 [  0.          11.83033621   1.47      ]
 [  0.           4.7926213    1.47      ]
 [  0.           0.           1.47      ]
 [  0.          -4.7926213    1.47      ]
 [  0.         -11.83033621   1.47      ]
 [  0.         -19.69547046   1.47      ]
 [  0.         -27.9632508    1.47      ]
 [  0.         -36.43710696   1.47      ]
 [  0.         -45.           1.47      ]
 [  0.         -53.56289304   1.47      ]
 [  0.         -62.0367492    1.47      ]
 [  0.         -70.30452954   1.47      ]
 [  0.         -78.16966379   1.47      ]
 [  0.         -85.2073787    1.47      ]
 [  0.         -90.           1.47      ]
 [  1.64008341  -1.6394119    1.47      ]
 [  1.64008341   1.6394119    1.47      ]
 [  2.37160039   8.01881788   1.47      ]
 [  2.37160039  -8.01881788   1.47      ]
 [  2.80356493 -15.58649429   1.47      ]
 [  2.80356493  15.58649429   1.47      ]
 [  3.16999007  23.70233802   1.47      ]
 [  3.16999007 -23.70233802   1.47      ]
 [  3.56208248 -32.09871039   1.47      ]
 [  3.56208248  32.09871039   1.47      ]
 [  4.04606896  40.63141594   1.47      ]
 [  4.04606896 -40.63141594   1.47      ]
 [  4.1063771   -4.09587122   1.47      ]
ecc...

NB: I am omitting the full data matrix on purpose since it contains quite a lot of data. If needed/required to make the problem fully reproducible, I will provide the full data.

This matrix represents something like this image:

Given an arbitrary point, I would like to compute the 3 closest points in the dataset that "contain" the input point.

My code so far is the following:

def compute_three_closest_positions(self, azimuth_angle, elevation):
    requested_position = np.array([azimuth_angle, elevation, 0])

    # computing the absolute difference between the requested angles and the available one in the dataset
    result = abs(self.sourcePositions - requested_position) #subtracting between dataset and requested point
    result = np.delete(result, 2, 1) # removing the r data 
    result = result.sum(axis=1) #summing the overall difference for each point

    # returning index of the closest points
    indexes = result.argsort()[:3] 

    closest_points = self.sourcePositions[indexes]
    
    return closest_points 

Basically I am subtracting the requested azimuth and elevation from all the points in the matrix dataset (self.sourcePositions), then i sum these differences for each point, compute the top 3 minimum indexes and then I use these indexes to access the points in my dataset.

The code works fine, the problem is that sometimes I get 3 closest points that do not contain the requested point.

Examples:

Wrong one:

  Requested point: azimut, elevation, distance
    [200   0   1.47]
    # As you might notice, the requested point is not inside the triangle created by the 3 closest points
    Three closest points: azimut, elevation, distance
    [[199.69547046   0.           1.47      ]
     [199.40203659   5.61508214   1.47      ]
     [199.40203659  -5.61508214   1.47      ]]

Good one:

Requested position:
[190   0   1.47]
# As you can notice, in this case the requested point is inside the triangle generated by the closest 3 points
Three closest points:
[[191.83033621   0.           1.47      ]
 [188.02560265   2.34839855   1.47      ]
 [188.02560265  -2.34839855   1.47      ]]

How could I fix this problem? I would like to get the 3 closest points which "triangle" (I am on a spherical surface, so it is not a real triangle) contain my requested point.

Answer 1

From starters azimut+elevation does not sound right for this I would prefer latitude+longitude therms because the ones you used are something else.

Now as mentioned in the comments if your points form a regular grid you can make an equation that maps topology into point and back where topology is described by integer index of your point in array (or 2 indexes). However if you are not capable of infering this or the grid is not regular then you can do this:

1. reorder your array so its sorted by latitude and long

   so longitude is in <0,2*PI> and latitude in <-PI/2,+PI/2> so I would sort your array by both while latitude has bigger priority (weight). Lets call this new array pnt[]

2. mapping point p to closest vertex of sphere

   simply binary search pnt[] until found point index ix that has bigest smaler or equal latitude than p

   then linearly search from ix (binary search could be used if you reorder your pnt[] to slices or remember how many points per latitude there is) until found biggest latitude that is smaller or equal to p.

   Now pnt[ix] is the "closest" point on sphere to p

3. list closest neighbors to pnt[ix]

   so simply pnt[ix] is from left side of longitude so pnt[ix+1] should be the next point (in case you cross the array size you or pole which needs to check the points with brute force but just the last few in array)

   now we just need to find corresponding points below or above these 2 points (depends on which side your p is). So find 2 closest points to p in the same way as #2 but for smaller and bigger latitude (one slice before,after). this will get you 3*2 points which forms (using always the 2 points found first) 4 potential triangles.

4. test possible triangles

   so you have potential triangle p0,p1,p2 that is "parallel" to sphere surface. So simply project your point onto its plane:

   u = p1-p0
   v = p2-p0
   v = cross(u,v)
   v = cross(u,v)
   p' = p0 + dot(p-p0,u)*u + dot(p-p0,v)*v
   

   so u,v are basis vectors and p0 is the origin of a plane... Now just test of p' is inside triangle so either use 2D and barycentric coordinates or exploit cross product and check for CW/CCW like:

    if (dot(cross(p1-p0,p'-p0),p')>=0)
     if (dot(cross(p2-p1,p'-p1),p')>=0)
      if (dot(cross(p0-p2,p'-p2),p')>=0)
       point_is_inside;
    if (dot(cross(p1-p0,p'-p0),p')<=0)
     if (dot(cross(p2-p1,p'-p1),p')<=0)
      if (dot(cross(p0-p2,p'-p2),p')<=0)
       point_is_inside;
   

   so if all 3 sides of triangle have the same CW/CCW ness to p' you found your triangle.