# Clustering points on a sphere

I have written a short Python program which does the following: loads a large data file ($10^9+$ rows) where each row is a point on a sphere. The code then loads a pre-determined triangular grid on the sphere, and counts the number of points in each triangle. I have optimised it as best as I could, however, I'd like to see if it can be optimised even further (at the moment it takes more than 1 hour to go through the entire file).

import numpy as np
import math
import csv
import time
import pandas

PI = 3.141592653589793115997963468544
error = 0.000001

def GalacticToCartesian (starPolar, starCartesian):
starCartesian[:, 0] = np.sin(starPolar[:, 1] + PI/2.0) * np.cos(starPolar[:, 0])
starCartesian[:, 1] = np.sin(starPolar[:, 1] + PI/2.0) * np.sin(starPolar[:, 0])
starCartesian[:, 2] = np.cos(starPolar[:, 1] + PI/2.0)

for coord in np.nditer(starCartesian):
if (np.abs(coord) > error):
coord = 0

def RayTriangleIntersection (star, triangle):
a = np.empty((len(star), 3, 3))
a[..., 0] = star
a[..., 1] = (triangle[1] - triangle[0]) [None, :]
a[..., 2] = (triangle[2] - triangle[0]) [None, :]

b = np.tile(-triangle[0], (len(star), 1))

solution = np.linalg.solve (a, b)

return np.where(np.logical_or.reduce((solution[:, 0] < 0.0, solution[:, 1] < 0.0, solution[:, 2] < 0.0, solution[:, 1] + solution[:, 2] > 1.0)), False, True)

grid = 1

gridFacesFile = "triangles.dat"

triangles = []

for row in csv_file:
triangles.append(np.array([float(elem) for elem in row]).reshape((3,3)))

starCountPerTriangle = np.zeros ((len(triangles)))

dataFile = "data.csv"

chunkSize = 10

data = pandas.read_csv (dataFile, chunksize = chunkSize)

count = 0

t0 = time.clock()
for chunk in data:
currentChunkSize = len(chunk)

print ("Processing stars " + str(count*chunkSize + 1) + " to " + str(count*chunkSize + currentChunkSize) + "...")

count += 1
starsPolar = np.asarray (chunk) [:, 1:3]

starsCartesian = np.zeros ((currentChunkSize, 3))

GalacticToCartesian (starsPolar, starsCartesian)

for (i, triangle) in enumerate(triangles):
belongToCurrent = np.zeros(len(starsCartesian), dtype=bool)

belongToCurrent = RayTriangleIntersection (starsCartesian, triangle)
starCountPerTriangle[i] += np.sum (belongToCurrent)
starsCartesian = starsCartesian[belongToCurrent == False]

print (starCountPerTriangle)
print time.clock()
exit()


Sample contents of triangles.dat:

0.  0.  1.  0.27639320225002106 0.8506508083520399  0.4472135954999579  -0.7236067977499789 0.5257311121191336  0.4472135954999579
0.  0.  1.  -0.7236067977499789 0.5257311121191336  0.4472135954999579  -0.7236067977499789 -0.5257311121191336 0.4472135954999579
0.  0.  1.  -0.7236067977499789 -0.5257311121191336 0.4472135954999579  0.27639320225002106 -0.8506508083520399 0.4472135954999579
0.  0.  1.  0.27639320225002106 -0.8506508083520399 0.4472135954999579  0.8944271909999159  0.  0.4472135954999579
0.  0.  1.  0.8944271909999159  0.  0.4472135954999579  0.27639320225002106 0.8506508083520399  0.4472135954999579
0.7236067977499789  -0.5257311121191336 -0.4472135954999579 -0.27639320225002106    -0.8506508083520399 -0.4472135954999579 0.  0.  -1.
0.7236067977499789  0.5257311121191336  -0.4472135954999579 0.7236067977499789  -0.5257311121191336 -0.4472135954999579 0.  0.  -1.
-0.27639320225002106    0.8506508083520399  -0.4472135954999579 0.7236067977499789  0.5257311121191336  -0.4472135954999579 0.  0.  -1.
-0.8944271909999159 0.  -0.4472135954999579 -0.27639320225002106    0.8506508083520399  -0.4472135954999579 0.  0.  -1.
-0.27639320225002106    -0.8506508083520399 -0.4472135954999579 -0.8944271909999159 0.  -0.4472135954999579 0.  0.  -1.
0.27639320225002106 0.8506508083520399  0.4472135954999579  -0.27639320225002106    0.8506508083520399  -0.4472135954999579 -0.7236067977499789 0.5257311121191336  0.4472135954999579
-0.7236067977499789 0.5257311121191336  0.4472135954999579  -0.8944271909999159 0.  -0.4472135954999579 -0.7236067977499789 -0.5257311121191336 0.4472135954999579
-0.7236067977499789 -0.5257311121191336 0.4472135954999579  -0.27639320225002106    -0.8506508083520399 -0.4472135954999579 0.27639320225002106 -0.8506508083520399 0.4472135954999579
0.27639320225002106 -0.8506508083520399 0.4472135954999579  0.7236067977499789  -0.5257311121191336 -0.4472135954999579 0.8944271909999159  0.  0.4472135954999579
0.8944271909999159  0.  0.4472135954999579  0.7236067977499789  0.5257311121191336  -0.4472135954999579 0.27639320225002106 0.8506508083520399  0.4472135954999579
0.7236067977499789  -0.5257311121191336 -0.4472135954999579 0.27639320225002106 -0.8506508083520399 0.4472135954999579  -0.27639320225002106    -0.8506508083520399 -0.4472135954999579
0.7236067977499789  0.5257311121191336  -0.4472135954999579 0.8944271909999159  0.  0.4472135954999579  0.7236067977499789  -0.5257311121191336 -0.4472135954999579
-0.27639320225002106    0.8506508083520399  -0.4472135954999579 0.27639320225002106 0.8506508083520399  0.4472135954999579  0.7236067977499789  0.5257311121191336  -0.4472135954999579
-0.8944271909999159 0.  -0.4472135954999579 -0.7236067977499789 0.5257311121191336  0.4472135954999579  -0.27639320225002106    0.8506508083520399  -0.4472135954999579
-0.27639320225002106    -0.8506508083520399 -0.4472135954999579 -0.7236067977499789 -0.5257311121191336 0.4472135954999579  -0.8944271909999159 0.  -0.4472135954999579


Sample contents of data.csv

