- added sdf folder ( doesnt work via pip or git=)

.idea/fluency.iml

@@ -5,7 +5,7 @@
       <sourceFolder url="file://$MODULE_DIR$/sdfcad" isTestSource="false" />
       <excludeFolder url="file://$MODULE_DIR$/.venv" />
     </content>
-    <orderEntry type="inheritedJdk" />
+    <orderEntry type="jdk" jdkName="Python 3.12 (fluency)" jdkType="Python SDK" />
     <orderEntry type="sourceFolder" forTests="false" />
   </component>
 </module>

.idea/misc.xml

@@ -3,7 +3,7 @@
   <component name="Black">
     <option name="sdkName" value="Python 3.11 (fluency)" />
   </component>
-  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.11 (fluency)" project-jdk-type="Python SDK" />
+  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.12 (fluency)" project-jdk-type="Python SDK" />
   <component name="PyCharmProfessionalAdvertiser">
     <option name="shown" value="true" />
   </component>

.idea/vcs.xml

@@ -2,6 +2,5 @@
 <project version="4">
   <component name="VcsDirectoryMappings">
     <mapping directory="$PROJECT_DIR$" vcs="Git" />
-    <mapping directory="$PROJECT_DIR$/sdfcad" vcs="Git" />
   </component>
 </project>

.idea/workspace.xml

@@ -4,23 +4,14 @@
     <option name="autoReloadType" value="SELECTIVE" />
   </component>
   <component name="ChangeListManager">
-    <list default="true" id="8f0bafd6-58a0-4b20-aa2b-ddc3ba278873" name="Changes" comment="init">
+    <list default="true" id="8f0bafd6-58a0-4b20-aa2b-ddc3ba278873" name="Changes" comment="- added screenshot">
-      <change afterPath="$PROJECT_DIR$/.idea/fluency.iml" afterDir="false" />
+      <change beforePath="$PROJECT_DIR$/.idea/fluency.iml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/fluency.iml" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/.idea/inspectionProfiles/Project_Default.xml" afterDir="false" />
+      <change beforePath="$PROJECT_DIR$/.idea/misc.xml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/misc.xml" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/.idea/inspectionProfiles/profiles_settings.xml" afterDir="false" />
+      <change beforePath="$PROJECT_DIR$/.idea/vcs.xml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/vcs.xml" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/.idea/misc.xml" afterDir="false" />
+      <change beforePath="$PROJECT_DIR$/.idea/workspace.xml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/workspace.xml" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/.idea/modules.xml" afterDir="false" />
+      <change beforePath="$PROJECT_DIR$/main.py" beforeDir="false" afterPath="$PROJECT_DIR$/main.py" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/.idea/vcs.xml" afterDir="false" />
+      <change beforePath="$PROJECT_DIR$/mesh_modules/simple_mesh.py" beforeDir="false" afterPath="$PROJECT_DIR$/mesh_modules/simple_mesh.py" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/.idea/workspace.xml" afterDir="false" />
+      <change beforePath="$PROJECT_DIR$/requirements.txt" beforeDir="false" afterPath="$PROJECT_DIR$/requirements.txt" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/2dtest.py" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/fluencyb.py" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/main.py" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/meshtest.py" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/modules/gl_widget.py" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/modules/out.stl" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/side_fluency.py" afterDir="false" />
-      <change afterPath="$PROJECT_DIR$/vulkan.py" afterDir="false" />
-      <change beforePath="$PROJECT_DIR$/fluency.py" beforeDir="false" afterPath="$PROJECT_DIR$/fluency.py" afterDir="false" />
     </list>
     <option name="SHOW_DIALOG" value="false" />
     <option name="HIGHLIGHT_CONFLICTS" value="true" />
@@ -35,6 +26,11 @@
     </option>
   </component>
   <component name="Git.Settings">
+    <option name="RECENT_BRANCH_BY_REPOSITORY">
+      <map>
+        <entry key="$PROJECT_DIR$" value="structure" />
+      </map>
+    </option>
     <option name="RECENT_GIT_ROOT_PATH" value="$PROJECT_DIR$" />
     <option name="ROOT_SYNC" value="DONT_SYNC" />
   </component>
@@ -48,7 +44,6 @@
   &quot;associatedIndex&quot;: 6
 }</component>
   <component name="ProjectId" id="2aDywQvESFCKbJK4JUVHIhkN4S6" />
-  <component name="ProjectLevelVcsManager" settingsEditedManually="true" />
   <component name="ProjectViewState">
     <option name="hideEmptyMiddlePackages" value="true" />
     <option name="showLibraryContents" value="true" />
@@ -57,17 +52,23 @@
   "keyToString": {
     "Python.2dtest.executor": "Run",
     "Python.3d_windows.executor": "Run",
+    "Python.Unnamed.executor": "Run",
+    "Python.draw_widget2d.executor": "Run",
+    "Python.draw_widget_solve.executor": "Run",
     "Python.fluency.executor": "Run",
     "Python.fluencyb.executor": "Run",
     "Python.gl_widget.executor": "Run",
     "Python.main.executor": "Run",
     "Python.meshtest.executor": "Run",
     "Python.side_fluency.executor": "Run",
+    "Python.simple_mesh.executor": "Run",
+    "Python.vtk_widget.executor": "Run",
     "Python.vulkan.executor": "Run",
     "RunOnceActivity.OpenProjectViewOnStart": "true",
     "RunOnceActivity.ShowReadmeOnStart": "true",
+    "RunOnceActivity.git.unshallow": "true",
     "git-widget-placeholder": "master",
-    "last_opened_file_path": "/Volumes/Data_drive/Programming/fluency/modules",
+    "last_opened_file_path": "/Volumes/Data_drive/Programming/fluency",
     "settings.editor.selected.configurable": "project.propVCSSupport.DirectoryMappings"
   }
 }]]></component>
@@ -78,6 +79,8 @@
   </component>
   <component name="RecentsManager">
     <key name="CopyFile.RECENT_KEYS">
+      <recent name="$PROJECT_DIR$" />
+      <recent name="$PROJECT_DIR$/drawing_modules" />
       <recent name="$PROJECT_DIR$/modules" />
     </key>
     <key name="MoveFile.RECENT_KEYS">
@@ -87,7 +90,7 @@
   <component name="SharedIndexes">
     <attachedChunks>
       <set>
-        <option value="bundled-python-sdk-09665e90c3a7-d3b881c8e49f-com.jetbrains.pycharm.community.sharedIndexes.bundled-PC-233.15026.15" />
+        <option value="bundled-python-sdk-4c141bd692a7-e2d783800521-com.jetbrains.pycharm.community.sharedIndexes.bundled-PC-251.26927.90" />
       </set>
     </attachedChunks>
   </component>
@@ -108,7 +111,143 @@
       <option name="project" value="LOCAL" />
       <updated>1703951701948</updated>
     </task>
-    <option name="localTasksCounter" value="2" />
+    <task id="LOCAL-00002" summary="- Basic oop sketch widget implement">
+      <option name="closed" value="true" />
+      <created>1729958532384</created>
+      <option name="number" value="00002" />
+      <option name="presentableId" value="LOCAL-00002" />
+      <option name="project" value="LOCAL" />
+      <updated>1729958532384</updated>
+    </task>
+    <task id="LOCAL-00003" summary="- Sketch projection partly works again :)">
+      <option name="closed" value="true" />
+      <created>1735563255455</created>
+      <option name="number" value="00003" />
+      <option name="presentableId" value="LOCAL-00003" />
+      <option name="project" value="LOCAL" />
+      <updated>1735563255455</updated>
+    </task>
+    <task id="LOCAL-00004" summary="- Sketch projection partly works again :)">
+      <option name="closed" value="true" />
+      <created>1735585968733</created>
+      <option name="number" value="00004" />
+      <option name="presentableId" value="LOCAL-00004" />
+      <option name="project" value="LOCAL" />
+      <updated>1735585968733</updated>
+    </task>
+    <task id="LOCAL-00005" summary="- Added new componnt controls">
+      <option name="closed" value="true" />
+      <created>1735601610504</created>
+      <option name="number" value="00005" />
+      <option name="presentableId" value="LOCAL-00005" />
+      <option name="project" value="LOCAL" />
+      <updated>1735601610504</updated>
+    </task>
+    <task id="LOCAL-00006" summary="- Added new componnt controls">
+      <option name="closed" value="true" />
+      <created>1735601786207</created>
+      <option name="number" value="00006" />
+      <option name="presentableId" value="LOCAL-00006" />
+      <option name="project" value="LOCAL" />
+      <updated>1735601786207</updated>
+    </task>
+    <task id="LOCAL-00007" summary="- changing compos for sketches works">
+      <option name="closed" value="true" />
+      <created>1735652081552</created>
+      <option name="number" value="00007" />
+      <option name="presentableId" value="LOCAL-00007" />
+      <option name="project" value="LOCAL" />
+      <updated>1735652081552</updated>
+    </task>
+    <task id="LOCAL-00008" summary="- changing compos including sketches and bodies">
+      <option name="closed" value="true" />
+      <created>1735662119176</created>
+      <option name="number" value="00008" />
+      <option name="presentableId" value="LOCAL-00008" />
+      <option name="project" value="LOCAL" />
+      <updated>1735662119176</updated>
+    </task>
+    <task id="LOCAL-00009" summary="- Drawing bodys depending on the selected compo&#10;- Cut working&#10;- Edit sketch working">
+      <option name="closed" value="true" />
+      <created>1735685300102</created>
+      <option name="number" value="00009" />
+      <option name="presentableId" value="LOCAL-00009" />
+      <option name="project" value="LOCAL" />
+      <updated>1735685300102</updated>
+    </task>
+    <task id="LOCAL-00010" summary="- delete sketch working&#10;- added mid point snap&#10;- added hovering line with distance">
+      <option name="closed" value="true" />
+      <created>1735763743346</created>
+      <option name="number" value="00010" />
+      <option name="presentableId" value="LOCAL-00010" />
+      <option name="project" value="LOCAL" />
+      <updated>1735763743346</updated>
+    </task>
+    <task id="LOCAL-00011" summary="- Added new buttons and settings">
+      <option name="closed" value="true" />
+      <created>1735825176611</created>
+      <option name="number" value="00011" />
+      <option name="presentableId" value="LOCAL-00011" />
+      <option name="project" value="LOCAL" />
+      <updated>1735825176611</updated>
+    </task>
+    <task id="LOCAL-00012" summary="- Added new buttons and settings">
+      <option name="closed" value="true" />
+      <created>1735842523870</created>
+      <option name="number" value="00012" />
+      <option name="presentableId" value="LOCAL-00012" />
+      <option name="project" value="LOCAL" />
+      <updated>1735842523870</updated>
+    </task>
+    <task id="LOCAL-00013" summary="- Added construction lines switching&#10;- Moved callbacks into sketchwidget from main.&#10;- Changed reset on right click">
+      <option name="closed" value="true" />
+      <created>1739739664763</created>
+      <option name="number" value="00013" />
+      <option name="presentableId" value="LOCAL-00013" />
+      <option name="project" value="LOCAL" />
+      <updated>1739739664763</updated>
+    </task>
+    <task id="LOCAL-00014" summary="- Added contrain displayed next to line&#10;- Slight change to point check from solver.">
+      <option name="closed" value="true" />
+      <created>1743193041868</created>
+      <option name="number" value="00014" />
+      <option name="presentableId" value="LOCAL-00014" />
+      <option name="project" value="LOCAL" />
+      <updated>1743193041868</updated>
+    </task>
+    <task id="LOCAL-00015" summary="- Added enabling of midpsnap and prepared others&#10;- Show dimesnion on hover">
+      <option name="closed" value="true" />
+      <created>1743284173326</created>
+      <option name="number" value="00015" />
+      <option name="presentableId" value="LOCAL-00015" />
+      <option name="project" value="LOCAL" />
+      <updated>1743284173326</updated>
+    </task>
+    <task id="LOCAL-00016" summary="- Fixed redraw when component changed">
+      <option name="closed" value="true" />
+      <created>1744555255868</created>
+      <option name="number" value="00016" />
+      <option name="presentableId" value="LOCAL-00016" />
+      <option name="project" value="LOCAL" />
+      <updated>1744555255868</updated>
+    </task>
+    <task id="LOCAL-00017" summary="- added MIT license">
+      <option name="closed" value="true" />
+      <created>1748764814845</created>
+      <option name="number" value="00017" />
+      <option name="presentableId" value="LOCAL-00017" />
+      <option name="project" value="LOCAL" />
+      <updated>1748764814845</updated>
+    </task>
+    <task id="LOCAL-00018" summary="- added screenshot">
+      <option name="closed" value="true" />
+      <created>1748765318267</created>
+      <option name="number" value="00018" />
+      <option name="presentableId" value="LOCAL-00018" />
+      <option name="project" value="LOCAL" />
+      <updated>1748765318267</updated>
+    </task>
+    <option name="localTasksCounter" value="19" />
     <servers />
   </component>
   <component name="Vcs.Log.Tabs.Properties">
@@ -127,6 +266,21 @@
       <path value="$PROJECT_DIR$/pythonProject" />
     </ignored-roots>
     <MESSAGE value="init" />
