From 3b915dd935b65a176554303bf3f12d5c9e9bf1ec Mon Sep 17 00:00:00 2001 From: mutablehurdle Date: Thu, 27 Jan 2022 07:59:36 -0800 Subject: [PATCH 1/9] Update d3.py --- sdf/d3.py | 78 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) diff --git a/sdf/d3.py b/sdf/d3.py index 237c1cd6..d7cd9457 100644 --- a/sdf/d3.py +++ b/sdf/d3.py @@ -280,6 +280,84 @@ def f(p): return np.sqrt((d2 + qz * qz) / m2) * np.sign(_max(qz, -py)) return f +#MetaCORE surfaces + +@sdf3 +def MO(h,slant,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.sin(z)+np.cos(x+slant*np.sin(y)))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +def EB(h,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.cos(x)+np.cos(y)*np.cos(z))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +#Minimal surfaces + +@sdf3 +def schwarzP(h,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.cos(x)+np.cos(y)+np.cos(z))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +def schwarzD(h,size, center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.sin(x)*np.sin(y)*np.sin(z)+np.sin(x)*np.cos(y)*np.cos(z)+np.cos(x)*np.sin(y)*np.cos(z))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + + +@sdf3 +def gyroid(h,t,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-t)-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +#note -- careful with bounds on this one +def scherkSecond(h,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.sin(z)-np.sinh(x)*np.sinh(y))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + # Platonic Solids @sdf3 From 5be0515b74d898d02651e19d1aba0f2dd2aa7470 Mon Sep 17 00:00:00 2001 From: mutablehurdle Date: Fri, 28 Jan 2022 14:28:24 -0800 Subject: [PATCH 2/9] linear graded gyroid graded_gyroid linearly transitions thickness in x direction, linearly transitions level in the y direction. Ultimate goal is to take thickness and level as scalar field inputs. WIP --- sdf/d3.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/sdf/d3.py b/sdf/d3.py index d7cd9457..44d9ed3e 100644 --- a/sdf/d3.py +++ b/sdf/d3.py @@ -345,6 +345,22 @@ def f(p): return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f +@sdf3 +def graded_gyroid(h_min,h_max,t_min,t_max,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + h = h_min+(h_max-h_min)*(x+size[0]/2)/size[0] + t = t_min+(t_max-t_min)*(y+size[1]/2)/size[1] + d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-t)-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + + + @sdf3 #note -- careful with bounds on this one def scherkSecond(h,size,center = ORIGIN): @@ -533,6 +549,7 @@ def f(p): return other(_vec(x, y, z)) return f + @op3 def transition_linear(f0, f1, p0=-Z, p1=Z, e=ease.linear): p0 = np.array(p0) From c14d092dac169949c283c87d0262006be8684d77 Mon Sep 17 00:00:00 2001 From: mutablehurdle Date: Sat, 5 Feb 2022 09:46:08 -0800 Subject: [PATCH 3/9] testing functionally graded gyroid --- build/lib/sdf/__init__.py | 25 ++ build/lib/sdf/d2.py | 286 +++++++++++++++++ build/lib/sdf/d3.py | 640 ++++++++++++++++++++++++++++++++++++++ build/lib/sdf/dn.py | 114 +++++++ build/lib/sdf/ease.py | 188 +++++++++++ build/lib/sdf/mesh.py | 243 +++++++++++++++ build/lib/sdf/progress.py | 82 +++++ build/lib/sdf/stl.py | 24 ++ build/lib/sdf/text.py | 153 +++++++++ build/lib/sdf/util.py | 6 + sdf/d3.py | 16 +- 11 files changed, 1776 insertions(+), 1 deletion(-) create mode 100644 build/lib/sdf/__init__.py create mode 100644 build/lib/sdf/d2.py create mode 100644 build/lib/sdf/d3.py create mode 100644 build/lib/sdf/dn.py create mode 100644 build/lib/sdf/ease.py create mode 100644 build/lib/sdf/mesh.py create mode 100644 build/lib/sdf/progress.py create mode 100644 build/lib/sdf/stl.py create mode 100644 build/lib/sdf/text.py create mode 100644 build/lib/sdf/util.py diff --git a/build/lib/sdf/__init__.py b/build/lib/sdf/__init__.py new file mode 100644 index 00000000..8f6055b5 --- /dev/null +++ b/build/lib/sdf/__init__.py @@ -0,0 +1,25 @@ +from . import d2, d3, ease + +from .util import * + +from .d2 import * + +from .d3 import * + +from .text import ( + measure_image, + measure_text, + image, + text, +) + +from .mesh import ( + generate, + save, + sample_slice, + show_slice, +) + +from .stl import ( + write_binary_stl, +) diff --git a/build/lib/sdf/d2.py b/build/lib/sdf/d2.py new file mode 100644 index 00000000..f16a9d42 --- /dev/null +++ b/build/lib/sdf/d2.py @@ -0,0 +1,286 @@ +import functools +import numpy as np +import operator + +from . import dn, d3, ease + +# Constants + +ORIGIN = np.array((0, 0)) + +X = np.array((1, 0)) +Y = np.array((0, 1)) + +UP = Y + +# SDF Class + +_ops = {} + +class SDF2: + def __init__(self, f): + self.f = f + def __call__(self, p): + return self.f(p).reshape((-1, 1)) + def __getattr__(self, name): + if name in _ops: + f = _ops[name] + return functools.partial(f, self) + raise AttributeError + def __or__(self, other): + return union(self, other) + def __and__(self, other): + return intersection(self, other) + def __sub__(self, other): + return difference(self, other) + def k(self, k=None): + self._k = k + return self + +def sdf2(f): + def wrapper(*args, **kwargs): + return SDF2(f(*args, **kwargs)) + return wrapper + +def op2(f): + def wrapper(*args, **kwargs): + return SDF2(f(*args, **kwargs)) + _ops[f.__name__] = wrapper + return wrapper + +def op23(f): + def wrapper(*args, **kwargs): + return d3.SDF3(f(*args, **kwargs)) + _ops[f.__name__] = wrapper + return wrapper + +# Helpers + +def _length(a): + return np.linalg.norm(a, axis=1) + +def _normalize(a): + return a / np.linalg.norm(a) + +def _dot(a, b): + return np.sum(a * b, axis=1) + +def _vec(*arrs): + return np.stack(arrs, axis=-1) + +_min = np.minimum +_max = np.maximum + +# Primitives + +@sdf2 +def circle(radius=1, center=ORIGIN): + def f(p): + return _length(p - center) - radius + return f + +@sdf2 +def line(normal=UP, point=ORIGIN): + normal = _normalize(normal) + def f(p): + return np.dot(point - p, normal) + return f + +@sdf2 +def slab(x0=None, y0=None, x1=None, y1=None, k=None): + fs = [] + if x0 is not None: + fs.append(line(X, (x0, 0))) + if x1 is not None: + fs.append(line(-X, (x1, 0))) + if y0 is not None: + fs.append(line(Y, (0, y0))) + if y1 is not None: + fs.append(line(-Y, (0, y1))) + return intersection(*fs, k=k) + +@sdf2 +def rectangle(size=1, center=ORIGIN, a=None, b=None): + if a is not None and b is not None: + a = np.array(a) + b = np.array(b) + size = b - a + center = a + size / 2 + return rectangle(size, center) + size = np.array(size) + def f(p): + q = np.abs(p - center) - size / 2 + return _length(_max(q, 0)) + _min(np.amax(q, axis=1), 0) + return f + +@sdf2 +def rounded_rectangle(size, radius, center=ORIGIN): + try: + r0, r1, r2, r3 = radius + except TypeError: + r0 = r1 = r2 = r3 = radius + def f(p): + x = p[:,0] + y = p[:,1] + r = np.zeros(len(p)).reshape((-1, 1)) + r[np.logical_and(x > 0, y > 0)] = r0 + r[np.logical_and(x > 0, y <= 0)] = r1 + r[np.logical_and(x <= 0, y <= 0)] = r2 + r[np.logical_and(x <= 0, y > 0)] = r3 + q = np.abs(p) - size / 2 + r + return ( + _min(_max(q[:,0], q[:,1]), 0).reshape((-1, 1)) + + _length(_max(q, 0)).reshape((-1, 1)) - r) + return f + +@sdf2 +def equilateral_triangle(): + def f(p): + k = 3 ** 0.5 + p = _vec( + np.abs(p[:,0]) - 1, + p[:,1] + 1 / k) + w = p[:,0] + k * p[:,1] > 0 + q = _vec( + p[:,0] - k * p[:,1], + -k * p[:,0] - p[:,1]) / 2 + p = np.where(w.reshape((-1, 1)), q, p) + p = _vec( + p[:,0] - np.clip(p[:,0], -2, 0), + p[:,1]) + return -_length(p) * np.sign(p[:,1]) + return f + +@sdf2 +def hexagon(r): + def f(p): + k = np.array((3 ** 0.5 / -2, 0.5, np.tan(np.pi / 6))) + p = np.abs(p) + p -= 2 * k[:2] * _min(_dot(k[:2], p), 0).reshape((-1, 1)) + p -= _vec( + np.clip(p[:,0], -k[2] * r, k[2] * r), + np.zeros(len(p)) + r) + return _length(p) * np.sign(p[:,1]) + return f + +@sdf2 +def rounded_x(w, r): + def f(p): + p = np.abs(p) + q = (_min(p[:,0] + p[:,1], w) * 0.5).reshape((-1, 1)) + return _length(p - q) - r + return f + +@sdf2 +def polygon(points): + points = [np.array(p) for p in points] + def f(p): + n = len(points) + d = _dot(p - points[0], p - points[0]) + s = np.ones(len(p)) + for i in range(n): + j = (i + n - 1) % n + vi = points[i] + vj = points[j] + e = vj - vi + w = p - vi + b = w - e * np.clip(np.dot(w, e) / np.dot(e, e), 0, 1).reshape((-1, 1)) + d = _min(d, _dot(b, b)) + c1 = p[:,1] >= vi[1] + c2 = p[:,1] < vj[1] + c3 = e[0] * w[:,1] > e[1] * w[:,0] + c = _vec(c1, c2, c3) + s = np.where(np.all(c, axis=1) | np.all(~c, axis=1), -s, s) + return s * np.sqrt(d) + return f + +# Positioning + +@op2 +def translate(other, offset): + def f(p): + return other(p - offset) + return f + +@op2 +def scale(other, factor): + try: + x, y = factor + except TypeError: + x = y = factor + s = (x, y) + m = min(x, y) + def f(p): + return other(p / s) * m + return f + +@op2 +def rotate(other, angle): + s = np.sin(angle) + c = np.cos(angle) + m = 1 - c + matrix = np.array([ + [c, -s], + [s, c], + ]).T + def f(p): + return other(np.dot(p, matrix)) + return f + +@op2 +def circular_array(other, count): + angles = [i / count * 2 * np.pi for i in range(count)] + return union(*[other.rotate(a) for a in angles]) + +# Alterations + +@op2 +def elongate(other, size): + def f(p): + q = np.abs(p) - size + x = q[:,0].reshape((-1, 1)) + y = q[:,1].reshape((-1, 1)) + w = _min(_max(x, y), 0) + return other(_max(q, 0)) + w + return f + +# 2D => 3D Operations + +@op23 +def extrude(other, h): + def f(p): + d = other(p[:,[0,1]]) + w = _vec(d.reshape(-1), np.abs(p[:,2]) - h / 2) + return _min(_max(w[:,0], w[:,1]), 0) + _length(_max(w, 0)) + return f + +@op23 +def extrude_to(a, b, h, e=ease.linear): + def f(p): + d1 = a(p[:,[0,1]]) + d2 = b(p[:,[0,1]]) + t = e(np.clip(p[:,2] / h, -0.5, 0.5) + 0.5) + d = d1 + (d2 - d1) * t.reshape((-1, 1)) + w = _vec(d.reshape(-1), np.abs(p[:,2]) - h / 2) + return _min(_max(w[:,0], w[:,1]), 0) + _length(_max(w, 0)) + return f + +@op23 +def revolve(other, offset=0): + def f(p): + xy = p[:,[0,1]] + q = _vec(_length(xy) - offset, p[:,2]) + return other(q) + return f + +# Common + +union = op2(dn.union) +difference = op2(dn.difference) +intersection = op2(dn.intersection) +blend = op2(dn.blend) +negate = op2(dn.negate) +dilate = op2(dn.dilate) +erode = op2(dn.erode) +shell = op2(dn.shell) +repeat = op2(dn.repeat) diff --git a/build/lib/sdf/d3.py b/build/lib/sdf/d3.py new file mode 100644 index 00000000..991de192 --- /dev/null +++ b/build/lib/sdf/d3.py @@ -0,0 +1,640 @@ +import functools +import numpy as np +import operator + +from . import dn, d2, ease, mesh + +# Constants + +ORIGIN = np.array((0, 0, 0)) + +X = np.array((1, 0, 0)) +Y = np.array((0, 1, 0)) +Z = np.array((0, 0, 1)) + +UP = Z + +# SDF Class + +_ops = {} + +class SDF3: + def __init__(self, f): + self.f = f + def __call__(self, p): + return self.f(p).reshape((-1, 1)) + def __getattr__(self, name): + if name in _ops: + f = _ops[name] + return functools.partial(f, self) + raise AttributeError + def __or__(self, other): + return union(self, other) + def __and__(self, other): + return intersection(self, other) + def __sub__(self, other): + return difference(self, other) + def k(self, k=None): + self._k = k + return self + def generate(self, *args, **kwargs): + return mesh.generate(self, *args, **kwargs) + def save(self, path, *args, **kwargs): + return mesh.save(path, self, *args, **kwargs) + def show_slice(self, *args, **kwargs): + return mesh.show_slice(self, *args, **kwargs) + +def sdf3(f): + def wrapper(*args, **kwargs): + return SDF3(f(*args, **kwargs)) + return wrapper + +def op3(f): + def wrapper(*args, **kwargs): + return SDF3(f(*args, **kwargs)) + _ops[f.__name__] = wrapper + return wrapper + +def op32(f): + def wrapper(*args, **kwargs): + return d2.SDF2(f(*args, **kwargs)) + _ops[f.__name__] = wrapper + return wrapper + +# Helpers + +def _length(a): + return np.linalg.norm(a, axis=1) + +def _normalize(a): + return a / np.linalg.norm(a) + +def _dot(a, b): + return np.sum(a * b, axis=1) + +def _vec(*arrs): + return np.stack(arrs, axis=-1) + +def _perpendicular(v): + if v[1] == 0 and v[2] == 0: + if v[0] == 0: + raise ValueError('zero vector') + else: + return np.cross(v, [0, 1, 0]) + return np.cross(v, [1, 0, 0]) + +_min = np.minimum +_max = np.maximum + +# Primitives + +@sdf3 +def sphere(radius=1, center=ORIGIN): + def f(p): + return _length(p - center) - radius + return f + +@sdf3 +def plane(normal=UP, point=ORIGIN): + normal = _normalize(normal) + def f(p): + return np.dot(point - p, normal) + return f + +@sdf3 +def slab(x0=None, y0=None, z0=None, x1=None, y1=None, z1=None, k=None): + fs = [] + if x0 is not None: + fs.append(plane(X, (x0, 0, 0))) + if x1 is not None: + fs.append(plane(-X, (x1, 0, 0))) + if y0 is not None: + fs.append(plane(Y, (0, y0, 0))) + if y1 is not None: + fs.append(plane(-Y, (0, y1, 0))) + if z0 is not None: + fs.append(plane(Z, (0, 0, z0))) + if z1 is not None: + fs.append(plane(-Z, (0, 0, z1))) + return intersection(*fs, k=k) + +@sdf3 +def box(size=1, center=ORIGIN, a=None, b=None): + if a is not None and b is not None: + a = np.array(a) + b = np.array(b) + size = b - a + center = a + size / 2 + return box(size, center) + size = np.array(size) + def f(p): + q = np.abs(p - center) - size / 2 + return _length(_max(q, 0)) + _min(np.amax(q, axis=1), 0) + return f + +@sdf3 +def rounded_box(size, radius): + size = np.array(size) + def f(p): + q = np.abs(p) - size / 2 + radius + return _length(_max(q, 0)) + _min(np.amax(q, axis=1), 0) - radius + return f + +@sdf3 +def wireframe_box(size, thickness): + size = np.array(size) + def g(a, b, c): + return _length(_max(_vec(a, b, c), 0)) + _min(_max(a, _max(b, c)), 0) + def f(p): + p = np.abs(p) - size / 2 - thickness / 2 + q = np.abs(p + thickness / 2) - thickness / 2 + px, py, pz = p[:,0], p[:,1], p[:,2] + qx, qy, qz = q[:,0], q[:,1], q[:,2] + return _min(_min(g(px, qy, qz), g(qx, py, qz)), g(qx, qy, pz)) + return f + +@sdf3 +def torus(r1, r2): + def f(p): + xy = p[:,[0,1]] + z = p[:,2] + a = _length(xy) - r1 + b = _length(_vec(a, z)) - r2 + return b + return f + +@sdf3 +def capsule(a, b, radius): + a = np.array(a) + b = np.array(b) + def f(p): + pa = p - a + ba = b - a + h = np.clip(np.dot(pa, ba) / np.dot(ba, ba), 0, 1).reshape((-1, 1)) + return _length(pa - np.multiply(ba, h)) - radius + return f + +@sdf3 +def cylinder(radius): + def f(p): + return _length(p[:,[0,1]]) - radius; + return f + +@sdf3 +def capped_cylinder(a, b, radius): + a = np.array(a) + b = np.array(b) + def f(p): + ba = b - a + pa = p - a + baba = np.dot(ba, ba) + paba = np.dot(pa, ba).reshape((-1, 1)) + x = _length(pa * baba - ba * paba) - radius * baba + y = np.abs(paba - baba * 0.5) - baba * 0.5 + x = x.reshape((-1, 1)) + y = y.reshape((-1, 1)) + x2 = x * x + y2 = y * y * baba + d = np.where( + _max(x, y) < 0, + -_min(x2, y2), + np.where(x > 0, x2, 0) + np.where(y > 0, y2, 0)) + return np.sign(d) * np.sqrt(np.abs(d)) / baba + return f + +@sdf3 +def rounded_cylinder(ra, rb, h): + def f(p): + d = _vec( + _length(p[:,[0,1]]) - ra + rb, + np.abs(p[:,2]) - h / 2 + rb) + return ( + _min(_max(d[:,0], d[:,1]), 0) + + _length(_max(d, 0)) - rb) + return f + +@sdf3 +def capped_cone(a, b, ra, rb): + a = np.array(a) + b = np.array(b) + def f(p): + rba = rb - ra + baba = np.dot(b - a, b - a) + papa = _dot(p - a, p - a) + paba = np.dot(p - a, b - a) / baba + x = np.sqrt(papa - paba * paba * baba) + cax = _max(0, x - np.where(paba < 0.5, ra, rb)) + cay = np.abs(paba - 0.5) - 0.5 + k = rba * rba + baba + f = np.clip((rba * (x - ra) + paba * baba) / k, 0, 1) + cbx = x - ra - f * rba + cby = paba - f + s = np.where(np.logical_and(cbx < 0, cay < 0), -1, 1) + return s * np.sqrt(_min( + cax * cax + cay * cay * baba, + cbx * cbx + cby * cby * baba)) + return f + +@sdf3 +def rounded_cone(r1, r2, h): + def f(p): + q = _vec(_length(p[:,[0,1]]), p[:,2]) + b = (r1 - r2) / h + a = np.sqrt(1 - b * b) + k = np.dot(q, _vec(-b, a)) + c1 = _length(q) - r1 + c2 = _length(q - _vec(0, h)) - r2 + c3 = np.dot(q, _vec(a, b)) - r1 + return np.where(k < 0, c1, np.where(k > a * h, c2, c3)) + return f + +@sdf3 +def ellipsoid(size): + size = np.array(size) + def f(p): + k0 = _length(p / size) + k1 = _length(p / (size * size)) + return k0 * (k0 - 1) / k1 + return f + +@sdf3 +def pyramid(h): + def f(p): + a = np.abs(p[:,[0,1]]) - 0.5 + w = a[:,1] > a[:,0] + a[w] = a[:,[1,0]][w] + px = a[:,0] + py = p[:,2] + pz = a[:,1] + m2 = h * h + 0.25 + qx = pz + qy = h * py - 0.5 * px + qz = h * px + 0.5 * py + s = _max(-qx, 0) + t = np.clip((qy - 0.5 * pz) / (m2 + 0.25), 0, 1) + a = m2 * (qx + s) ** 2 + qy * qy + b = m2 * (qx + 0.5 * t) ** 2 + (qy - m2 * t) ** 2 + d2 = np.where( + _min(qy, -qx * m2 - qy * 0.5) > 0, + 0, _min(a, b)) + return np.sqrt((d2 + qz * qz) / m2) * np.sign(_max(qz, -py)) + return f + +#MetaCORE surfaces + +@sdf3 +def MO(h,slant,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.sin(z)+np.cos(x+slant*np.sin(y)))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +def EB(h,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.cos(x)+np.cos(y)*np.cos(z))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +#Minimal surfaces + +@sdf3 +def schwarzP(h,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.cos(x)+np.cos(y)+np.cos(z))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +def schwarzD(h,size, center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.sin(x)*np.sin(y)*np.sin(z)+np.sin(x)*np.cos(y)*np.cos(z)+np.cos(x)*np.sin(y)*np.cos(z))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + + +@sdf3 +def gyroid(h,t,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-t)-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +def FG_gyroid(h,t,size,center = ORIGIN): + size = np.array(size) + def f(p): + tt = t(p) + hh = h(p) + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-tt)-hh + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +def graded_gyroid(h_min,h_max,t_min,t_max,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + h = h_min+(h_max-h_min)*(x+size[0]/2)/size[0] + t = t_min+(t_max-t_min)*(y+size[1]/2)/size[1] + d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-t)-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + + +@sdf3 +#note -- careful with bounds on this one +def scherkSecond(h,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.sin(z)-np.sinh(x)*np.sinh(y))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + + +# Platonic Solids + +@sdf3 +def tetrahedron(r): + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + return (_max(np.abs(x + y) - z, np.abs(x - y) + z) - 1) / np.sqrt(3) + return f + +@sdf3 +def octahedron(r): + def f(p): + return (np.sum(np.abs(p), axis=1) - r) * np.tan(np.radians(30)) + return f + +@sdf3 +def dodecahedron(r): + x, y, z = _normalize(((1 + np.sqrt(5)) / 2, 1, 0)) + def f(p): + p = np.abs(p / r) + a = np.dot(p, (x, y, z)) + b = np.dot(p, (z, x, y)) + c = np.dot(p, (y, z, x)) + q = (_max(_max(a, b), c) - x) * r + return q + return f + +@sdf3 +def icosahedron(r): + r *= 0.8506507174597755 + x, y, z = _normalize(((np.sqrt(5) + 3) / 2, 1, 0)) + w = np.sqrt(3) / 3 + def f(p): + p = np.abs(p / r) + a = np.dot(p, (x, y, z)) + b = np.dot(p, (z, x, y)) + c = np.dot(p, (y, z, x)) + d = np.dot(p, (w, w, w)) - x + return _max(_max(_max(a, b), c) - x, d) * r + return f + +# Positioning + +@op3 +def translate(other, offset): + def f(p): + return other(p - offset) + return f + +@op3 +def scale(other, factor): + try: + x, y, z = factor + except TypeError: + x = y = z = factor + s = (x, y, z) + m = min(x, min(y, z)) + def f(p): + return other(p / s) * m + return f + +@op3 +def rotate(other, angle, vector=Z): + x, y, z = _normalize(vector) + s = np.sin(angle) + c = np.cos(angle) + m = 1 - c + matrix = np.array([ + [m*x*x + c, m*x*y + z*s, m*z*x - y*s], + [m*x*y - z*s, m*y*y + c, m*y*z + x*s], + [m*z*x + y*s, m*y*z - x*s, m*z*z + c], + ]).T + def f(p): + return other(np.dot(p, matrix)) + return f + +@op3 +def rotate_to(other, a, b): + a = _normalize(np.array(a)) + b = _normalize(np.array(b)) + dot = np.dot(b, a) + if dot == 1: + return other + if dot == -1: + return rotate(other, np.pi, _perpendicular(a)) + angle = np.arccos(dot) + v = _normalize(np.cross(b, a)) + return rotate(other, angle, v) + +@op3 +def orient(other, axis): + return rotate_to(other, UP, axis) + +@op3 +def circular_array(other, count, offset=0): + other = other.translate(X * offset) + da = 2 * np.pi / count + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.hypot(x, y) + a = np.arctan2(y, x) % da + d1 = other(_vec(np.cos(a - da) * d, np.sin(a - da) * d, z)) + d2 = other(_vec(np.cos(a) * d, np.sin(a) * d, z)) + return _min(d1, d2) + return f + +# Alterations + +@op3 +def elongate(other, size): + def f(p): + q = np.abs(p) - size + x = q[:,0].reshape((-1, 1)) + y = q[:,1].reshape((-1, 1)) + z = q[:,2].reshape((-1, 1)) + w = _min(_max(x, _max(y, z)), 0) + return other(_max(q, 0)) + w + return f + +@op3 +def twist(other, k): + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + c = np.cos(k * z) + s = np.sin(k * z) + x2 = c * x - s * y + y2 = s * x + c * y + z2 = z + return other(_vec(x2, y2, z2)) + return f + +@op3 +def bend(other, k): + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + c = np.cos(k * x) + s = np.sin(k * x) + x2 = c * x - s * y + y2 = s * x + c * y + z2 = z + return other(_vec(x2, y2, z2)) + return f + +@op3 +def bend_linear(other, p0, p1, v, e=ease.linear): + p0 = np.array(p0) + p1 = np.array(p1) + v = -np.array(v) + ab = p1 - p0 + def f(p): + t = np.clip(np.dot(p - p0, ab) / np.dot(ab, ab), 0, 1) + t = e(t).reshape((-1, 1)) + return other(p + t * v) + return f + +@op3 +def bend_radial(other, r0, r1, dz, e=ease.linear): + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + r = np.hypot(x, y) + t = np.clip((r - r0) / (r1 - r0), 0, 1) + z = z - dz * e(t) + return other(_vec(x, y, z)) + return f + + +@op3 +def transition_linear(f0, f1, p0=-Z, p1=Z, e=ease.linear): + p0 = np.array(p0) + p1 = np.array(p1) + ab = p1 - p0 + def f(p): + d1 = f0(p) + d2 = f1(p) + t = np.clip(np.dot(p - p0, ab) / np.dot(ab, ab), 0, 1) + t = e(t).reshape((-1, 1)) + return t * d2 + (1 - t) * d1 + return f + +@op3 +def transition_radial(f0, f1, r0=0, r1=1, e=ease.linear): + def f(p): + d1 = f0(p) + d2 = f1(p) + r = np.hypot(p[:,0], p[:,1]) + t = np.clip((r - r0) / (r1 - r0), 0, 1) + t = e(t).reshape((-1, 1)) + return t * d2 + (1 - t) * d1 + return f + +@op3 +def wrap_around(other, x0, x1, r=None, e=ease.linear): + p0 = X * x0 + p1 = X * x1 + v = -Y + if r is None: + r = np.linalg.norm(p1 - p0) / (2 * np.pi) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.hypot(x, y) - r + d = d.reshape((-1, 1)) + a = np.arctan2(y, x) + t = (a + np.pi) / (2 * np.pi) + t = e(t).reshape((-1, 1)) + q = p0 + (p1 - p0) * t + v * d + q[:,2] = z + return other(q) + return f + +# 3D => 2D Operations + +@op32 +def slice(other): + # TODO: support specifying a slice plane + # TODO: probably a better way to do this + s = slab(z0=-1e-9, z1=1e-9) + a = other & s + b = other.negate() & s + def f(p): + p = _vec(p[:,0], p[:,1], np.zeros(len(p))) + A = a(p).reshape(-1) + B = -b(p).reshape(-1) + w = A <= 0 + A[w] = B[w] + return A + return f + +# Common + +union = op3(dn.union) +difference = op3(dn.difference) +intersection = op3(dn.intersection) +blend = op3(dn.blend) +negate = op3(dn.negate) +dilate = op3(dn.dilate) +erode = op3(dn.erode) +shell = op3(dn.shell) +repeat = op3(dn.repeat) diff --git a/build/lib/sdf/dn.py b/build/lib/sdf/dn.py new file mode 100644 index 00000000..9c245e82 --- /dev/null +++ b/build/lib/sdf/dn.py @@ -0,0 +1,114 @@ +import itertools +import numpy as np + +_min = np.minimum +_max = np.maximum + +def union(a, *bs, k=None): + def f(p): + d1 = a(p) + for b in bs: + d2 = b(p) + K = k or getattr(b, '_k', None) + if K is None: + d1 = _min(d1, d2) + else: + h = np.clip(0.5 + 0.5 * (d2 - d1) / K, 0, 1) + m = d2 + (d1 - d2) * h + d1 = m - K * h * (1 - h) + return d1 + return f + +def difference(a, *bs, k=None): + def f(p): + d1 = a(p) + for b in bs: + d2 = b(p) + K = k or getattr(b, '_k', None) + if K is None: + d1 = _max(d1, -d2) + else: + h = np.clip(0.5 - 0.5 * (d2 + d1) / K, 0, 1) + m = d1 + (-d2 - d1) * h + d1 = m + K * h * (1 - h) + return d1 + return f + +def intersection(a, *bs, k=None): + def f(p): + d1 = a(p) + for b in bs: + d2 = b(p) + K = k or getattr(b, '_k', None) + if K is None: + d1 = _max(d1, d2) + else: + h = np.clip(0.5 - 0.5 * (d2 - d1) / K, 0, 1) + m = d2 + (d1 - d2) * h + d1 = m + K * h * (1 - h) + return d1 + return f + +def blend(a, *bs, k=0.5): + def f(p): + d1 = a(p) + for b in bs: + d2 = b(p) + K = k or getattr(b, '_k', None) + d1 = K * d2 + (1 - K) * d1 + return d1 + return f + +def negate(other): + def f(p): + return -other(p) + return f + +def dilate(other, r): + def f(p): + return other(p) - r + return f + +def erode(other, r): + def f(p): + return other(p) + r + return f + +def shell(other, thickness): + def f(p): + return np.abs(other(p)) - thickness / 2 + return f + +def repeat(other, spacing, count=None, padding=0): + count = np.array(count) if count is not None else None + spacing = np.array(spacing) + + def neighbors(dim, padding, spacing): + try: + padding = [padding[i] for i in range(dim)] + except Exception: + padding = [padding] * dim + try: + spacing = [spacing[i] for i in range(dim)] + except Exception: + spacing = [spacing] * dim + for i, s in enumerate(spacing): + if s == 0: + padding[i] = 0 + axes = [list(range(-p, p + 1)) for p in padding] + return list(itertools.product(*axes)) + + def f(p): + q = np.divide(p, spacing, out=np.zeros_like(p), where=spacing != 0) + if count is None: + index = np.round(q) + else: + index = np.clip(np.round(q), -count, count) + + indexes = [index + n for n in neighbors(p.shape[-1], padding, spacing)] + A = [other(p - spacing * i) for i in indexes] + a = A[0] + for b in A[1:]: + a = _min(a, b) + return a + return f diff --git a/build/lib/sdf/ease.py b/build/lib/sdf/ease.py new file mode 100644 index 00000000..9c7fc5ab --- /dev/null +++ b/build/lib/sdf/ease.py @@ -0,0 +1,188 @@ +import numpy as np + +def linear(t): + return t + +def in_quad(t): + return t * t + +def out_quad(t): + return -t * (t - 2) + +def in_out_quad(t): + u = 2 * t - 1 + a = 2 * t * t + b = -0.5 * (u * (u - 2) - 1) + return np.where(t < 0.5, a, b) + +def in_cubic(t): + return t * t * t + +def out_cubic(t): + u = t - 1 + return u * u * u + 1 + +def in_out_cubic(t): + u = t * 2 + v = u - 2 + a = 0.5 * u * u * u + b = 0.5 * (v * v * v + 2) + return np.where(u < 1, a, b) + +def in_quart(t): + return t * t * t * t + +def out_quart(t): + u = t - 1 + return -(u * u * u * u - 1) + +def in_out_quart(t): + u = t * 2 + v = u - 2 + a = 0.5 * u * u * u * u + b = -0.5 * (v * v * v * v - 2) + return np.where(u < 1, a, b) + +def in_quint(t): + return t * t * t * t * t + +def out_quint(t): + u = t - 1 + return u * u * u * u * u + 1 + +def in_out_quint(t): + u = t * 2 + v = u - 2 + a = 0.5 * u * u * u * u * u + b = 0.5 * (v * v * v * v * v + 2) + return np.where(u < 1, a, b) + +def in_sine(t): + return -np.cos(t * np.pi / 2) + 1 + +def out_sine(t): + return np.sin(t * np.pi / 2) + +def in_out_sine(t): + return -0.5 * (np.cos(np.pi * t) - 1) + +def in_expo(t): + a = np.zeros(len(t)) + b = 2 ** (10 * (t - 1)) + return np.where(t == 0, a, b) + +def out_expo(t): + a = np.zeros(len(t)) + 1 + b = 1 - 2 ** (-10 * t) + return np.where(t == 1, a, b) + +def in_out_expo(t): + zero = np.zeros(len(t)) + one = zero + 1 + a = 0.5 * 2 ** (20 * t - 10) + b = 1 - 0.5 * 2 ** (-20 * t + 10) + return np.where(t == 0, zero, np.where(t == 1, one, np.where(t < 0.5, a, b))) + +def in_circ(t): + return -1 * (np.sqrt(1 - t * t) - 1) + +def out_circ(t): + u = t - 1 + return np.sqrt(1 - u * u) + +def in_out_circ(t): + u = t * 2 + v = u - 2 + a = -0.5 * (np.sqrt(1 - u * u) - 1) + b = 0.5 * (np.sqrt(1 - v * v) + 1) + return np.where(u < 1, a, b) + +def in_elastic(t, k=0.5): + u = t - 1 + return -1 * (2 ** (10 * u) * np.sin((u - k / 4) * (2 * np.pi) / k)) + +def out_elastic(t, k=0.5): + return 2 ** (-10 * t) * np.sin((t - k / 4) * (2 * np.pi / k)) + 1 + +def in_out_elastic(t, k=0.5): + u = t * 2 + v = u - 1 + a = -0.5 * (2 ** (10 * v) * np.sin((v - k / 4) * 2 * np.pi / k)) + b = 2 ** (-10 * v) * np.sin((v - k / 4) * 2 * np.pi / k) * 0.5 + 1 + return np.where(u < 1, a, b) + +def in_back(t): + k = 1.70158 + return t * t * ((k + 1) * t - k) + +def out_back(t): + k = 1.70158 + u = t - 1 + return u * u * ((k + 1) * u + k) + 1 + +def in_out_back(t): + k = 1.70158 * 1.525 + u = t * 2 + v = u - 2 + a = 0.5 * (u * u * ((k + 1) * u - k)) + b = 0.5 * (v * v * ((k + 1) * v + k) + 2) + return np.where(u < 1, a, b) + +def in_bounce(t): + return 1 - out_bounce(1 - t) + +def out_bounce(t): + a = (121 * t * t) / 16 + b = (363 / 40 * t * t) - (99 / 10 * t) + 17 / 5 + c = (4356 / 361 * t * t) - (35442 / 1805 * t) + 16061 / 1805 + d = (54 / 5 * t * t) - (513 / 25 * t) + 268 / 25 + return np.where( + t < 4 / 11, a, np.where( + t < 8 / 11, b, np.where( + t < 9 / 10, c, d))) + +def in_out_bounce(t): + a = in_bounce(2 * t) * 0.5 + b = out_bounce(2 * t - 1) * 0.5 + 0.5 + return np.where(t < 0.5, a, b) + +def in_square(t): + a = np.zeros(len(t)) + b = a + 1 + return np.where(t < 1, a, b) + +def out_square(t): + a = np.zeros(len(t)) + b = a + 1 + return np.where(t > 0, b, a) + +def in_out_square(t): + a = np.zeros(len(t)) + b = a + 1 + return np.where(t < 0.5, a, b) + +def _main(): + import matplotlib.pyplot as plt + fs = [ + linear, + in_quad, out_quad, in_out_quad, + in_cubic, out_cubic, in_out_cubic, + in_quart, out_quart, in_out_quart, + in_quint, out_quint, in_out_quint, + in_sine, out_sine, in_out_sine, + in_expo, out_expo, in_out_expo, + in_circ, out_circ, in_out_circ, + in_elastic, out_elastic, in_out_elastic, + in_back, out_back, in_out_back, + in_bounce, out_bounce, in_out_bounce, + in_square, out_square, in_out_square, + ] + x = np.linspace(0, 1, 1000) + for f in fs: + y = f(x) + plt.plot(x, y, label=f.__name__) + plt.legend() + plt.show() + +if __name__ == '__main__': + _main() diff --git a/build/lib/sdf/mesh.py b/build/lib/sdf/mesh.py new file mode 100644 index 00000000..8895d23c --- /dev/null +++ b/build/lib/sdf/mesh.py @@ -0,0 +1,243 @@ +from functools import partial +from multiprocessing.pool import ThreadPool +from skimage import measure + +import multiprocessing +import itertools +import numpy as np +import time + +from . import progress, stl + +WORKERS = multiprocessing.cpu_count() +SAMPLES = 2 ** 22 +BATCH_SIZE = 32 + +def _marching_cubes(volume, level=0): + verts, faces, _, _ = measure.marching_cubes(volume, level) + return verts[faces].reshape((-1, 3)) + +def _cartesian_product(*arrays): + la = len(arrays) + dtype = np.result_type(*arrays) + arr = np.empty([len(a) for a in arrays] + [la], dtype=dtype) + for i, a in enumerate(np.ix_(*arrays)): + arr[...,i] = a + return arr.reshape(-1, la) + +def _skip(sdf, job): + X, Y, Z = job + x0, x1 = X[0], X[-1] + y0, y1 = Y[0], Y[-1] + z0, z1 = Z[0], Z[-1] + x = (x0 + x1) / 2 + y = (y0 + y1) / 2 + z = (z0 + z1) / 2 + r = abs(sdf(np.array([(x, y, z)])).reshape(-1)[0]) + d = np.linalg.norm(np.array((x-x0, y-y0, z-z0))) + if r <= d: + return False + corners = np.array(list(itertools.product((x0, x1), (y0, y1), (z0, z1)))) + values = sdf(corners).reshape(-1) + same = np.all(values > 0) if values[0] > 0 else np.all(values < 0) + return same + +def _worker(sdf, job, sparse): + X, Y, Z = job + if sparse and _skip(sdf, job): + return None + # return _debug_triangles(X, Y, Z) + P = _cartesian_product(X, Y, Z) + volume = sdf(P).reshape((len(X), len(Y), len(Z))) + try: + points = _marching_cubes(volume) + except Exception: + return [] + # return _debug_triangles(X, Y, Z) + scale = np.array([X[1] - X[0], Y[1] - Y[0], Z[1] - Z[0]]) + offset = np.array([X[0], Y[0], Z[0]]) + return points * scale + offset + +def _estimate_bounds(sdf): + # TODO: raise exception if bound estimation fails + s = 16 + x0 = y0 = z0 = -1e9 + x1 = y1 = z1 = 1e9 + prev = None + for i in range(32): + X = np.linspace(x0, x1, s) + Y = np.linspace(y0, y1, s) + Z = np.linspace(z0, z1, s) + d = np.array([X[1] - X[0], Y[1] - Y[0], Z[1] - Z[0]]) + threshold = np.linalg.norm(d) / 2 + if threshold == prev: + break + prev = threshold + P = _cartesian_product(X, Y, Z) + volume = sdf(P).reshape((len(X), len(Y), len(Z))) + where = np.argwhere(np.abs(volume) <= threshold) + x1, y1, z1 = (x0, y0, z0) + where.max(axis=0) * d + d / 2 + x0, y0, z0 = (x0, y0, z0) + where.min(axis=0) * d - d / 2 + return ((x0, y0, z0), (x1, y1, z1)) + +def generate( + sdf, + step=None, bounds=None, samples=SAMPLES, + workers=WORKERS, batch_size=BATCH_SIZE, + verbose=True, sparse=True): + + start = time.time() + + if bounds is None: + bounds = _estimate_bounds(sdf) + (x0, y0, z0), (x1, y1, z1) = bounds + + if step is None and samples is not None: + volume = (x1 - x0) * (y1 - y0) * (z1 - z0) + step = (volume / samples) ** (1 / 3) + + try: + dx, dy, dz = step + except TypeError: + dx = dy = dz = step + + if verbose: + print('min %g, %g, %g' % (x0, y0, z0)) + print('max %g, %g, %g' % (x1, y1, z1)) + print('step %g, %g, %g' % (dx, dy, dz)) + + X = np.arange(x0, x1, dx) + Y = np.arange(y0, y1, dy) + Z = np.arange(z0, z1, dz) + + s = batch_size + Xs = [X[i:i+s+1] for i in range(0, len(X), s)] + Ys = [Y[i:i+s+1] for i in range(0, len(Y), s)] + Zs = [Z[i:i+s+1] for i in range(0, len(Z), s)] + + batches = list(itertools.product(Xs, Ys, Zs)) + num_batches = len(batches) + num_samples = sum(len(xs) * len(ys) * len(zs) + for xs, ys, zs in batches) + + if verbose: + print('%d samples in %d batches with %d workers' % + (num_samples, num_batches, workers)) + + points = [] + skipped = empty = nonempty = 0 + bar = progress.Bar(num_batches, enabled=verbose) + pool = ThreadPool(workers) + f = partial(_worker, sdf, sparse=sparse) + for result in pool.imap(f, batches): + bar.increment(1) + if result is None: + skipped += 1 + elif len(result) == 0: + empty += 1 + else: + nonempty += 1 + points.extend(result) + bar.done() + + if verbose: + print('%d skipped, %d empty, %d nonempty' % (skipped, empty, nonempty)) + triangles = len(points) // 3 + seconds = time.time() - start + print('%d triangles in %g seconds' % (triangles, seconds)) + + return points + +def save(path, *args, **kwargs): + points = generate(*args, **kwargs) + if path.lower().endswith('.stl'): + stl.write_binary_stl(path, points) + else: + mesh = _mesh(points) + mesh.write(path) + +def _mesh(points): + import meshio + points, cells = np.unique(points, axis=0, return_inverse=True) + cells = [('triangle', cells.reshape((-1, 3)))] + return meshio.Mesh(points, cells) + +def _debug_triangles(X, Y, Z): + x0, x1 = X[0], X[-1] + y0, y1 = Y[0], Y[-1] + z0, z1 = Z[0], Z[-1] + + p = 0.25 + x0, x1 = x0 + (x1 - x0) * p, x1 - (x1 - x0) * p + y0, y1 = y0 + (y1 - y0) * p, y1 - (y1 - y0) * p + z0, z1 = z0 + (z1 - z0) * p, z1 - (z1 - z0) * p + + v = [ + (x0, y0, z0), + (x0, y0, z1), + (x0, y1, z0), + (x0, y1, z1), + (x1, y0, z0), + (x1, y0, z1), + (x1, y1, z0), + (x1, y1, z1), + ] + + return [ + v[3], v[5], v[7], + v[5], v[3], v[1], + v[0], v[6], v[4], + v[6], v[0], v[2], + v[0], v[5], v[1], + v[5], v[0], v[4], + v[5], v[6], v[7], + v[6], v[5], v[4], + v[6], v[3], v[7], + v[3], v[6], v[2], + v[0], v[3], v[2], + v[3], v[0], v[1], + ] + +def sample_slice( + sdf, w=1024, h=1024, + x=None, y=None, z=None, bounds=None): + + if bounds is None: + bounds = _estimate_bounds(sdf) + (x0, y0, z0), (x1, y1, z1) = bounds + + if x is not None: + X = np.array([x]) + Y = np.linspace(y0, y1, w) + Z = np.linspace(z0, z1, h) + extent = (Z[0], Z[-1], Y[0], Y[-1]) + axes = 'ZY' + elif y is not None: + Y = np.array([y]) + X = np.linspace(x0, x1, w) + Z = np.linspace(z0, z1, h) + extent = (Z[0], Z[-1], X[0], X[-1]) + axes = 'ZX' + elif z is not None: + Z = np.array([z]) + X = np.linspace(x0, x1, w) + Y = np.linspace(y0, y1, h) + extent = (Y[0], Y[-1], X[0], X[-1]) + axes = 'YX' + else: + raise Exception('x, y, or z position must be specified') + + P = _cartesian_product(X, Y, Z) + return sdf(P).reshape((w, h)), extent, axes + +def show_slice(*args, **kwargs): + import matplotlib.pyplot as plt + show_abs = kwargs.pop('abs', False) + a, extent, axes = sample_slice(*args, **kwargs) + if show_abs: + a = np.abs(a) + im = plt.imshow(a, extent=extent, origin='lower') + plt.xlabel(axes[0]) + plt.ylabel(axes[1]) + plt.colorbar(im) + plt.show() diff --git a/build/lib/sdf/progress.py b/build/lib/sdf/progress.py new file mode 100644 index 00000000..54ad1610 --- /dev/null +++ b/build/lib/sdf/progress.py @@ -0,0 +1,82 @@ +import sys +import time + +def pretty_time(seconds): + seconds = int(round(seconds)) + s = seconds % 60 + m = (seconds // 60) % 60 + h = (seconds // 3600) + return '%d:%02d:%02d' % (h, m, s) + +class Bar(object): + + def __init__(self, max_value=100, min_value=0, enabled=True): + self.min_value = min_value + self.max_value = max_value + self.value = min_value + self.start_time = time.time() + self.enabled = enabled + + @property + def percent_complete(self): + t = (self.value - self.min_value) / (self.max_value - self.min_value) + return t * 100 + + @property + def elapsed_time(self): + return time.time() - self.start_time + + @property + def eta(self): + t = self.percent_complete / 100 + if t == 0: + return 0 + return (1 - t) * self.elapsed_time / t + + def increment(self, delta): + self.update(self.value + delta) + + def update(self, value): + self.value = value + if self.enabled: + sys.stdout.write(' %s \r' % self.render()) + sys.stdout.flush() + + def done(self): + self.update(self.max_value) + self.stop() + + def stop(self): + if self.enabled: + sys.stdout.write('\n') + sys.stdout.flush() + + def render(self): + items = [ + self.render_percent_complete(), + self.render_value(), + self.render_bar(), + self.render_elapsed_time(), + self.render_eta(), + ] + return ' '.join(items) + + def render_percent_complete(self): + return '%3.0f%%' % self.percent_complete + + def render_value(self): + if self.min_value == 0: + return '(%g of %g)' % (self.value, self.max_value) + else: + return '(%g)' % (self.value) + + def render_bar(self, size=30): + a = int(round(self.percent_complete / 100.0 * size)) + b = size - a + return '[' + '#' * a + '-' * b + ']' + + def render_elapsed_time(self): + return pretty_time(self.elapsed_time) + + def render_eta(self): + return pretty_time(self.eta) diff --git a/build/lib/sdf/stl.py b/build/lib/sdf/stl.py new file mode 100644 index 00000000..f8e733dd --- /dev/null +++ b/build/lib/sdf/stl.py @@ -0,0 +1,24 @@ +import numpy as np +import struct + +def write_binary_stl(path, points): + n = len(points) // 3 + + points = np.array(points, dtype='float32').reshape((-1, 3, 3)) + normals = np.cross(points[:,1] - points[:,0], points[:,2] - points[:,0]) + normals /= np.linalg.norm(normals, axis=1).reshape((-1, 1)) + + dtype = np.dtype([ + ('normal', ('= tw-1) | (j < 0) | (j >= th-1) + d[outside] = q[outside] + return d + + return f + +def _bilinear_interpolate(a, x, y): + x0 = np.floor(x).astype(int) + x1 = x0 + 1 + y0 = np.floor(y).astype(int) + y1 = y0 + 1 + + x0 = np.clip(x0, 0, a.shape[1] - 1) + x1 = np.clip(x1, 0, a.shape[1] - 1) + y0 = np.clip(y0, 0, a.shape[0] - 1) + y1 = np.clip(y1, 0, a.shape[0] - 1) + + pa = a[y0, x0] + pb = a[y1, x0] + pc = a[y0, x1] + pd = a[y1, x1] + + wa = (x1 - x) * (y1 - y) + wb = (x1 - x) * (y - y0) + wc = (x - x0) * (y1 - y) + wd = (x - x0) * (y - y0) + + return wa * pa + wb * pb + wc * pc + wd * pd diff --git a/build/lib/sdf/util.py b/build/lib/sdf/util.py new file mode 100644 index 00000000..2d17b778 --- /dev/null +++ b/build/lib/sdf/util.py @@ -0,0 +1,6 @@ +import math + +pi = math.pi + +degrees = math.degrees +radians = math.radians diff --git a/sdf/d3.py b/sdf/d3.py index 44d9ed3e..991de192 100644 --- a/sdf/d3.py +++ b/sdf/d3.py @@ -345,6 +345,20 @@ def f(p): return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f +@sdf3 +def FG_gyroid(h,t,size,center = ORIGIN): + size = np.array(size) + def f(p): + tt = t(p) + hh = h(p) + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-tt)-hh + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + @sdf3 def graded_gyroid(h_min,h_max,t_min,t_max,size,center = ORIGIN): size = np.array(size) @@ -360,7 +374,6 @@ def f(p): return f - @sdf3 #note -- careful with bounds on this one def scherkSecond(h,size,center = ORIGIN): @@ -374,6 +387,7 @@ def f(p): return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f + # Platonic Solids @sdf3 From 8f25aec4a9c3bf1582c19d951269afef0920b333 Mon Sep 17 00:00:00 2001 From: mutablehurdle <80230725+mutablehurdle@users.noreply.github.com> Date: Thu, 3 Mar 2022 13:26:47 -0800 Subject: [PATCH 4/9] Update d3.py made some length changes to FG_gyroid --- sdf/d3.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sdf/d3.py b/sdf/d3.py index 991de192..ad3c80e3 100644 --- a/sdf/d3.py +++ b/sdf/d3.py @@ -349,8 +349,8 @@ def f(p): def FG_gyroid(h,t,size,center = ORIGIN): size = np.array(size) def f(p): - tt = t(p) - hh = h(p) + tt = _length(t(p)) + hh = _length(h(p)) x = p[:,0] y = p[:,1] z = p[:,2] From 4c0cc4cb7c618ebcaaa7c4052e5b864e52b28115 Mon Sep 17 00:00:00 2001 From: mutablehurdle <80230725+mutablehurdle@users.noreply.github.com> Date: Thu, 15 Sep 2022 08:08:34 -0700 Subject: [PATCH 5/9] Update d3.py --- sdf/d3.py | 58 +++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 54 insertions(+), 4 deletions(-) diff --git a/sdf/d3.py b/sdf/d3.py index ad3c80e3..108f773f 100644 --- a/sdf/d3.py +++ b/sdf/d3.py @@ -345,20 +345,48 @@ def f(p): return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f +# @sdf3 +# def FG_gyroid(h,t,size,center = ORIGIN): +# size = np.array(size) +# def f(p): +# tt = _length(t(p)) +# hh = _length(h(p)) +# x = p[:,0] +# y = p[:,1] +# z = p[:,2] +# d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-tt)-hh +# q = np.abs(p - center) - size / 2 +# return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) +# return f + @sdf3 -def FG_gyroid(h,t,size,center = ORIGIN): +def FG_gyroid(h_min,h_max,fh,t_min,t_max,ft,size,k=1.0,center = ORIGIN, e = ease.linear): size = np.array(size) def f(p): - tt = _length(t(p)) - hh = _length(h(p)) + Gh = _length(fh(p)) + Gt = _length(ft(p)) x = p[:,0] y = p[:,1] z = p[:,2] - d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-tt)-hh + h = h_min+(h_max-h_min)*np.clip(Gh,0,1) + h = e(h).reshape((-1, 1)) + t = t_min+(t_max-t_min)*np.clip(Gt,0,1) + t = e(t).reshape((-1, 1)) + d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-t)-h q = np.abs(p - center) - size / 2 return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f +# def transition_general(f0, f1, f2, k=1.0, e=ease.linear): +# def f(p): +# d1 = f0(p) +# d2 = f1(p) +# G = f2(p) +# t = np.clip(1./(1.+np.exp(k*G)), 0, 1) +# t = e(t).reshape((-1, 1)) +# return t * d2 + (1 - t) * d1 +# return f + @sdf3 def graded_gyroid(h_min,h_max,t_min,t_max,size,center = ORIGIN): size = np.array(size) @@ -588,6 +616,28 @@ def f(p): return t * d2 + (1 - t) * d1 return f +@op3 +def transition_sigmoid(f0, f1, k=1.0, e=ease.linear): + def f(p): + d1 = f0(p) + d2 = f1(p) + G = p[:,2] + t = np.clip(1./(1.+np.exp(k*G)), 0, 1) + t = e(t).reshape((-1, 1)) + return t * d2 + (1 - t) * d1 + return f + +@op3 +def transition_general(f0, f1, f2, k=1.0, e=ease.linear): + def f(p): + d1 = f0(p) + d2 = f1(p) + G = f2(p) + t = np.clip(1./(1.+np.exp(k*G)), 0, 1) + t = e(t).reshape((-1, 1)) + return t * d2 + (1 - t) * d1 + return f + @op3 def wrap_around(other, x0, x1, r=None, e=ease.linear): p0 = X * x0 From 7ab68e3336e98b4eba5457155f1d7badfcd82f9c Mon Sep 17 00:00:00 2001 From: mutablehurdle Date: Thu, 15 Sep 2022 08:09:03 -0700 Subject: [PATCH 6/9] Update d3.py --- sdf/d3.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/sdf/d3.py b/sdf/d3.py index ad3c80e3..7e5deb11 100644 --- a/sdf/d3.py +++ b/sdf/d3.py @@ -577,6 +577,29 @@ def f(p): return t * d2 + (1 - t) * d1 return f +@op3 +def transition_spherical(f0, f1,r0=5, x0=0,y0=0,z0=0,k=1.0, e=ease.linear): + def f(p): + d1 = f0(p) + d2 = f1(p) + r = (p[:,0]-x0)**2+(p[:,1]-y0)**2+(p[:,2]-z0)**2 - r0**2 + t = 1./(1.0+np.exp(k*r)) + t = e(t).reshape((-1, 1)) + return t * d2 + (1 - t) * d1 + return f + +@op3 +def transition_sdf(f0, f1,f2,k=0.25,h=1.0, e=ease.linear): + def f(p): + d1 = f0(p) + d2 = f1(p) + d3 = f2(p) + r = _length(d3-h) + t = 1./(1.0+np.exp(k*r)) + t = e(t).reshape((-1, 1)) + return t * d2 + (1 - t) * d1 + return f + @op3 def transition_radial(f0, f1, r0=0, r1=1, e=ease.linear): def f(p): From 8f440b3bab23dfe363ac3b94c1d8596a75ffb05a Mon Sep 17 00:00:00 2001 From: mutablehurdle Date: Thu, 15 Sep 2022 08:17:46 -0700 Subject: [PATCH 7/9] Update d3.py --- sdf/d3.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sdf/d3.py b/sdf/d3.py index 3d6ff0f6..9aef1ea6 100644 --- a/sdf/d3.py +++ b/sdf/d3.py @@ -334,13 +334,13 @@ def f(p): @sdf3 -def gyroid(h,t,size,center = ORIGIN): +def gyroid(thickness,topology,n,size = (1,1,1),center = ORIGIN): size = np.array(size) def f(p): x = p[:,0] y = p[:,1] z = p[:,2] - d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-t)-h + d = np.abs(np.cos(n*x)*np.sin(n*y)+np.cos(n*y)*np.sin(n*z)+np.cos(n*z)*np.sin(n*x)-topology)-thickness q = np.abs(p - center) - size / 2 return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f From 96b8daeefa74874e768ba3df28cb49fd501f83a0 Mon Sep 17 00:00:00 2001 From: mutablehurdle Date: Thu, 22 Sep 2022 20:25:19 -0700 Subject: [PATCH 8/9] Update d3.py --- sdf/d3.py | 154 +++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 148 insertions(+), 6 deletions(-) diff --git a/sdf/d3.py b/sdf/d3.py index 9aef1ea6..7c1e4cb1 100644 --- a/sdf/d3.py +++ b/sdf/d3.py @@ -294,6 +294,23 @@ def f(p): return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f +@sdf3 +def cylindrical_MO(thickness,m,n,slant,mode = "vertical",size = 2*np.pi,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + rho = np.sqrt(x**2+y**2) + theta = np.arctan2(y,x) + if mode == "vertical": + d = np.abs(np.sin(z)+np.cos(m*theta+slant*np.sin(n*rho)))-thickness + elif mode == "horizontal": + d = np.abs(np.sin(z)+np.cos(m*rho+slant*np.sin(n*theta)))-thickness + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + @sdf3 def EB(h,size,center = ORIGIN): size = np.array(size) @@ -306,41 +323,166 @@ def f(p): return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f +@sdf3 +def cylindrical_EB(h,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + rho = np.sqrt(x**2+y**2) + theta = np.arctan2(x,y) + d = np.abs(np.cos(rho)+np.cos(theta)*np.cos(z))-h + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f #Minimal surfaces +# @sdf3 +# def schwarzP(h,size,center = ORIGIN): +# size = np.array(size) +# def f(p): +# x = p[:,0] +# y = p[:,1] +# z = p[:,2] +# d = np.abs(np.cos(x)+np.cos(y)+np.cos(z))-h +# q = np.abs(p - center) - size / 2 +# return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) +# return f + @sdf3 -def schwarzP(h,size,center = ORIGIN): +def schwarzP(thickness,topology,size,center = ORIGIN): size = np.array(size) def f(p): x = p[:,0] y = p[:,1] z = p[:,2] - d = np.abs(np.cos(x)+np.cos(y)+np.cos(z))-h + d = np.abs(np.cos(x)+np.cos(y)+np.cos(z)-topology)-thickness q = np.abs(p - center) - size / 2 return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f @sdf3 -def schwarzD(h,size, center = ORIGIN): +def cylindrical_schwarzP(thickness,topology,size,center = ORIGIN): size = np.array(size) def f(p): x = p[:,0] y = p[:,1] z = p[:,2] - d = np.abs(np.sin(x)*np.sin(y)*np.sin(z)+np.sin(x)*np.cos(y)*np.cos(z)+np.cos(x)*np.sin(y)*np.cos(z))-h + rho = np.sqrt(x**2+y**2) + theta = np.arctan2(x,y) + d = np.abs(np.cos(rho)+np.cos(theta)+np.cos(z)-topology)-thickness q = np.abs(p - center) - size / 2 return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f +# @sdf3 +# def schwarzD(h,size, center = ORIGIN): +# size = np.array(size) +# def f(p): +# x = p[:,0] +# y = p[:,1] +# z = p[:,2] +# d = np.abs(np.sin(x)*np.sin(y)*np.sin(z)+np.sin(x)*np.cos(y)*np.cos(z)+np.cos(x)*np.sin(y)*np.cos(z))-h +# q = np.abs(p - center) - size / 2 +# return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) +# return f +@sdf3 +def schwarzD(thickness,topology,size, center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.sin(x)*np.sin(y)*np.sin(z)+np.sin(x)*np.cos(y)*np.cos(z)+np.cos(x)*np.sin(y)*np.cos(z)-topology)-thickness + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +def cylindrical_schwarzD(thickness,topology,size, center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + rho = np.sqrt(x**2+y**2) + theta = np.arctan2(x,y) + d = np.abs(np.sin(rho)*np.sin(theta)*np.sin(z)+np.sin(rho)*np.cos(theta)*np.cos(z)+np.cos(rho)*np.sin(theta)*np.cos(z)-topology)-thickness + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +def fischer_koch(thickness,topology,size, center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.cos(2*x)*np.sin(y)*np.cos(z)+np.cos(2*y)*np.cos(x)*np.sin(z)+np.cos(y)*np.sin(x)*np.cos(2*z)-topology)-thickness + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +def lidinoid(thickness,topology,size, center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.sin(2*x)*np.cos(y)*np.sin(z)+np.sin(2*y)*np.cos(z)*np.sin(x)+np.sin(2*z)*np.cos(x)*np.sin(y)-np.cos(2*x)*np.cos(2*y)-np.cos(2*y)*np.cos(2*z)-np.cos(2*z)*np.cos(2*x)-topology)-thickness + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + +@sdf3 +def neovius(thickness,topology,size, center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(3*(np.cos(x)+np.cos(y)+np.cos(z))+4*np.cos(x)*np.cos(y)*np.cos(z)-topology)-thickness + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f +# @sdf3 +# def gyroid(thickness,topology,n,size = (2*np.pi,2*np.pi,2*np.pi),center = ORIGIN): +# size = np.array(size) +# def f(p): +# x = p[:,0] +# y = p[:,1] +# z = p[:,2] +# d = np.abs(np.cos(n*x)*np.sin(n*y)+np.cos(n*y)*np.sin(n*z)+np.cos(n*z)*np.sin(n*x)-topology)-thickness +# q = np.abs(p - center) - size / 2 +# return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) +# return f + +@sdf3 +def gyroid(thickness,topology,size,center = ORIGIN): + size = np.array(size) + def f(p): + x = p[:,0] + y = p[:,1] + z = p[:,2] + d = np.abs(np.cos(x)*np.sin(y)+np.cos(y)*np.sin(z)+np.cos(z)*np.sin(x)-topology)-thickness + q = np.abs(p - center) - size / 2 + return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) + return f + @sdf3 -def gyroid(thickness,topology,n,size = (1,1,1),center = ORIGIN): +def cylindrical_gyroid(thickness,topology,n,size = (2*np.pi,2*np.pi,2*np.pi),center = ORIGIN): size = np.array(size) def f(p): x = p[:,0] y = p[:,1] z = p[:,2] - d = np.abs(np.cos(n*x)*np.sin(n*y)+np.cos(n*y)*np.sin(n*z)+np.cos(n*z)*np.sin(n*x)-topology)-thickness + rho = np.sqrt(x**2+y**2) + theta = np.arctan2(x,y) + d = np.abs(np.cos(n*rho)*np.sin(n*theta)+np.cos(n*theta)*np.sin(n*z)+np.cos(n*z)*np.sin(n*rho)-topology)-thickness q = np.abs(p - center) - size / 2 return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0)) return f From a30359865da9054e0fdfc0107a4376a1eaf84709 Mon Sep 17 00:00:00 2001 From: mutablehurdle <80230725+mutablehurdle@users.noreply.github.com> Date: Wed, 5 Oct 2022 15:54:45 -0700 Subject: [PATCH 9/9] Update d3.py --- sdf/d3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sdf/d3.py b/sdf/d3.py index 7c1e4cb1..fdfd6a6b 100644 --- a/sdf/d3.py +++ b/sdf/d3.py @@ -481,7 +481,7 @@ def f(p): y = p[:,1] z = p[:,2] rho = np.sqrt(x**2+y**2) - theta = np.arctan2(x,y) + theta = np.arctan2(y,x) d = np.abs(np.cos(n*rho)*np.sin(n*theta)+np.cos(n*theta)*np.sin(n*z)+np.cos(n*z)*np.sin(n*rho)-topology)-thickness q = np.abs(p - center) - size / 2 return _max(d,_length(_max(q, 0)) + _min(np.amax(q, axis=1), 0))