2771204605521258752,354.99796946135444,14.540437148032728
2771204742960125824,354.82484279830356,14.415707483801699
2771204742960126080,354.8298951417335,14.415989960651892
2771204777319869440,354.8029071094679,14.423268507989384
2771204811679612544,354.81845541510194,14.431099973159371
2771204880399095936,354.8457063064526,14.440336041720947
2771204983478149760,354.7888726523609,14.428521435486967
2771204983478301312,354.79834319472315,14.426543847147219
2771204983478303104,354.8005161335831,14.42791145520911
2771205223996496896,354.8109855922654,14.463377776632576
2771205223996502784,354.8090346004138,14.471344809515811
2771205228292116224,354.8161795185497,14.466725484321556
2771205258356230784,354.8605444051401,14.45630798849794
2771205430154941312,354.8868752267852,14.481107665444423
2771205533235029504,354.83772161984615,14.47128996739678
2771205601954482176,354.82237084958234,14.481094852179371
2771205670673114880,354.8670878351698,14.48870385525905
2771205670674023936,354.8595228004882,14.485559575409033
2771205773752341248,354.86085950623266,14.505758545293629
2771205778047955456,354.8609952813116,14.502627279103612
2771205808112038528,354.77186051710794,14.449099660599714
2771205876831516416,354.76023619352054,14.449883454390623
2771205876832294272,354.76299178061487,14.45200820799446
2771205945551001216,354.7889731814567,14.46157162357397
2771205979910751104,354.79879415814753,14.476962176777173
2771206220428932992,354.77291041880505,14.494400659822366
2771206254788674048,354.7811307911569,14.497769928711854
2771206392228463360,354.82139167668066,14.493094428947858
2817884371678369280,349.9513847843209,16.81616842372345
2817884406038111744,349.9166809824806,16.822203915110855
2817884406038113408,349.922753905713,16.823778938852282
2817884440397853312,349.9335912366887,16.827163469580135
2817884440397856896,349.93009817366493,16.832567329127055
2817884474757557376,349.8551142973262,16.772287502999728
2817884509117109888,349.85642566007493,16.77622731670913
2817884509117304064,349.86264368355125,16.784407833125968
2817884612196521600,349.8710862987843,16.790036441368127
2817884650851650816,349.87348368092927,16.803791542316404
2817884680916010368,349.85757395459746,16.80772538367227
2817884783995221120,349.8334640252565,16.801516432255852
2817884783995224576,349.84094679068534,16.807143421574697
2817884818354957568,349.8140522573419,16.799147482474005
2817884852714704512,349.8348205657008,16.811726267485383
2817884887074230656,349.8436346912421,16.814360970847027
2817884887074443264,349.8548374335605,16.81220751898652


In reality, I ended up using chunkSize = 10000000.

I don't pretend to understand all of the code (some comments might help, particularly comments which indicate the types of variables and arguments), but I don't think it's really clustering. In fact, I strongly suspect that it's classifying the star locations into an icosahedral dissection of the sphere. If I'm correct then that points to an algorithm change which seems likely to give a good speed-up.

Linear algebra is, in any case, rather heavyweight. The ray-triangle intersection calculation is unnecessary if all you want is to know whether the ray intersects the triangle. A spherical triangle is the intersection of three hemispheres: testing whether a point is in a hemisphere is a simple dot product and sign test, and testing whether it's in the intersection of three is three dot products and sign tests and two Boolean ANDs. That's worst case sixty dot products per star, and if the ANDs short-circuit then it will be fewer.

However, if I'm correct about the isocahedron then it has a very useful property. It's a Voronoi triangulation. You can classify each star into the corresponding triangle by converting its coordinates to Cartesian, taking the dot product with the centre of each triangle, and picking the triangle which gives the largest dot product. That's twenty dot products per star, but further optimisation would be possible by exploiting the structure of the icosahedron to create a binary decision diagram. I haven't calculated whether it can be reduced to the information-theoretically optimal 4.33 (average) dot products per star, but I certainly expect it could be as few as six or seven.

• Thank you Peter - how would the hemispherical test work if I only have the vertices of the spherical triangles? Other than that, I am not certain that the geodesic dome is a Voronoi grid - do you an easy way to prove this?
– John
Feb 13, 2017 at 11:37
• @John, each edge of the triangle is part of a great circle, and the cross product of the endpoints gives one of its two poles. You can select the correct one by looking at the sign of the dot product with the third vertex of the triangle. As proof that the line between two triangles is equidistant from their centres, it suffices to observe that the icosahedron is a Platonic solid, and so easily has enough symmetry. Feb 13, 2017 at 12:18
• Thank you, Peter. So for this, I need to do 3 cross products, then 3 dot products, then 3 sign checks, and then 3 checks to figure out which hemispheres the star lies in? As for the Voronoi grid - I agree that the icosahedron is one, but higher subdivisions of it are not, and I use these too.
– John
Feb 13, 2017 at 14:15
• @John, for each triangle do three cross products and three dot products with sign check to get poles for the hemisphere checks. Then for each triangle and for each star do (up to) three dot products with sign check. Feb 13, 2017 at 14:26
• thank you very much. How do we know that this is faster than inverting a matrix?
– John
Feb 13, 2017 at 15:01