-    <option name="LAST_COMMIT_MESSAGE" value="init" />
+    <MESSAGE value="- Basic oop sketch widget implement" />
+    <MESSAGE value="- Renabled extrusion with new object system" />
+    <MESSAGE value="- Sketch projection partly works again :)" />
+    <MESSAGE value="- Added new componnt controls" />
+    <MESSAGE value="- changing compos for sketches works" />
+    <MESSAGE value="- changing compos including sketches and bodies" />
+    <MESSAGE value="- Drawing bodys depending on the selected compo&#10;- Cut working&#10;- Edit sketch working" />
+    <MESSAGE value="- delete sketch working&#10;- added mid point snap&#10;- added hovering line with distance" />
+    <MESSAGE value="- Added new buttons and settings" />
+    <MESSAGE value="- Added construction lines switching&#10;- Moved callbacks into sketchwidget from main.&#10;- Changed reset on right click" />
+    <MESSAGE value="- Added contrain displayed next to line&#10;- Slight change to point check from solver." />
+    <MESSAGE value="- Added enabling of midpsnap and prepared others&#10;- Show dimesnion on hover" />
+    <MESSAGE value="- Fixed redraw when component changed" />
+    <MESSAGE value="- added MIT license" />
+    <MESSAGE value="- added screenshot" />
+    <option name="LAST_COMMIT_MESSAGE" value="- added screenshot" />
   </component>
 </project>

main.py

@@ -9,7 +9,8 @@ from PySide6.QtCore import Qt, QPoint, Signal, QSize
 from PySide6.QtWidgets import QApplication, QMainWindow, QSizePolicy, QInputDialog, QDialog, QVBoxLayout, QHBoxLayout, QLabel, QDoubleSpinBox, QCheckBox, QPushButton, QButtonGroup
 from Gui import Ui_fluencyCAD  # Import the generated GUI module
 from drawing_modules.vtk_widget import VTKWidget
-from drawing_modules.vysta_widget import PyVistaWidget
+import numpy as np
+
 from drawing_modules.draw_widget_solve import SketchWidget
 from sdf import *
 from python_solvespace import SolverSystem, ResultFlag

mesh_modules/simple_mesh.py

@@ -1,6 +1,6 @@
 import numpy as np
 from scipy.spatial import Delaunay, ConvexHull
-from shapely.geometry import Polygon, Point
+#from shapely.geometry import Polygon, Point
 
 
 def alpha_shape(points, alpha):

requirements.txt

@@ -1,60 +1,61 @@
-certifi==2025.1.31
+asttokens==3.0.0
-cffi==1.17.1
+attrs==25.3.0
-charset-normalizer==3.4.1
+black==24.10.0
-contourpy==1.3.1
+click==8.2.1
+contourpy==1.3.2
 cycler==0.12.1
+decorator==5.2.1
+executing==2.2.0
 flexcache==0.3
 flexparser==0.4
-fonttools==4.56.0
+fonttools==4.58.1
-freetype-py==2.5.1
+h5py==3.13.0
-hsluv==5.0.4
-idna==3.10
 imageio==2.37.0
+ipython==9.3.0
+ipython_pygments_lexers==1.1.1
+jedi==0.19.2
 kiwisolver==1.4.8
 lazy_loader==0.4
 markdown-it-py==3.0.0
-matplotlib==3.10.1
+matplotlib==3.10.3
+matplotlib-inline==0.1.7
 mdurl==0.1.2
 meshio==5.3.5
+mypy_extensions==1.1.0
 names==0.3.0
-networkx==3.4.2
+networkx==3.5
-Nuitka==2.6.9
+Nuitka==2.7.10
-numpy==2.2.4
+numpy==2.2.6
-numpy-stl==3.2.0
 ordered-set==4.1.0
-packaging==24.2
+packaging==25.0
-panda3d-gltf==1.3.0
+parso==0.8.4
-panda3d-simplepbr==0.13.0
+pathspec==0.12.1
-pillow==11.1.0
+pexpect==4.9.0
+pillow==11.2.1
 Pint==0.24.4
-platformdirs==4.3.7
+platformdirs==4.3.8
-pooch==1.8.2
+prompt_toolkit==3.0.51
-pycparser==2.22
+ptyprocess==0.7.0
-pygame==2.6.1
+pure_eval==0.2.3
 Pygments==2.19.1
 pyparsing==3.2.3
-PySide6_Essentials==6.8.3
+PySide6==6.9.0
+PySide6_Addons==6.9.0
+PySide6_Essentials==6.9.0
 python-dateutil==2.9.0.post0
-python-solvespace==3.0.8
+python_solvespace==3.0.8
-python-utils==3.9.1
-pyvista==0.44.2
-pyvistaqt==0.11.2
-QtPy==2.4.3
-requests==2.32.3
 rich==13.9.4
 scikit-image==0.25.2
-scipy==1.15.2
+scipy==1.15.3
-scooby==0.10.0
+sdfcad @ git+https://gitlab.com/nobodyinperson/sdfCAD@42505b5181c88dda2fd66ac9d387533fbe4145f3
-sdfcad @ git+https://gitlab.com/nobodyinperson/sdfCAD@9bd4e9021c6ee7e685ee28e8a3a5d2d2c028190c
+shiboken6==6.9.0
-shapely==2.1.0rc1
-shiboken6==6.8.3
 six==1.17.0
-tifffile==2025.3.13
+stack-data==0.6.3
-trimesh==4.6.5
+tifffile==2025.5.26
-tripy==1.0.0
+tokenize_rt==6.2.0
-typing_extensions==4.13.0
+traitlets==5.14.3
-urllib3==2.3.0
+typing_extensions==4.13.2
-vispy==0.14.3
+vtk==9.4.2
-vtk==9.4.1
+wcwidth==0.2.13
-vulkan==1.3.275.1
+xlrd==2.0.2
 zstandard==0.23.0

sdf/__init__.py

@@ -0,0 +1,26 @@
+from . import d2, d3, ease
+
+from .util import *
+from .units import units
+
+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,
+)

sdf/d2.py

