# Find triangle which contains point on the sphere

Suppose I have a mesh of the sphere (points on the sphere and a triangulation).

What is a good and efficient way to find the triangle which contains a point on the sphere (the point does not need to be a triangulation node)? It would be nice if the operation could be vectorized, to improve Matlab speed.

The purpose of this is to construct a simple spherical (linear) interpolation operator on the sphere.

• I don't see the difference between this problem and finding a point in a triangle in the plane. Can you elaborate? – Wolfgang Bangerth Oct 9 '15 at 11:36
• I don't have Matlab right now, but have you tried using triangulation.pointlocation directly? Not sure if this works for 3D points but "2D" simplices (triangles). I believe triangulation.cartesiantobarycentric does work in this case (the 3D point is 1st projected onto the triangle's plane), so the pointlocation() function may also ... – GeoMatt22 Oct 9 '15 at 16:51
• @WolfgangBangerth: I thought that the surface problem was more complicated than the plane one and needed other methods. – Beni Bogosel Oct 9 '15 at 20:45
• @GeoMatt22: I'll try to see if I can use the functionalities of the functions you mention. – Beni Bogosel Oct 9 '15 at 20:45
• The sphere problem is as easy as the 2d problem because projecting an arbitrary point onto the sphere is so simple. In essence, you only have to work in polar coordinates, which by and large are two dimensional if you ensure that no vertex lies at either pole. – Wolfgang Bangerth Oct 11 '15 at 1:47

A generic way of looking up the element in which a point of given coordinates may be is to sort them into a quadtree (octree in 3D). The leaves of the tree will contain only element having an intersection with a given range of coordinates (in your case this could be ranges in $(\theta, \phi)$), the number of elements per leaf is kept smaller than a small enough constant $n$ (around 10 to 20), and the test of whether the point belongs to an element is performed on all of these. This gives a cost of order $n + log (N/n)$ for an $N$-element triangulation.
In the case of the sphere there is a simple and computationally efficient method in the case the triangles are not far from being equilateral (thanks for the correction in the comments). Suppose we are working on the unit sphere, so the center of the sphere is at the origin. Then in order to find the triangle which contains a certain point $P$ we construct the exterior normals to each triangle and search the one which makes the minimal angle with the vector $\vec{OP}$. The triangle corresponding to the closest normal will contain the point.
• Besides, the cost here is in $N$, number of elements in the mesh, so you should locate only a constant number of point to ensure that your algorithm has linear cost -- so unsuitable for interpolation, which requires $N$ look-ups. – Joce Feb 8 '16 at 16:34