Created
August 15, 2016 17:33
-
-
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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.