@@ -0,0 +1,390 @@
+import functools
+import numpy as np
+import operator
+import copy
+
+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 fillet(self, r):
+        newSelf = copy.deepcopy(self)
+        newSelf._r = r
+        return newSelf
+
+    def radius(self, *args, **kwargs):
+        return self.fillet(*args, **kwargs)
+
+    def k(self, *args, **kwargs):
+        return self.fillet(*args, **kwargs)
+
+    def r(self, *args, **kwargs):
+        return self.fillet(*args, **kwargs)
+
+    def chamfer(self, c):
+        newSelf = copy.deepcopy(self)
+        newSelf._c = c
+        return newSelf
+
+    def c(self, *args, **kwargs):
+        return self.chamfer(*args, **kwargs)
+
+
+def sdf2(f):
+    @functools.wraps(f)
+    def wrapper(*args, **kwargs):
+        return SDF2(f(*args, **kwargs))
+
+    return wrapper
+
+
+def op2(f):
+    @functools.wraps(f)
+    def wrapper(*args, **kwargs):
+        return SDF2(f(*args, **kwargs))
+
+    _ops[f.__name__] = wrapper
+    return wrapper
+
+
+def op23(f):
+    @functools.wraps(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=None, diameter=None, center=ORIGIN):
+    if (radius is not None) == (diameter is not None):
+        raise ValueError(f"Specify either radius or diameter")
+    if radius is None:
+        radius = diameter / 2
+
+    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, r=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, r=r)
+
+
+@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(radius=None, diameter=None):
+    if (radius is not None) == (diameter is not None):
+        raise ValueError(f"Specify either radius or diameter")
+    if radius is None:
+        radius = diameter / 2
+    radius *= 3**0.5 / 2
+
+    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] * radius, k[2] * radius), np.zeros(len(p)) + radius
+        )
+        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
+
+
+def RegularPolygon(n, r=1):
+    ri = r * np.cos(np.pi / n)
+    return intersection(
+        *[slab(y0=-ri).rotate(a) for a in np.arange(0, 2 * np.pi, 2 * np.pi / n)]
+    )
+
+
+@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=np.inf):
+    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]]
+        # use horizontal distance to Z axis as X coordinate in 2D shape
+        # use Z coordinate as Y coordinate in 2D shape
+        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)
+mirror = op2(dn.mirror)
+modulate_between = op2(dn.modulate_between)
+stretch = op2(dn.stretch)

sdf/dn.py

@@ -0,0 +1,366 @@
+import itertools
+from functools import reduce, partial
+import warnings
+
+from . import ease
+
+import numpy as np
+
+_min = np.minimum
+_max = np.maximum
+
+
+def distance_to_plane(p, origin, normal):
+    """
+    Calculate the distance of a point ``p`` to the plane around ``origin`` with
+    normal ``normal``. This is dimension-independent, so e.g. the z-coordinate
+    can be omitted.
+
+    Args:
+        p (array): either [x,y,z] or [[x,y,z],[x,y,z],...]
+        origin (vector): a point on the plane
+        normal (vector): normal vector of the plane
+
+    Returns:
+        int: distance to plane
+    """
+    normal = normal / np.linalg.norm(normal)
+    return abs((p - origin) @ normal)
+
+
+def minimum(a, b, r=0):
+    if r:
+        Δ = b - a
+        h = np.clip(0.5 + 0.5 * Δ / r, 0, 1)
+        return b - Δ * h - r * h * (1 - h)
+    else:
+        return np.minimum(a, b)
+
+
+def maximum(a, b, r=0):
+    if r:
+        Δ = b - a
+        h = np.clip(0.5 - 0.5 * Δ / r, 0, 1)
+        return b - Δ * h + r * h * (1 - h)
+    else:
+        return np.maximum(a, b)
+
+
+def union(*sdfs, chamfer=0, c=0, radius=0, r=0, fillet=0, f=0):
+    c = max(chamfer, c)
+    r = max(radius, r, fillet, f)
+    sqrt05 = np.sqrt(0.5)
+
+    def f(p):
+        sdfs_ = iter(sdfs)
+        d1 = next(sdfs_)(p)
+        for sdf in sdfs_:
+            d2 = sdf(p)
+            R = r or getattr(sdf, "_r", 0)
+            C = c or getattr(sdf, "_c", 0)
+            parts = (d1, d2)
+            if C:
+                parts = (minimum(d1, d2), (d1 + d2 - C) * sqrt05)
+            d1 = minimum(*parts, R)
+
+        return d1
+
+    return f
+
+
+def intersection(*sdfs, chamfer=0, c=0, radius=0, r=0, fillet=0, f=0):
+    c = max(chamfer, c)
+    r = max(radius, r, fillet, f)
+    sqrt05 = np.sqrt(0.5)
+
+    def f(p):
+        sdfs_ = iter(sdfs)
+        d1 = next(sdfs_)(p)
+        for sdf in sdfs_:
+            d2 = sdf(p)
+            R = r or getattr(sdf, "_r", 0)
+            C = c or getattr(sdf, "_c", 0)
+            parts = (d1, d2)
+            if C:
+                parts = (maximum(d1, d2), (d1 + d2 + C) * sqrt05)
+            d1 = maximum(*parts, R)
+
+        return d1
+
+    return f
+
+
+def difference(*sdfs, chamfer=0, c=0, radius=0, r=0, fillet=0, f=0):
+    c = max(chamfer, c)
+    r = max(radius, r, fillet, f)
+    sqrt05 = np.sqrt(0.5)
+
+    def f(p):
+        sdfs_ = iter(sdfs)
+        d1 = next(sdfs_)(p)
+        for sdf in sdfs_:
+            d2 = sdf(p)
+            R = r or getattr(sdf, "_r", 0)
+            C = c or getattr(sdf, "_c", 0)
+            parts = (d1, -d2)
+            if C:
+                parts = (maximum(d1, -d2), (d1 - d2 + C) * sqrt05)
+            d1 = maximum(*parts, R)
+
+        return d1
+
+    return f
+
+
+def union_legacy(a, *bs, r=None):
+    def f(p):
+        d1 = a(p)
+        for b in bs:
+            d2 = b(p)
+            K = k or getattr(b, "_r", 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_legacy(a, *bs, r=None):
+    def f(p):
+        d1 = a(p)
+        for b in bs:
+            d2 = b(p)
+            K = k or getattr(b, "_r", 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_legacy(a, *bs, r=None):
+    def f(p):
+        d1 = a(p)
+        for b in bs:
+            d2 = b(p)
+            K = k or getattr(b, "_r", 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, r=0.5):
+    def f(p):
+        d1 = a(p)
+        for b in bs:
+            d2 = b(p)
+            K = k or getattr(b, "_r", 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=1, type="center"):
+    """
+    Keep only a margin of a given thickness around the object's boundary.
+
+    Args:
+        thickness (float): the resulting thickness
+        type (str): what kind of shell to generate.
+
+            ``"center"`` (default)
+                shell is spaced symmetrically around boundary
+            ``"outer"``
+                the resulting shell will be ``thickness`` larger than before
+            ``"inner"``
+                the resulting shell will be as large as before
+    """
+    return dict(
+        center=lambda p: np.abs(other(p)) - thickness / 2,
+        inner=other - other.erode(thickness),
+        outer=other.dilate(thickness) - other,
+    )[type]
+
+
+def modulate_between(sdf, a, b, e=ease.in_out_cubic):
+    """
+    Apply a distance offset transition between two control points
+    (e.g. make a rod thicker or thinner at some point or add a bump)
+
+    Args:
+        a, b (vectors): the two control points
+        e (scalar function): the distance offset function, will be called with
+            values between 0 (at control point ``a``) and 1 (at control point
+            ``b``).  Its result will be subtracted from the given SDF, thus
+            enlarging the object by that value.
+    """
+
+    # unit vector from control point a to b
+    ab = (ab := b - a) / (L := np.linalg.norm(ab))
+
+    def f(p):
+        # project current point onto control direction, clip and apply easing
+        offset = e(np.clip((p - a) @ ab / L, 0, 1))
+        return (dist := sdf(p)) - offset.reshape(dist.shape)
+
+    return f
+
+
+def stretch(sdf, a, b, symmetric=False, e=ease.linear):
+    """
+    Grab the object at point ``a`` and stretch the entire plane to ``b``.
+
+    Args:
+        a, b (point vectors): the control points
+        symmetric (bool): also stretch the same into the other direction.
+        e (Easing): easing to apply
+
+    Examples
+    ========
+
+    .. code-block:: python
+
+        # make a capsule
+        sphere(5).stretch(ORIGIN, 10*Z).save() # same as capsule(ORIGIN, 10*Z, 5)
+        # make an egg
+        sphere(5).stretch(ORIGIN, 10*Z, e=ease.smoothstep[:0.44]).save()
+    """
+    ab = (ab := b - a) / (L := np.linalg.norm(ab))
+
+    def f(p):
+        # s = ”how far are we between a and b as fraction?”
+        # if symmetric=True this also goes into the negative direction
+        s = np.clip((p - a) @ ab / L, -1 if symmetric else 0, 1)
+        # we return the sdf at a point 'behind' (p minus ...)
+        # the current point, but we go only as far back as the stretch distance
+        # at max
+        return sdf(p - (np.sign(s) * e(abs(s)) * L * ab[:, np.newaxis]).T)
+
+    return f
+
+
+def shear(sdf, fix, grab, move, e=ease.linear):
+    """
+    Grab the object at point ``grab`` and shear the entire plane in direction
+    ``move``, keeping point ``fix`` in place. If ``move`` is orthogonal to the
+    direction ``fix``->``grab``, then this operation is a shear.
+
+    Args:
+        fix, grab (point vectors): the control points
+        move (point vector): direction to shear to
+        e (Easing): easing to apply
+
+    Examples
+    ========
+
+    .. code-block:: python
+
+        # make a capsule
+        box([20,10,50]).shear(fix=-15*Z, grab=15*Z, move=-5*X, e=ease.smoothstep)
+    """
+    ab = (ab := grab - fix) / (L := np.linalg.norm(ab))
+
+    def f(p):
+        # s = ”how far are we between a and b as fraction?”
+        s = (p - fix) @ ab / L
+        return sdf(p - move * np.expand_dims(e(np.clip(s, 0, 1)), axis=1))
+
+    return f
+
+
+def mirror(other, direction, at=0):
+    """
+    Mirror around a given plane defined by ``origin`` reference point and
+    ``direction``.
+
+    Args:
+        direction (vector): direction to mirror to (e.g. :any:`X` to mirror along X axis)
+        at (3D vector): point to mirror at. Default is the origin.
+    """
+    direction = direction / np.linalg.norm(direction)
+
+    def f(p):
+        projdir = np.expand_dims((p - at) @ direction, axis=1) * direction
+        # mirrored point:
+        # - project 'p' onto 'direction' (result goes into 'projdir' direction)
+        # - projected point is at   'at + projdir'
+        # - remember direction from projected point to the original point (p - (at + projdir))
+        # - from origin 'at' go backwards the projected direction (at - projdir)
+        # - from that target, move along the remembered direction (p - (at + projdir))
+        # - pmirr = at - projdir + (p - (at + projdir))
+        # - the 'at' cancels out, the projdir is subtracted twice from the point
+        return other(p - 2 * projdir)
+
+    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

sdf/ease.py

@@ -0,0 +1,637 @@
+# system modules
+from dataclasses import dataclass
+from typing import Callable
+import itertools
+import functools
+import warnings
+
+# external modules
+import numpy as np
+import scipy.optimize
+
+
+@dataclass
+@functools.total_ordering
+class Extremum:
+    """
+    Container for min and max in Easing
+    """
+
+    pos: float
+    value: float
+
+    def __eq__(self, other):
+        return self.value == other.value
+
+    def __lt__(self, other):
+        return self.value < other.value
+
+
+@dataclass
+@functools.total_ordering
+class Easing:
+    """
+    A function defined on the interval [0;1]
+    """
+
+    f: Callable[float, float]
+    name: str
+
+    def modifier(decorated_fun):
+        @functools.wraps(decorated_fun)
+        def wrapper(self, *args, **kwargs):
+            newfun = decorated_fun(self, *args, **kwargs)
+            arglist = ",".join(
+                itertools.chain(map(str, args), (f"{k}={v}" for k, v in kwargs.items()))
+            )
+            newfun.__name__ = f"{self.f.__name__}.{decorated_fun.__name__}({arglist})"
+            return type(self)(f=newfun, name=newfun.__name__)
+
+        return wrapper
+
+    def __repr__(self):
+        return self.name
+
+    def __str__(self):
+        return self.name
+
+    @functools.cached_property
+    def is_ascending(self):
+        return np.all(np.diff(self.f(np.linspace(0, 1, 100))) >= 0)
+
+    @functools.cached_property
+    def is_symmetric(self):
+        t = np.linspace(0, 0.5, 100)
+        return np.allclose(self.f(t), self.f(1 - t))
+
+    @property
+    @modifier
+    def reverse(self):
+        """
+        Revert the function so it goes the other way round (starts at the end)
+        """
+        return lambda t: self.f(1 - t)
+
+    @property
+    @modifier
+    def symmetric(self):
+        """
+        Mirror and squash function to make it symmetric
+        """
+        return lambda t: self.f(-2 * (np.abs(t - 0.5) - 0.5))
+
+    @modifier
+    def mirror(self, x=None, y=None, copy=False):
+        """
+        Mirror function around an x and/or y value.
+
+        Args:
+            x (float): x value to mirror around
+            y (float): y value to mirror around
+            copy (bool): when mirroring around x, do copy-mirror
+        """
+        if (x, y) == (None, None):
+            x = 0.5
+
+        def mirrored(t):
+            if x is not None:
+                t = 2 * x - t
+            if copy:
+                t = np.abs(-t)
+            if y is None:
+                return self.f(t)
+            else:
+                return y - self.f(t)
+
+        return mirrored
+
+    @modifier
+    def clip(self, min=None, max=None):
+        """
+        Clip function at low and/or high values
+        """
+        if min is None and max is None:
+            min = 0
+            max = 1
+        return lambda t: np.clip(self.f(t), min, max)
+
+    @modifier
+    def clip_input(self, min=None, max=None):
+        """
+        Clip input parameter, i.e. extrapolate constantly outside the interval.
+        """
+        if min is None and max is None:
+            min = 0
+            max = 1
+        return lambda t: self.f(np.clip(t, min, max))
+
+    @property
+    @modifier
+    def clipped(self):
+        """
+        Clipped parameter and result to [0;1]
+        """
+        return lambda t: np.clip(self(np.clip(t, 0, 1)), 0, 1)
+
+    @modifier
+    def append(self, other, e=None):
+        """
+        Append another easing function and squish both into the [0;1] interval
+        """
+        if e is None:
+            e = in_out_square
+
+        def f(t):
+            mix = e(t)
+            return self.f(t * 2) * (1 - mix) + other((t - 0.5) * 2) * mix
+
+        return f
+
+    @modifier
+    def prepend(self, other, e=None):
+        """
+        Prepend another easing function and squish both into the [0;1] interval
+        """
+        if e is None:
+            e = in_out_square
+
+        def f(t):
+            mix = e(t)
+            return other(t * 2) * (1 - mix) + self.f((t - 0.5) * 2) * mix
+
+        return f
+
+    @modifier
+    def shift(self, offset):
+        """
+        Shift function on x-axis into positive direction by ``offset``.
+        """
+        return lambda t: self.f(t - offset)
+
+    @modifier
+    def repeat(self, n=2):
+        """
+        Repeat the function a total of n times in the interval [0;1].
+        """
+        return lambda t: self.f(t % (1 / n) * n)
+
+    @modifier
+    def multiply(self, factor):
+        """
+        Scale function by ``factor``
+        """
+        if isinstance(factor, Easing):
+            return lambda t: self(t) * factor(t)
+        else:
+            return lambda t: factor * self.f(t)
+
+    @modifier
+    def add(self, offset):
+        """
+        Add ``offset`` to function
+        """
+        if isinstance(offset, Easing):
+            return lambda t: self(t) + offset(t)
+        else:
+            return lambda t: self.f(t) + offset
+
+    def __add__(self, offset):
+        return self.add(offset)
+
+    def __radd__(self, offset):
+        return self.add(offset)
+
+    def __sub__(self, offset):
+        return self.add(-offset)
+
+    def __rsub__(self, offset):
+        return self.add(-offset)
+
+    def __mul__(self, factor):
+        return self.multiply(factor)
+
+    def __rmul__(self, factor):
+        return self.multiply(factor)
+
+    def __neg__(self):
+        return self.multiply(-1)
+
+    def __truediv__(self, factor):
+        return self.multiply(1 / factor)
+
+    def __or__(self, other):
+        return self.transition(other)
+
+    def __rshift__(self, offset):
+        return self.shift(offset)
+
+    def __lshift__(self, offset):
+        return self.shift(-offset)
+
+    def __getitem__(self, index):
+        if isinstance(index, Easing):
+            return self.chain(index)
+        if isinstance(index, slice):
+            return self.zoom(
+                0 if index.start is None else index.start,
+                1 if index.stop is None else index.stop,
+            )
+        else:
+            raise ValueError(
+                f"{index = } has to be slice of floats or an easing function"
+            )
+
+    @modifier
+    def chain(self, f=None):
+        """
+        Feed parameter through the given function before evaluating this function.
+        """
+        if f is None:
+            f = self.f
+        return lambda t: self.f(f(t))
+
+    @modifier
+    def zoom(self, left, right=None):
+        """
+        Arrange so that the interval [left;right] is moved into [0;1]
+        If only one argument is given, zoom in/out by moving edges that far.
+        """
+        if left is not None and right is None:
+            if left >= 0.5:
+                raise ValueError(
+                    f"{left = } is > 0.5 which doesn't make sense (bounds would cross)"
+                )
+            left = left
+            right = 1 - left
+        if left >= right:
+            raise ValueError(f"{right = } bound must be greater than {left = }")
+        return self.chain(linear.between(left, right)).f
+
+    @modifier
+    def between(self, left=0, right=1, e=None):
+        """
+        Arrange so ``f(0)==a`` and ``f(1)==b``.
+        """
+        f0, f1 = self.f(np.array([0, 1]))
+        la = f0 - left
+        lb = f1 - right
+        if e is None:  # linear is defined later
+            e = (
+                self  # use ourself as transition when we're ascending within [0;1]
+                if (self.is_ascending and np.allclose(self.f(np.array([0, 1])), [0, 1]))
+                else linear
+            )
+
+        def f(t):
+            t_ = e(t)
+            return self.f(t_) - (la * (1 - t_)) - lb * t_
+
+        return f
+
+    @modifier
+    def transition(self, other, e=None):
+        """
+        Transiton from one easing to another
+        """
+        if e is None:
+            e = linear
+
+        def f(t):
+            t_ = e(t)
+            return self.f(t) * (1 - t_) + other(t) * t_
+
+        return f
+
+    @classmethod
+    def function(cls, decorated_fun):
+        return cls(f=decorated_fun, name=decorated_fun.__name__)
+
+    def plot(self, *others, xlim=(0, 1), ax=None):
+        import matplotlib.pyplot as plt  # lazy import for speed
+        from cycler import cycler
+
+        if ax is None:
+            fig, ax_ = plt.subplots()
+        else:
+            ax_ = ax
+
+        try:
+            ax_.set_prop_cycle(
+                cycler(linestyle=["solid", "dashed", "dotted"], linewidth=[1, 1, 2])
+                * plt.rcParams["axes.prop_cycle"]
+            )
+        except ValueError as e:
+            pass
+
+        t = np.linspace(*xlim, 1000)
+        funs = list(others or [])
+        if isinstance(self, Easing):
+            funs.insert(0, self)
+        for f in funs:
+            ax_.plot(t, f(t), label=getattr(f, "name", getattr(f, "__name__", str(f))))
+        ax_.legend(ncol=int(np.ceil(len(ax_.get_lines()) / 10)))
+        if ax is None:
+            plt.show()
+        return ax_
+
+    @functools.cached_property
+    def min(self):
+        v = self.f(t := np.linspace(0, 1, 1000))
+        approxmin = Extremum(pos=t[i := np.argmin(v)], value=v[i])
+        opt = scipy.optimize.minimize(self, x0=[approxmin.pos], bounds=[(0, 1)])
+        optmin = Extremum(pos=opt.x[0], value=opt.fun)
+        return min(approxmin, optmin)
+
+    @functools.cached_property
+    def max(self):
+        """
+        Determine the maximum value
+        """
+        v = self.f(t := np.linspace(0, 1, 1000))
+        approxmax = Extremum(pos=t[i := np.argmax(v)], value=v[i])
+        opt = scipy.optimize.minimize(-self, x0=[approxmax.pos], bounds=[(0, 1)])
+        optmax = Extremum(pos=opt.x[0], value=-opt.fun)
+        return max(approxmax, optmax)
+
+    @functools.cached_property
+    def mean(self):
+        return np.mean(self.f(np.linspace(0, 1, 1000)))
+
+    def __lt__(self, e):
+        return np.all(self.f(t := np.linspace(0, 1, 50)) < e.f(t))
+
+    def __eq__(self, e):
+        return np.allclose(self.f(t := np.linspace(0, 1, 50)), e.f(t))
+
+    def __call__(self, t):
+        return self.f(t)
+
+
+@Easing.function
+def linear(t):
+    return t
+
+
+@Easing.function
+def in_quad(t):
+    return t * t
+
+
+@Easing.function
+def out_quad(t):
+    return -t * (t - 2)
+
+
+@Easing.function
+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)
+
+
+@Easing.function
+def in_cubic(t):
+    return t * t * t
+
+
+@Easing.function
+def out_cubic(t):
+    u = t - 1
+    return u * u * u + 1
+
+
+@Easing.function
+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)
+
+
+@Easing.function
+def in_quart(t):
+    return t * t * t * t
+
+
+@Easing.function
+def out_quart(t):
+    u = t - 1
+    return -(u * u * u * u - 1)
+
+
+@Easing.function
+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)
+
+
+@Easing.function
+def in_quint(t):
+    return t * t * t * t * t
+
+
+@Easing.function
+def out_quint(t):
+    u = t - 1
+    return u * u * u * u * u + 1
+
+
+@Easing.function
+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)
+
+
+@Easing.function
+def in_sine(t):
+    return -np.cos(t * np.pi / 2) + 1
+
+
+@Easing.function
+def out_sine(t):
+    return np.sin(t * np.pi / 2)
+
+
+@Easing.function
+def in_out_sine(t):
+    return -0.5 * (np.cos(np.pi * t) - 1)
+
+
+@Easing.function
+def in_expo(t):
+    a = np.zeros(len(t))
+    b = 2 ** (10 * (t - 1))
+    return np.where(t == 0, a, b)
+
+
+@Easing.function
+def out_expo(t):
+    a = np.zeros(len(t)) + 1
+    b = 1 - 2 ** (-10 * t)
+    return np.where(t == 1, a, b)
+
+
+@Easing.function
+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)))
+
+
+@Easing.function
+def in_circ(t):
+    return -1 * (np.sqrt(1 - t * t) - 1)
+
+
+@Easing.function
+def out_circ(t):
+    u = t - 1
+    return np.sqrt(1 - u * u)
+
+
+@Easing.function
+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)
+
+
+@Easing.function
+def in_elastic(t, k=0.5):
+    u = t - 1
+    return -1 * (2 ** (10.0 * u) * np.sin((u - k / 4) * (2 * np.pi) / k))
+
+
+@Easing.function
+def out_elastic(t, k=0.5):
+    return 2 ** (-10.0 * t) * np.sin((t - k / 4) * (2 * np.pi / k)) + 1
+
+
+@Easing.function
+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)
+
+
+@Easing.function
+def in_back(t):
+    k = 1.70158
+    return t * t * ((k + 1) * t - k)
+
+
+@Easing.function
+def out_back(t):
+    k = 1.70158
+    u = t - 1
+    return u * u * ((k + 1) * u + k) + 1
+
+
+@Easing.function
+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)
+
+
+@Easing.function
+def in_bounce(t):
+    return 1 - out_bounce(1 - t)
+
+
+@Easing.function
+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)))
+
+
+@Easing.function
+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)
+
+
+@Easing.function
+def in_square(t):
+    return np.heaviside(t - 1, 0)
+
+
+@Easing.function
+def out_square(t):
+    return np.heaviside(t + 1, 0)
+
+
+@Easing.function
+def in_out_square(t):
+    return np.heaviside(t - 0.5, 0)
+
+
+def constant(x):
+    return Easing(f=lambda t: np.full_like(t, x), name=f"constant({x})")
+
+
+zero = constant(0)
+one = constant(1)
+
+
+@Easing.function
+def smoothstep(t):
+    t = np.clip(t, 0, 1)
+    return 3 * t * t - 2 * t * t * t
+
+
+def _main():
+    import matplotlib.pyplot as plt
+    from cycler import cycler
+
+    plt.rcParams["axes.prop_cycle"] *= cycler(
+        linestyle=["solid", "dashed", "dotted"], linewidth=[1, 2, 3]
+    )
+    plt.rcParams["figure.autolayout"] = True
+    plt.rcParams["axes.grid"] = True
+    plt.rcParams["axes.axisbelow"] = True
+    plt.rcParams["legend.fontsize"] = "small"
+    LOCALS = globals()
+    print(f"{LOCALS = }")
+    fig, axes = plt.subplots(nrows=2)
+    Easing.plot(
+        *sorted((obj for n, obj in LOCALS.items() if isinstance(obj, Easing)), key=str),
+        ax=axes[0],
+    )
+    Easing.plot(
+        in_sine.symmetric,
+        in_out_sine.symmetric.multiply(-0.6),
+        linear.symmetric.multiply(-0.7),
+        in_out_sine.multiply(-0.6).symmetric,
+        out_sine.multiply(-0.6).reverse.symmetric.multiply(2),
+        out_bounce.add(-0.5),
+        ax=axes[1],
+    )
+    axes[0].set_title("Standard")
+    axes[1].set_title("Derived")
+    plt.show()
+
+
+if __name__ == "__main__":
+    _main()

sdf/errors.py

@@ -0,0 +1,42 @@
+import warnings
+import functools
+
+
+class SDFCADError(Exception):
+    pass
+
+
+class SDFCADInfiniteObjectError(Exception):
+    """
+    Error raised when an infinite object is encountered where not suitable.
+    """
+
+    pass
+
+
+class SDFCADWarning(Warning):
+    pass
+
+
+class SDFCADAlphaQualityWarning(SDFCADWarning):
+    show = True
+
+
+def alpha_quality(decorated_fun):
+    @functools.wraps(decorated_fun)
+    def wrapper(*args, **kwargs):
+        if SDFCADAlphaQualityWarning.show:
+            warnings.warn(
+                f"{decorated_fun.__name__}() is alpha quality "
+                f"and might give wrong results. Use with care. "
+                f"Hide this warning by setting sdf.errors.SDFCADAlphaQualityWarning.show=False.",
+                SDFCADAlphaQualityWarning,
+            )
+            with warnings.catch_warnings():
+                # Don't reissue nested alpha quality warnings
+                warnings.simplefilter("ignore", SDFCADAlphaQualityWarning)
+                return decorated_fun(*args, **kwargs)
+        else:
+            return decorated_fun(*args, **kwargs)
+
+    return wrapper

sdf/mesh.py

@@ -0,0 +1,282 @@
+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**18
+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)
+    f = partial(_worker, sdf, sparse=sparse)
+    with ThreadPool(workers) as pool:
+        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 str(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()

sdf/progress.py

@@ -0,0 +1,83 @@
+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)

sdf/stl.py

@@ -0,0 +1,27 @@
+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", ("<f", 3)),
+            ("points", ("<f", (3, 3))),
+            ("attr", "<H"),
+        ]
+    )
+
+    a = np.zeros(n, dtype=dtype)
+    a["points"] = points
+    a["normal"] = normals
+
+    with open(path, "wb") as fp:
+        fp.write(b"\x00" * 80)
+        fp.write(struct.pack("<I", n))
+        fp.write(a.tobytes())

sdf/text.py

@@ -0,0 +1,160 @@
+from PIL import Image, ImageFont, ImageDraw
+import scipy.ndimage as nd
+import numpy as np
+
+from . import d2
+
+# TODO: add support for newlines?
+
+PIXELS = 2**22
+
+
+def _load_image(thing):
+    if isinstance(thing, str):
+        return Image.open(thing)
+    elif isinstance(thing, (np.ndarray, np.generic)):
+        return Image.fromarray(thing)
+    return Image.fromarray(np.array(thing))
+
+
+def measure_text(name, text, width=None, height=None):
+    font = ImageFont.truetype(name, 96)
+    x0, y0, x1, y1 = font.getbbox(text)
+    aspect = (x1 - x0) / (y1 - y0)
+    if width is None and height is None:
+        height = 1
+    if width is None:
+        width = height * aspect
+    if height is None:
+        height = width / aspect
+    return (width, height)
+
+
+def measure_image(thing, width=None, height=None):
+    im = _load_image(thing)
+    w, h = im.size
+    aspect = w / h
+    if width is None and height is None:
+        height = 1
+    if width is None:
+        width = height * aspect
+    if height is None:
+        height = width / aspect
+    return (width, height)
+
+
+@d2.sdf2
+def text(font_name, text, width=None, height=None, pixels=PIXELS, points=512):
+    # load font file
+    font = ImageFont.truetype(font_name, points)
+
+    # compute texture bounds
+    p = 0.2
+    x0, y0, x1, y1 = font.getbbox(text)
+    px = int((x1 - x0) * p)
+    py = int((y1 - y0) * p)
+    tw = x1 - x0 + 1 + px * 2
+    th = y1 - y0 + 1 + py * 2
+
+    # render text to image
+    im = Image.new("L", (tw, th))
+    draw = ImageDraw.Draw(im)
+    draw.text((px - x0, py - y0), text, font=font, fill=255)
+
+    return _sdf(width, height, pixels, px, py, im)
+
+
+@d2.sdf2
+def image(thing, width=None, height=None, pixels=PIXELS):
+    im = _load_image(thing).convert("L")
+    return _sdf(width, height, pixels, 0, 0, im)
+
+
+def _sdf(width, height, pixels, px, py, im):
+    tw, th = im.size
+
+    # downscale image if necessary
+    factor = (pixels / (tw * th)) ** 0.5
+    if factor < 1:
+        tw, th = int(round(tw * factor)), int(round(th * factor))
+        px, py = int(round(px * factor)), int(round(py * factor))
+        im = im.resize((tw, th))
+
+    # convert to numpy array and apply distance transform
+    im = im.convert("1")
+    a = np.array(im)
+    inside = -nd.distance_transform_edt(a)
+    outside = nd.distance_transform_edt(~a)
+    texture = np.zeros(a.shape)
+    texture[a] = inside[a]
+    texture[~a] = outside[~a]
+
+    # save debug image
+    # a = np.abs(texture)
+    # lo, hi = a.min(), a.max()
+    # a = (a - lo) / (hi - lo) * 255
+    # im = Image.fromarray(a.astype('uint8'))
+    # im.save('debug.png')
+
+    # compute world bounds
+    pw = tw - px * 2
+    ph = th - py * 2
+    aspect = pw / ph
+    if width is None and height is None:
+        height = 1
+    if width is None:
+        width = height * aspect
+    if height is None:
+        height = width / aspect
+    x0 = -width / 2
+    y0 = -height / 2
+    x1 = width / 2
+    y1 = height / 2
+
+    # scale texture distances
+    scale = width / tw
+    texture *= scale
+
+    # prepare fallback rectangle
+    # TODO: reduce size based on mesh resolution instead of dividing by 2
+    rectangle = d2.rectangle((width / 2, height / 2))
+
+    def f(p):
+        x = p[:, 0]
+        y = p[:, 1]
+        u = (x - x0) / (x1 - x0)
+        v = (y - y0) / (y1 - y0)
+        v = 1 - v
+        i = u * pw + px
+        j = v * ph + py
+        d = _bilinear_interpolate(texture, i, j)
+        q = rectangle(p).reshape(-1)
+        outside = (i < 0) | (i >= 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

sdf/units.py

@@ -0,0 +1,3 @@
+import pint
+
+units = pint.UnitRegistry()

sdf/util.py

@@ -0,0 +1,32 @@
+import math
+import functools
+import inspect
+import numpy as np
+
+pi = math.pi
+
+degrees = math.degrees
+radians = math.radians
+
+
+def n_trailing_ascending_positive(d):
+    """
+    Determine how many elements in a given sequence are positive and ascending.
+
+    Args:
+        d (sequence of numbers): the sequence to check
+
+    Returns:
+        int : the amount of trailing ascending positive elements
+    """
+    d = np.array(d).flatten()
+    # is the next element larger than previous and positive?
+    order = (d[1:] > d[:-1]) & (d[:-1] > 0)
+    # TODO: Not happy at all with this if/else mess. Is there no easier way to find the
+    #       index in a numpy array after which the values are only ascending? 🤔
+    if np.all(order):  # all ascending
+        return d.size
+    elif np.all(~order):  # none ascending
+        return 0
+    else:  # count from end how many are ascending
+        return np.argmin(order[::-1]) + 1