Spherical Fibonacci Coordinates

spherical_fibonacci.py

import bpy
from numpy import (sum, dot, inf, log, sin, cos, pi,
                   arccos, arctan2, floor, sqrt, trunc,
                   )
from mathutils import Vector, Matrix
import bmesh
from bpy import context

PHI = (1 + sqrt(5)) / 2 # golden ratio

def F(k):    
    val = (PHI**k - (1 - PHI)**k) / sqrt(5)
    return int(round(val))

def clamp(val, minval, maxval):
    return max(min(val, maxval), minval)

def mad(a, b, c):
    return dot(a, b) + c

def madfrac(a, b):
    return mad(a, b, -floor(a * b)) 

class SphericalCoordsPoint:
    @property
    def xyz(self):
        theta = self.theta
        zi = self.zi
        x = cos(theta) * sin(zi)
        y = sin(theta) * sin(zi)
        z = cos(zi)
        return Vector((x,y,z))
    
    def __init__(self, theta, zi):
        self.theta = theta
        self.zi = zi
        #self.xyz = self.point(theta, zi)
        
    def __repr__(self):
        return "SphericalCoord(%.4f, %.4f)" % (self.theta, self.zi)

class SphericalFibonacciPoint(SphericalCoordsPoint):
    @property
    def k(self):
        theta = self.xyz.z
        k = int(max(2, floor(log(n * pi * sqrt(5) * (1 - theta * theta)) / log(PHI * PHI)))) 
        return k    
    def __init__(self, index, n):
        self.index = i = index
        self.theta = 2 * pi * i / PHI
        self.zi = arccos(1 - (2 * i + 1) / n)
        
    def __repr__(self):
        return "%d" % self.i       

class SF():
    '''
    Spherical Fibonacci Coordinates
    '''
    _points = []
    
    @property
    def points(self):
        return [p for p in self._points]
    
    def basis_vectors(self, k):
            bv = []
            n = self.n
            for k in [k, k+1]:
                b = SphericalCoordsPoint(-pow(-1, k) * 2 * pi * pow(PHI, -k),
                         -2 * F(k) / n)
                bv.append(b)
            return bv[0], bv[1]
        
       
    def inverse(self, p):
        n = self.n
        phi = min(arctan2(p.y, p.x), pi)
        costheta = p.z
        
        k = max(2, floor(log(n * pi * sqrt(5) * (1 - costheta * costheta)) / log(PHI * PHI)))
        fk = pow(PHI, k) / sqrt(5)
        f0 = round(fk)
        f1 = round(fk * PHI)
        
        B = Matrix([[2 * pi * madfrac(f0 + 1, PHI -1) - 2 * pi * (PHI - 1),
                     2 * pi * madfrac(f1 + 1, PHI -1) - 2 * pi * (PHI - 1)],
                     [-2 * f0 / n, -2 * f1 / n]])

        invB = B.inverted()
        
        c = invB * Vector([phi, costheta - (1 - 1/n)])
        c = Vector([floor(round(axis, 3)) for axis in c])
        d = inf
        j = 0
        closest_neighbours = []
        for s in range(4):
            costheta = B[1].dot(Vector([s % 2, s // 2]) + c) + (1 - 1 / n)
            costheta = clamp(costheta, -1, 1) * 2 - costheta
            i = floor(n * 0.5 * (1 - costheta))
            phi = 2 * pi * madfrac(i, PHI - 1)
            costheta = 1 - (2 * i + 1) * (1 / n) #rcp(n) = 1 / n # fast reciprocal
            if abs(costheta) > 1:
                costheta = 1
            sintheta = sqrt(1 - costheta * costheta)       
            closest_neighbours.append(int(i))
            q = Vector((cos(phi) * sintheta,
                       sin(phi) * sintheta,
                       costheta))
            squaredistance = (q - p).length_squared
            if squaredistance < d:
                d = squaredistance
                j = i
        return int(j), int(k), closest_neighbours


    def __init__(self, n):
        self.n = n
        for i in range(n):
            self._points.append(SphericalFibonacciPoint(i, n))
        pass

# MODULE ------------------- 
# Make an Spherical Fibonacci Point Set

class SFBMesh():
    _mesh = None
    
    @property    
    def mesh(self):
        if self._mesh is None:
            self._mesh = bpy.data.meshes.new("SFSphere")

        self.bm.to_mesh(self._mesh)
        return self._mesh        
        
    def edge(self):   
        bm = self.bm
        SF = self.SF
        edge_list = list()
        for p in SF.points:
            for bv in SF.basis_vectors(p.k):
                verts = []

                v = SphericalCoordsPoint(p.theta + bv.theta,
                                       p.zi + bv.zi).xyz   

                j, K, cl = SF.inverse(v)
                if j != p.index:
                    e = [p.index, j]
                    e.sort()
                    if e not in edge_list:
                        edge_list.append(e)
                else:
                    print(j, "closest point is self")
                    
            continue
        for e in edge_list:
            bm.edges.new([self.verts[i] for i in e])

        bmesh.ops.contextual_create(bm, geom=bm.edges)

    def __init__(self, n, mesh=None):
        self.n = n
        self.SF = SF(n)
        self.bm = bmesh.new()

        self.verts = [self.bm.verts.new(p.xyz) for p in self.SF.points]
        self._mesh = mesh 

Here's where I got by looking at nearest neighbours as given by inverse method, and restricting edges to 3 and 4 and dropping off extras by either max length or edge angle.