Subdivided icosahedron - how to find the nearest vertex to an arbitrary point

Question

I have an application that creates an approximation to sphere by subdividing an icosahedron. The Cartesian vertex coordinates are converted to spherical coordinates so that all vertices sit on the surface of a unit sphere.

What I need to do next is find the nearest vertex to an arbitrary point on the surface of the sphere. I have come up with two simple algorithms...

1. Brute force search - will be OK for a small number of vertices, but will be excessive for finer subdivisions.

2. Sorted / Indexed search - sort the vertices into some form of order by azimuth and inclination and then create a rough index to speed up a brute force search by limiting its scope.

I was wondering if there was a more subtle, and hopefully higher performing algorithm that I can use instead of one of the two above.

Update 1: I have just recalled that for another part of the application the vertices store information about their neighbours. My new algorithm is

1. Pick an arbitrary start vertex. Find which of its neighbours has a smaller distance to the point to locate. Use this neighbour as the new start vertex. Repeat until none of the vertex's neighbours has a smaller distance to the point. This vertex is the closest to the point.

Answer 1

I have done quite a bit in this area - and will share the repo when I am finished. But for the specific question above, I think a good approach is to use a KDTree.

Here is the set of face-centres of the Dymaxion-based icosahedron, loaded into a KDTree, and then a set of random points are tested on it.

This is a very fast means of establishing the correct triangle of an icosahedron. The following includes a trivial 3D point -> decimal lat/long converter xyz_ll (for sphere normalised projections only) but still reasonably suited for placing points onto google maps.

Because most of these icosahedron mappings do not use a geoid, it's important (when mapping latitude/longitude) to have already normalised the real world values to a unit sphere such as Plate Carrée.

import numpy as np
import scipy as sp
import scipy.spatial
import math


def xyz_ll(p: tuple) -> tuple:
    x, y, z = p
    return math.degrees(math.atan2(z, math.sqrt(x * x + y * y))), math.degrees(math.atan2(y, x))

def ll_xyz(ll: tuple) -> tuple:
    phi, theta = math.radians(ll[0]), math.radians(ll[1])
    x = math.cos(phi) * math.cos(theta)
    y = math.cos(phi) * math.sin(theta)
    z = math.sin(phi)  # z is 'up'
    return x, y, z


if __name__ == '__main__':
# This is a rendering (in unit sphere coordinates)
# of the centres of the 20 faces of a dymaxion icosahedron.
# suitable for equilateral triangles!

    cx = np.array([
        (0.17476898, 0.5231760133333333, 0.5720300666666667),
        (-0.08387563000000002, 0.77832093, 0.13659113999999997),
        (-0.5903144233333334, 0.28559418333333336, -0.44882131666666664),
        (-0.5884910200000001, 0.5302967366666667, 0.0627648033333333),
        (0.6446662, 0.27407261, 0.37518718666666667),
        (0.5903144233333334, -0.28559418333333336, 0.44882131666666664),
        (-0.16999525000000001, 0.11746359000000002, 0.7673197833333334),
        (0.08682595666666666, -0.38238388000000006, 0.69117259),
        (0.5884910200000001, -0.5302967366666667, -0.0627648033333333),
        (-0.22617043, -0.6869057566666666, 0.3293677933333333),
        (0.08387563000000002, -0.77832093, -0.13659113999999997),
        (-0.17476898, -0.5231760133333333, -0.5720300666666667),
        (0.22617043, 0.6869057566666666, -0.3293677933333333),
        (-0.08682595666666666, 0.38238387999999995, -0.69117259),
        (0.6764340433333333, 0.3752631566666666, -0.18190732666666665),
        (0.6417158733333334, -0.12186444, -0.4525765433333333),
        (0.16999525000000001, -0.11746359, -0.7673197833333333),
        (-0.6417158733333334, 0.12186444, 0.4525765433333333),
        (-0.6764340433333333, -0.3752631566666666, 0.18190732666666665),
        (-0.6446662, -0.27407261, -0.37518718666666667)
    ])
    tree = sp.spatial.KDTree(cx)  # this was improved in SciPy 1.6

    # make 100 random unit vectors suitable for a sphere and normalise.
    n = 100
    v = np.random.normal(size=(n, 3))
    vn = [0, 0, 1] + v / np.linalg.norm(v, axis=1, keepdims=True)

    for pt in vn:
        dx, k = tree.query(pt, workers=-1, k=1)
        print(pt, tree.data[k], xyz_ll(tree.data[k]), dx)

While the above only uses 20 face centres, it's easily extended to as many normals as you wish.

In terms of conversions though, it is often better not be singular when calculating points - for example, I determine sub-triangles via midpoint intersection (which is highly reliable in both 2D and 3D for equilateral triangles), and store both corresponding values for each point. This means that there are (1) NO rotations needed, and (2) zero complexity in conversion from the unit sphere to the plane.

In this case, an np.linspace can be used for each dimension - and then one can just use that. eg

w = np.searchsorted(self.lon, lon) % self.width

This returns the index value for that dimension which allows a simple per-dimension lookup, as long as the indices match!

Below is a sample of my current project - which is 'just for fun'. This projects the Nasa 'blue marble' Plate Carrée (8km here)

Answer 2

Scanning through the responses, I think I may be off base, but what you're after is simple. I think.

Since you're dealing with just points that sit on the sphere, you can just drop a line from the vertex to the center of the sphere, drop another line from the arbitrary point to the center and solve for the angle created between them. Smaller is better. The easiest and cheapest way I think would be the dot product. The angle basically falls out of it. Here's a link about it: http://www.kynd.info/library/mathandphysics/dotProduct_01/ For testing them, I would suggest picking a vertex, testing it, then testing its neighbors. It SHOULD always be in the direction of the smallest neighbor (angle should always decrease as you get closer to the vertex you're after)

Anyhow, I hope that's what you're after. Oh, and I came across this page while looking for your subdivision algorithm. Hard to find; if you could post a link to it I think it would help out a lot more than just myself.

Answer 3

If the icosahedron has one vertex at the north pole and the opposite vertex at the south pole then there are 2 groups each of 5 vertices which are in planes parallel to the equator. With a little geometry I figure that these planes are at N/S 57.3056° (decimals, not dd.mmss). This divides your icosahedron into 4 latitude zones;

• anything north (south) of 28.6528° is closest to the vertex at the nearer pole;
• anything between the equator and north (south) 28.6528° is closer to one of the 5 vertices in that zone.

I'm working this as a navigator would, arcs measured in degrees and denoted north and south; if you prefer a more mathematical convention you can translate this all to your version of spherical coordinates quite easily.

I suspect, though I haven't coded it, that checking the distance to 5 vertices and selecting the nearest will be quicker than more sophisticated approaches based on partitioning the surface of the sphere into the projections of the faces of the icosahedron, or projecting the points on the sphere back onto the icosahedron and working the problem in that coordinate system.

For example, the approach you suggest in your update 1 will require the computation of the distance to 6 vertices (the first, arbitrarily chosen one and its 5 neighbours) at least.

It doesn't matter (if you only want to know which vertex is nearest) whether you calculate distances in Cartesian or spherical coordinates. However, calculation in Cartesian coordinates avoids a lot of calls to trigonometric functions.

If, on the other hand, you haven't arranged your icosahedron with vertices at the poles of your sphere, well, you should have !

Answer 4

One of possible solutions is to build BSP tree for vertices: http://en.wikipedia.org/wiki/Binary_space_partitioning