Instantly share code, notes, and snippets.

# batFINGER/spherical_fibonacci.py

Created August 15, 2016 17:33
Show Gist options
• Save batFINGER/64db074e95b716f839a71882b7efcc50 to your computer and use it in GitHub Desktop.
Spherical Fibonacci Coordinates
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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

### batFINGER commented Sep 20, 2016

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.