Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Spherical Fibonacci Coordinates
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
Copy link
Author

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.

screenshot from 2016-09-21 03 06 53

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment