Add note on parallel_transport

.gitignore

@@ -1,4 +1,5 @@
 *~
 #*#
-__pycache__/*
-.ipynb_checkpoints/*
+__pycache__/
+.ipynb_checkpoints/
+.ccls-cache/

README.md

@@ -1,84 +1,52 @@
-# To-do items, wanted features, bugs:
+# automata_scratch
 
-## Cool
-- More complicated: Examples of *merging*. I'm not sure on the theory
-  behind this.
-  
-## Annoying/boring
-- https://en.wikipedia.org/wiki/Polygon_triangulation - do this to
-  fix my wave example!
-  - http://www.polygontriangulation.com/2018/07/triangulation-algorithm.html
-- Clean up examples.ram_horn_branch(). The way I clean it up might
-  help inform some cleaner designs.
-- I really need to standardize some of the behavior of fundamental
-  operations (with regard to things like sizes they generate). This
-  is behavior that, if it changes, will change a lot of things that I'm
-  trying to keep consistent so that my examples still work.
-- Winding order.  It is consistent through seemingly
-  everything, except for reflection and close_boundary_simple.
-  (When there are two parallel boundaries joined with something like
-  join_boundary_simple, traversing these boundaries in their actual order
-  to generate triangles - like in close_boundary_simple - will produce
-  opposite winding order on each. Imagine a transparent clock: seen from the
-  front, it moves clockwise, but seen from the back, it moves
-  counter-clockwise.)
-- File that bug that I've seen in trimesh/three.js
-  (see trimesh_fail.ipynb)
-- Why do I get the weird zig-zag pattern on the triangles,
-  despite larger numbers of them? Is it something in how I
-  twist the frames?
-  - How can I compute the *torsion* on a quad? I think it
-    comes down to this: torsion applied across the quad I'm
-    triangulating leading to neither diagonal being a
-    particularly good choice.  Subdividing the boundary seems
-    to help, but other triangulation methods (e.g. turning a
-    quad to 4 triangles by adding the centroid) could be good
-    too.
-  - Facets/edges are just oriented the wrong way...
-  - Picking at random the diagonal on the quad to triangulate with
-    does seem to turn 'error' just to noise, and in its own way this
-    is preferable.
-- Integrate parallel_transport work and reuse what I can
-- /mnt/dev/graphics_misc/isosurfaces_2018_2019 - perhaps include my
-  spiral isosurface stuff from here
+This is repo has a few projects that are related in terms of
+high-level goal, but almost completely unrelated in their descent.
 
-## Abstractions
+- `python_isosurfaces_2018_2019` is some Python & Maxima code from
+  2018-2019 from me trying to turn the usual spiral isosurface into a
+  parametric formula of sorts in order to triangulate it more
+  effectively.
+- `parallel_transport` is some Python code from 2019 September which
+  implemented parallel frame transport, i.e.
+  [Parallel Transport Approach to Curve Framing](https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.8103).  It is mostly scratch
+  code / proof-of-concept.
+- `python_extrude_meshgen` is some Python code from around 2019
+  September which did a sort of extrusion-based code generation.
+  While this had some good results and some good ideas, the basic
+  model was too limited in terms of the topology it could express.
+- `libfive_subdiv` is a short project around 2021 July attempting to
+  use the Python bindings of [libfive](https://www.libfive.com/), and
+  automatic differentiation in
+  [autograd](https://github.com/HIPS/autograd), to turn implicit
+  surfaces to meshes which were suitable for subdivision via something
+  like
+  [OpenSubdiv](https://graphics.pixar.com/opensubdiv/overview.html)
+  (in turn so that I could render with them without having to use
+  insane numbers of triangles or somehow hide the obvious errors in
+  the geometry).  Briefly, the process was to use edges with crease
+  weights which were set based on the curvature of the implicit
+  surface.  While I accomplished this process, it didn't fulfill the
+  goal.  Shortly thereafter, I was re-reading
+  [Massively Parallel Rendering of Complex Closed-Form Implicit Surfaces](https://www.mattkeeter.com/research/mpr/) - which, like libfive, is by Matt Keeter -
+  and found a section I'd ignored on the difficulties of producing
+  good meshes from isosurfaces for the sake of rendering.  I kept
+  the code around because I figured it would be useful to refer to
+  later, particularly for the integration with Blender - but
+  otherwise shelved this effort.
+- `blender_scraps` contains some scraps of Python code meant to be
+  used inside of Blender's Python scripting - and it contains some
+  conversions from another project, Prosha, for procedural mesh
+  generation in Rust (itself based on learnings from
+  `python_extrude_meshgen`).  These examples were proof-of-concept of
+  generating meshes as control cages rather than as "final" meshes.
 
-- This has a lot of functions parametrized over a lot
-  of functions.  Need to work with this somehow.  (e.g. should
-  it subdivide this boundary? should it merge opening/closing
-  boundaries?)
-- Some generators produce boundaries that can be directly merged
-  and produce sensible geometry.  Some generators produce
-  boundaries that are only usable when they are further
-  transformed (and would produce degenerate geometry).  What sort
-  of nomenclature captures this?
+It would probably make sense to rename this repo to something with
+`procedural` in the name rather than `automata` since at some point it
+ceased to have much to do with automata.
 
-- How can I capture the idea of a group of parameters which, if
-  they are all scaled in the correct way (some linearly, others
-  inversely perhaps), generated geometry that is more or less
-  identical except that it is higher-resolution?
-- Use mixins to extend 3D transformations to things (matrices,
-  cages, meshes, existing transformations)
-- I can transform a Cage.  Why not a CageGen?
-  
-## ????
-- Embed this in Blender?
-  
-## Future thoughts
+## Projects not covered here
 
-- What if I had a function that could generate a Cage as if
-  from a parametric formula and smoothly vary its orientation?
-  My existing tools could easily turn this to a mesh. If I could vary
-  the detail of the Cage itself (if needed), then I could also
-  generate a mesh at an arbitrary level of detail simply by sampling at
-  finer and finer points on the parameter space.  (This might also tie
-  into the Parallel Transport work.)
-- What are the limitations of using Cages?
-- Current system is very "generative".  Could I do basically L-system
-  if I have rules for how a much is *refined*?  What about IFS?
-- Do this in Rust once I understand WTF I am doing
-
-## Other thoughts
-
-- Why do I never use the term "extruding" to describe what I'm doing?
+- curl-noise work (both in Clojure and in Python/vispy)
+- parallel transport
+- prosha, of course

blender_scraps/barbs.py

@@ -47,6 +47,8 @@ class Barbs(object):
         # be converted last-minute to tuples. (Why: we need to refer
         # to prior vertices and arithmetic is much easier from an
         # array, but Blender eventually needs tuples.)
+        self.creases_side = set()
+        self.creases_joint = set()
 
     def run(self, iters) -> (list, list):
         # 'iters' is ignored for now
@@ -87,6 +89,10 @@ class Barbs(object):
                   [bound[2], bound[3], a0 + 3, a0 + 2])
         self.barb(iters - 1, xform.compose(self.sides[3]),
                   [bound[3], bound[0], a0 + 0, a0 + 3])
+        for i in range(4):
+            j = (i + 1) % 4
+            self.creases_joint.add((a0 + i, a0 + j))
+            self.creases_side.add((bound[i], a0 + i))
 
     def barb(self, iters, xform, bound):
         if self.limit_check(xform, iters):

blender_scraps/scratch.py

@@ -0,0 +1,79 @@
+# This is a pile of patterns and snippets I've used in Blender in the
+# past; it isn't meant to be run on its own.
+
+import bmesh
+import bpy
+
+# Here is a good starting place (re: creating geometry):
+# https://docs.blender.org/api/current/info_gotcha.html#n-gons-and-tessellation
+
+# https://wiki.blender.org/wiki/Source/Modeling/BMesh/Design
+
+# Current object's data (must be selected):
+me = bpy.context.object.data
+
+def bmesh_set_creases(obj, vert_pairs, crease_val):
+    # Walk through the edges in 'obj'. For those *undirected* edges in
+    # 'vert_pairs' (a set of (vi, vj) tuples, where vi and vj are vertex
+    # indices, and tuple order is irrelevant), set the crease to 'crease_val'.
+    bm = bmesh.new()
+    bm.from_mesh(obj)
+    creaseLayer = bm.edges.layers.crease.verify()
+    for e in bm.edges:
+        idxs = tuple([v.index for v in e.verts])
+        print(idxs)
+        if idxs in vert_pairs or idxs[::-1] in vert_pairs:
+            e[creaseLayer] = crease_val
+    bm.to_mesh(obj)
+    bm.free()
+
+# My bpy.types.MeshPolygon objects:
+for i,poly in enumerate(me.polygons):
+    t = type(poly)
+    #print(f"poly {i}: {t}")
+    verts = list(poly.vertices)
+    print(f"poly {poly.index}: vertices={verts}")
+    #s = poly.loop_start
+    #n = poly.loop_total
+    #print(f"  loop_start={s} loop_total={n}")
+    #v = [l.vertex_index for l in me.loops[s:(s+n)]]
+    #print(f"  loop: {v}")
+    # Vector type:
+    v2 = [me.vertices[i].co for i in verts]
+    print(f"  verts: {v2}")
+    # Yes, this works:
+    #for i in verts:
+    #    me.vertices[i].co.x -= 1.0
+
+# Pattern for loading external module from Blender's Python (and
+# reloading as necessary):
+import sys
+ext_path = "/home/hodapp/source/automata_scratch/blender_scraps"
+if ext_path not in sys.path:
+    sys.path.append(ext_path)
+import whatever
+whatever = importlib.reload(whatever)
+# Note that if 'whatever' itself imports modules that may have changed
+# since the last import, you may need to do this same importlib
+# incantation!
+
+# Crease access - but the wrong way to change them:
+for edge in me.edges:
+    v = list(edge.vertices)
+    print(f"edge {edge.index}: crease={edge.crease} vertices={v}")
+    #edge.crease = 0.7
+
+# Creating a mesh with vertices & faces in Python via bpy with:
+# v - list of (x, y, z) tuples
+# f - list of (v0, v1, v2...) tuples, each with face's vertex indices
+mesh = bpy.data.meshes.new('mesh_thing')
+mesh.from_pydata(v, [], f)
+mesh.update(calc_edges=True)
+# set creases beforehand:
+#bmesh_set_creases(mesh, b.creases_joint, 0.7)
+obj = bpy.data.objects.new('obj_thing', mesh)
+# set obj's transform matrix:
+#obj.matrix_world = Matrix(...)
+# also acceptable to set creases:
+#bmesh_set_creases(obj.data, b.creases_joint, 0.7)
+bpy.context.scene.collection.objects.link(obj)

blender_scraps/tree_thing.py

@@ -0,0 +1,163 @@
+# Hasty conversion from the Rust in prosha/src/examples.rs & Barbs
+
+# This is mostly right, except:
+# - It doesn't yet do creases.
+
+import numpy as np
+import xform
+
+# Mnemonics:
+X = np.array([1.0, 0.0, 0.0])
+Y = np.array([0.0, 1.0, 0.0])
+Z = np.array([0.0, 0.0, 1.0])
+
+class TreeThing(object):
+    def __init__(self, f: float=0.6, depth: int=10, scale_min: float=0.02):
+        self.scale_min = scale_min
+        v = np.array([-1.0, 0.0, 1.0])
+        v /= np.linalg.norm(v)
+        self.incr = (xform.Transform().
+                     translate(0, 0, 0.9*f).
+                     rotate(v, 0.4*f).
+                     scale(1.0 - (1.0 - 0.95)*f))
+        # 'Base' vertices, used throughout:
+        self.base = np.array([
+            [-0.5, -0.5, 0.0],
+            [-0.5,  0.5, 0.0],
+            [ 0.5,  0.5, 0.0],
+            [ 0.5, -0.5, 0.0],
+        ])
+        # 'Transition' vertices:
+        self.trans = np.array([
+            # Top edge midpoints:
+            [-0.5,  0.0, 0.0],  # 0 - connects b0-b1
+            [ 0.0,  0.5, 0.0],  # 2 - connects b2-b3
+            [ 0.5,  0.0, 0.0],  # 1 - connects b1-b2
+            [ 0.0, -0.5, 0.0],  # 3 - connects b3-b0
+            # Top middle:
+            [ 0.0,  0.0, 0.0],  # 4 - midpoint/centroid of all
+        ])
+        # splits[i] gives transformation from a 'base' layer to the
+        # i'th split (0 to 3):
+        self.splits = [
+            xform.Transform().
+            rotate(Z, np.pi/2 * i).
+            translate(0.25, 0.25, 0.0).
+            scale(0.5)
+            for i in range(4)
+        ]
+        # Face & vertex accumulators:
+        self.faces = []
+        # self.faces will be a list of tuples (each one of length 4
+        # and containing indices into self.verts)
+        self.verts = []
+        # self.verts will be a list of np.array of shape (3,) but will
+        # be converted last-minute to tuples. (Why: we need to refer
+        # to prior vertices and arithmetic is much easier from an
+        # array, but Blender eventually needs tuples.)
+        self.creases_side = set()
+        self.creases_joint = set()
+        self.depth = depth
+
+    def run(self):
+        self.verts.extend(self.base)
+        self.faces.append((0, 1, 2, 3))
+        self.child(xform.Transform(), self.depth, [0, 1, 2, 3])
+        verts = [tuple(v) for v in self.verts]
+        faces = [tuple(f) for f in self.faces]
+        return verts, faces
+
+    def trunk(self, xf: xform.Transform, b):
+
+        if self.limit_check(xf):
+            # Note opposite winding order
+            verts = [b[i] for i in [3,2,1,0]]
+            self.faces.append(verts)
+            return
+
+        incr = (xform.Transform().
+                translate(0.0, 0.0, 1.0).
+                rotate(Z, 0.15).
+                rotate(X, 0.1).
+                scale(0.95))
+        sides = [
+            xform.Transform().
+            rotate(Z, -np.pi/2 * i).
+            rotate(Y, -np.pi/2).
+            translate(0.5, 0.0, 0.5)
+            for i in range(4)
+        ]
+        xf2 = xf.compose(incr)
+        g = xf2.apply_to(self.base)
+        a0 = len(self.verts)
+        self.verts.extend(g)
+
+        # TODO: Turn this to a cleaner loop?
+        self.main(iters - 1, xf2, [a0, a0 + 1, a0 + 2, a0 + 3])
+        self.child(iters - 1, xf.compose(self.sides[0]),
+                   [b[0], b[1], a0 + 1, a0 + 0])
+        self.child(iters - 1, xf.compose(self.sides[1]),
+                   [b[1], b[2], a0 + 2, a0 + 1])
+        self.child(iters - 1, xf.compose(self.sides[2]),
+                   [b[2], b[3], a0 + 3, a0 + 2])
+        self.child(iters - 1, xf.compose(self.sides[3]),
+                   [b[3], b[0], a0 + 0, a0 + 3])
+
+    def limit_check(self, xf: xform.Transform) -> bool:
+        # Assume all scales are the same (for now)
+        sx,_,_ = xf.get_scale()
+        return sx < self.scale_min
+
+    def child(self, xf: xform.Transform, depth, b):
+        if self.limit_check(xf):
+            # Note opposite winding order
+            verts = [b[i] for i in [3,2,1,0]]
+            self.faces.append(verts)
+            return
+
+        xf2 = xf.compose(self.incr)
+        if depth > 0:
+            # Just recurse on the current path:
+            n0 = len(self.verts)
+            self.verts.extend(xf2.apply_to(self.base))
+
+            # Connect parallel faces:
+            n = len(self.base)
+            for i, b0 in enumerate(b):
+                j = (i + 1) % n
+                b1 = b[j]
+                a0 = n0 + i
+                a1 = n0 + j
+                self.faces.append((a0, a1, b1, b0))
+
+            self.child(xf2, depth - 1, [n0, n0 + 1, n0 + 2, n0 + 3]);
+        else:
+            n = len(self.verts)
+            self.verts.extend(xf2.apply_to(self.base))
+            m01 = len(self.verts)
+            self.verts.extend(xf2.apply_to(self.trans))
+            m12, m23, m30, c = m01 + 1, m01 + 2, m01 + 3, m01 + 4
+            self.faces.extend([
+                # two faces straddling edge from vertex 0:
+                (b[0], n+0, m01),
+                (b[0], m30, n+0),
+                # two faces straddling edge from vertex 1:
+                (b[1], n+1, m12),
+                (b[1], m01, n+1),
+                # two faces straddling edge from vertex 2:
+                (b[2], n+2, m23),
+                (b[2], m12, n+2),
+                # two faces straddling edge from vertex 3:
+                (b[3], n+3, m30),
+                (b[3], m23, n+3),
+                # four faces from edge (0,1), (1,2), (2,3), (3,0):
+                (b[0], m01, b[1]),
+                (b[1], m12, b[2]),
+                (b[2], m23, b[3]),
+                (b[3], m30, b[0]),
+            ])
+
+            self.child(xf2.compose(self.splits[0]), self.depth, [c, m12, n+2, m23]);
+            self.child(xf2.compose(self.splits[1]), self.depth, [c, m01, n+1, m12]);
+            self.child(xf2.compose(self.splits[2]), self.depth, [c, m30, n+0, m01]);
+            self.child(xf2.compose(self.splits[3]), self.depth, [c, m23, n+3, m30]);

implicit.org

@@ -0,0 +1,22 @@
+* Slush-pile from implicit rant
+  (TODO: This was shelved from another post.  Do something with it.)
+  
+  None of my attempts have ever quite worked as soon as I start
+  dealing with implicit surfaces that have much sharpness to
+  them... which isn't saying much, since I am not exactly world-class
+  at dealing with meshes and mesh topology directly, but it's also
+  legitimately a very hard problem.  Even when I tried in CGAL, which
+  can handle implicit surfaces and does a far better job at generating
+  "good" meshes than anything I could hope to write myself, the
+  problem I ran into is that it was generating good meshes for
+  industrial CAD, not good meshes for /rendering/.  That is, their
+  inaccuracy to the true surface was bounded, and I could always
+  reduce it by just using denser meshes.  For CAD applications, this
+  seems fine. For rendering, though, I kept running into some major
+  inefficiencies in terms of number of triangles required (and as well
+  the amount of time required). I don't know much about the algorithms
+  they use, but I am fairly sure that for detail level N (inversely
+  proportional to minimum feature size) both triangle count, time, and
+  likely memory usage grow with ~O(N^3)~.
+
+  (TODO: Give an example of the above perhaps?)

libfive_subdiv/test_subdiv.py

@@ -0,0 +1,174 @@
+#!/usr/bin/env python3
+
+# Chris Hodapp, 2021-07-17
+#
+# This code is: yet another attempt at producing better meshes from
+# implicit surfaces / isosurfaces.  My paper notes from around the
+# same time period describe some more of why and how.
+#
+# This depends on the Python bindings for libfive (circa revision
+# 601730dc), on numpy, and on autograd from
+# https://github.com/HIPS/autograd for automatic differentiation.
+#
+# For an implicit surface expressed in a Python function, it:
+# - uses libfive to generate a mesh for this implicit surface,
+# - dumps this face-vertex data (numpy arrays) to disk in a form Blender
+#   can load pretty easily, (this is done only because exporting and
+#   loading an STL resulted in vertex and face indices being out of sync
+#   for some reason, perhaps libfive's meshing having randomness.)
+# - iterates over each edge from libfive's mesh data,
+# - for that edge, computes the curvature of the surface perpendicular
+#   to that edge,
+# - saves this curvature away in another file Blender can load.
+#
+# There are then some Blender routines for its Python API which load
+# the mesh, load the curvatures, and then try to turn these per-edge
+# curvature values to edge crease weights.  The hope was that this
+# would allow subdivision to work effectively on the resultant mesh in
+# sharper (higher-curvature) areas - lower crease weights should fit
+# lower-curvature areas better, and higher crease weights should keep
+# a sharper edge from being dulled too much by subdivision.
+# 
+# I tried with spiral_implicit, my same spiral isosurface function
+# from 2005 June yet again, as the implicit surface, but also yet
+# again, it proved a very difficult surface to work with.
+
+# Below is some elisp so that I can use the right environment in Emacs
+# and elpy:
+# 
+# (setq python-shell-interpreter "nix-shell" python-shell-interpreter-args " -I nixpkgs=/home/hodapp/nixpkgs -p python3Packages.libfive python3Packages.autograd python3Packages.numpy --command \"python3 -i\"")
+
+# This is a kludge to ensure libfive's bindings can be found:
+#import os, sys
+#os.environ["LIBFIVE_FRAMEWORK_DIR"]="/nix/store/gcxmz71b4i6bmsb1alafr4cqdnl19dn5-libfive-unstable-e93fef9d/lib/"
+#sys.path.insert(0, "/nix/store/gcxmz71b4i6bmsb1alafr4cqdnl19dn5-libfive-unstable-e93fef9d/lib/python3.8/site-packages/")
+
+import autograd.numpy as np
+from autograd import grad, elementwise_grad as egrad
+
+from libfive.shape import shape
+
+# The implicit surface is below.  It returns two functions that
+# compute the same thing: a vectorized version (f) that can handle
+# array inputs with (x,y,z) rows, and a version (g) that can also
+# handle individual x,y,z. f is needed for autograd, g is needed for
+# libfive.
+def spiral_implicit(outer, inner, freq, phase, thresh):
+    def g(x,y,z):
+        d1 = outer*y - inner*np.sin(freq*x + phase)
+        d2 = outer*z - inner*np.cos(freq*x + phase)
+        return d1*d1 + d2*d2 - thresh*thresh
+    def f(pt):
+        x,y,z = [pt[..., i] for i in range(3)]
+        return g(x,y,z)
+    return f, g
+
+def any_perpendicular(vecs):
+    # For 'vecs' of shape (..., 3), this returns an array of shape
+    # (..., 3) in which every corresponding vector is perpendicular
+    # (but nonzero).  'vecs' does not need to be normalized, and the
+    # returned vectors are not normalized.
+    x,y,z = [vecs[..., i] for i in range(3)]
+    a0 = np.zeros_like(x)
+    # The condition has the extra dimension added to make it (..., 1)
+    # so it broadcasts properly with the branches, which are (..., 3):
+    p = np.where((np.abs(z) < np.abs(x))[...,None],
+                 np.stack((y,  -x, a0), axis=-1),
+                 np.stack((a0, -z, y),  axis=-1))
+    return p
+
+def intersect_implicit(surface_fn):
+    # surface_fn(x,y,z)=0 is an implicit surface.  This returns a
+    # function f(s, t, pt, u, v) which - for f(s,t,...) = 0 is the
+    # implicit curve created by intersecting the surface with a plane
+    # passing through point 'pt' and with two perpendicular unit
+    # vectors 'u' and 'v' that lie on the plane.
+    def g(pts_2d, pt_center, u, v, **kw):
+        s,t = [pts_2d[..., i, None] for i in range(2)]
+        pt_3d = pt_center + s*u + t*v
+        return surface_fn(pt_3d, **kw)
+    return g
+
+def implicit_curvature_2d(curve_fn):
+    # Returns a function which computes curvature of an implicit
+    # curve, curve_fn(s,t)=0.  The resultant function takes two
+    # arguments as well.
+    #
+    # First derivatives:
+    _g1 = egrad(curve_fn)
+    # Second derivatives:
+    _g2s = egrad(lambda *a, **kw: _g1(*a, **kw)[...,0])
+    _g2t = egrad(lambda *a, **kw: _g1(*a, **kw)[...,1])
+    # Doing 'egrad' twice doesn't have the intended effect, so here I
+    # split up the first derivative manually.
+    def f(st, **kw):
+        g1  = _g1(st, **kw)
+        g2s = _g2s(st, **kw)
+        g2t = _g2t(st, **kw)
+        ds  = g1[..., 0]
+        dt  = g1[..., 1]
+        dss = g2s[..., 0]
+        dst = g2s[..., 1]
+        dtt = g2t[..., 1]
+        return (-dt*dt*dss + 2*ds*dt*dst - ds*ds*dtt) / ((ds*ds + dt*dt)**(3/2))
+    return f
+
+f_arr, f = spiral_implicit(2.0, 0.4, 20.0, 0.0, 0.3)
+fs = shape(f)
+print(fs)
+
+kw={
+    "xyz_min": (-0.5, -0.5, -0.5),
+    "xyz_max": (0.5, 0.5, 0.5),
+    "resolution": 20,
+}
+# To save directly as STL:
+# fs.save_stl("spiral.stl", **kw)
+
+print(f"letting libfive generate mesh...")
+verts, tris = fs.get_mesh(**kw)
+verts = np.array(verts, dtype=np.float32)
+tris = np.array(tris, dtype=np.uint32)
+
+print(f"Saving {len(verts)} vertices, {len(tris)} faces")
+np.save("spiral_verts.npy", verts)
+np.save("spiral_tris.npy", tris)
+
+print(f"Computing curvatures...")
+
+# Shape (N, 3, 3). Final axis is (x,y,z).
+tri_verts = verts[tris]
+# Compute all 3 midpoints (over each edge):
+v_pairs = [(tri_verts[:, i, :], tri_verts[:, (i+1)%3, :])
+           for i in range(3)]
+print(f"midpoints")
+tri_mids = np.stack([(vi+vj)/2 for vi,vj in v_pairs],
+                    axis=1)
+print(f"edge vectors")
+# Compute normalized edge vectors:
+diff = [vj-vi for vi,vj in v_pairs]
+edge_vecs = np.stack([d/np.linalg.norm(d, axis=1, keepdims=True) for d in diff],
+                     axis=1)
+print(f"perpendiculars")
+# Find perpendicular to all edge vectors:
+v1 = any_perpendicular(edge_vecs)
+v1 /= np.linalg.norm(v1, axis=-1, keepdims=True)
+# and perpendiculars to both:
+v2 = np.cross(edge_vecs, v1)
+
+print(f"implicit curves")
+isect_2d = intersect_implicit(f_arr)
+curv_fn = implicit_curvature_2d(isect_2d)
+print(f"gradients & curvature")
+k = curv_fn(np.zeros((tri_mids.shape[0], 3, 2)), pt_center=tri_mids, u=v1, v=v2)
+
+print(f"writing")
+np.save("spiral_curvature.npy", k)
+
+# for i,k_i in enumerate(k):
+#     for j in range(k.shape[1]):
+#         mid = tri_mids[i, j, :]
+#         k_ij = k[i,j]
+#         v1 = tris[i][j]
+#         v2 = tris[i][(j + 1) % 3]
+#         print(f"{i}: {v1} to {v2}, {k_ij:.3f}")

parallel_transport/parallel_transport.py

@@ -0,0 +1,84 @@
+import numpy
+import sympy
+from sympy.vector import CoordSys3D
+
+def mtx_rotate_by_vector(b, theta):
+    """Returns 3x3 matrix for rotating around some 3D vector."""
+    # Source:
+    # Appendix A of "Parallel Transport to Curve Framing"
+    c = numpy.cos(theta)
+    s = numpy.sin(theta)
+    rot = numpy.array([
+        [c+(b[0]**2)*(1-c), b[0]*b[1]*(1-c)-s*b[2], b[2]*b[0]*(1-c)+s*b[1]],
+        [b[0]*b[1]*(1-c)+s*b[2], c+(b[1]**2)*(1-c), b[2]*b[1]*(1-c)-s*b[0]],
+        [b[0]*b[2]*(1-c)-s*b[1], b[1]*b[2]*(1-c)+s*b[0], c+(b[2]**2)*(1-c)],
+    ])
+    return rot
+
+def gen_faces(v1a, v1b, v2a, v2b):
+    """Returns faces (as arrays of vertices) connecting two pairs of
+    vertices."""
+    # Keep winding order consistent!
+    f1 = numpy.array([v1b, v1a, v2a])
+    f2 = numpy.array([v2b, v1b, v2a])
+    return f1, f2
+
+def approx_tangents(points):
+    """Returns an array of approximate tangent vectors.  Assumes a
+    closed path, and approximates point I using neighbors I-1 and I+1 -
+    that is, treating the three points as a circle.
+    
+    Input:
+    points -- Array of shape (N,3). points[0,:] is assumed to wrap around
+              to points[-1,:].
+    
+    Output:
+    tangents -- Array of same shape as 'points'. Each row is normalized.
+    """
+    d = numpy.roll(points, -1, axis=0) - numpy.roll(points, +1, axis=0)
+    d = d/numpy.linalg.norm(d, axis=1)[:,numpy.newaxis]
+    return d
+
+def approx_arc_length(points):
+    p2 = numpy.roll(points, -1, axis=0)
+    return numpy.sum(numpy.linalg.norm(points - p2, axis=1))
+
+def torsion(v, arg):
+    """Returns an analytical SymPy expression for torsion of a 3D curve.
+    
+    Inputs:
+    v -- SymPy expression returning a 3D vector
+    arg -- SymPy symbol for v's variable
+    """
+    # https://en.wikipedia.org/wiki/Torsion_of_a_curve#Alternative_description
+    dv1 = v.diff(arg)
+    dv2 = dv1.diff(arg)
+    dv3 = dv2.diff(arg)
+    v1_x_v2 = dv1.cross(dv2)
+    # This calls for the square of the norm in denominator - but that
+    # is just dot product with itself:
+    return v1_x_v2.dot(dv3) / (v1_x_v2.dot(v1_x_v2))
+
+def torsion_integral(curve_expr, var, a, b):
+    # The line integral from section 3.1 of "Parallel Transport to Curve
+    # Framing".  This should work in theory, but with the functions I've
+    # actually tried, evalf() is ridiculously slow.
+    c = torsion(curve_expr, var)
+    return sympy.Integral(c * (sympy.diff(curve_expr, var).magnitude()), (var, a, b))
+
+def torsion_integral_approx(curve_expr, var, a, b, step):
+    # A numerical approximation of the line integral from section 3.1 of
+    # "Parallel Transport to Curve Framing"
+    N = CoordSys3D('N')
+    # Get a (callable) derivative function of the curve:
+    curve_diff = curve_expr.diff(var)
+    diff_fn = sympy.lambdify([var], N.origin.locate_new('na', curve_diff).express_coordinates(N))
+    # And a torsion function:
+    torsion_fn = sympy.lambdify([var], torsion(curve_expr, var))
+    # Generate values of 'var' to use:
+    vs = numpy.arange(a, b, step)
+    # Evaluate derivative function & torsion function over these:
+    d = numpy.linalg.norm(numpy.array(diff_fn(vs)).T, axis=1)
+    torsions = torsion_fn(vs)
+    # Turn this into basically a left Riemann sum (I'm lazy):
+    return -(d * torsions * step).sum()

parallel_transport/quat.py

@@ -0,0 +1,24 @@
+import numpy
+import quaternion
+
+def conjugate_by(vec, quat):
+    """Turn 'vec' to a quaternion, conjugate it by 'quat', and return it."""
+    q2 = quat * vec2quat(vec) * quat.conjugate()
+    return quaternion.as_float_array(q2)[:,1:]
+
+def rotation_quaternion(axis, angle):
+    """Returns a quaternion for rotating by some axis and angle.
+    
+    Inputs:
+    axis -- numpy array of shape (3,), with axis to rotate around
+    angle -- angle in radians by which to rotate
+    """
+    qc = numpy.cos(angle / 2)
+    qs = numpy.sin(angle / 2)
+    qv = qs * axis
+    return numpy.quaternion(qc, qv[0], qv[1], qv[2])
+
+def vec2quat(vs):
+    qs = numpy.zeros(vs.shape[0], dtype=numpy.quaternion)
+    quaternion.as_float_array(qs)[:,1:4] = vs
+    return qs

parallel_transport/spiral_parametric3.py

@@ -0,0 +1,232 @@
+#!/usr/bin/env python3
+
+import sys
+import numpy
+import stl.mesh
+
+# TODO:
+# - This still has some strange errors around high curvature.
+# They are plainly visible in the cross-section.
+# - Check rotation direction
+# - Fix phase, which only works if 0 (due to how I work with y)
+# Things don't seem to line up right.
+# - Why is there still a gap when using Array modifier?
+# Check beginning and ending vertices maybe
+# - Organize this so that it generates both meshes when run
+
+# This is all rather tightly-coupled.  Almost everything is specific
+# to the isosurface I was trying to generate. walk_curve may be able
+# to generalize to some other shapes.
+class ExplicitSurfaceThing(object):
+    def __init__(self, freq, phase, scale, inner, outer, rad, ext_phase):
+        self.freq = freq
+        self.phase = phase
+        self.scale = scale
+        self.inner = inner
+        self.outer = outer
+        self.rad = rad
+        self.ext_phase = ext_phase
+        
+    def angle(self, z):
+        return self.freq*z + self.phase
+
+    def max_z(self):
+        # This value is the largest |z| for which 'radical' >= 0
+        # (thus, for x_cross to have a valid solution)
+        return (numpy.arcsin(self.rad / self.inner) - self.phase) / self.freq
+
+    def radical(self, z):
+        return self.rad*self.rad - self.inner*self.inner * (numpy.sin(self.angle(z)))**2
+
+    # Implicit curve function
+    def F(self, x, z):
+        return (self.outer*x - self.inner*numpy.cos(self.angle(z)))**2 + (self.inner*numpy.sin(self.angle(z)))**2 - self.rad**2
+
+    # Partial 1st derivatives of F:
+    def F_x(self, x, z):
+        return 2 * self.outer * self.outer * x - 2 * self.inner * self.outer * numpy.cos(self.angle(z))
+
+    def F_z(self, x, z):
+        return 2 * self.freq * self.inner * self.outer * numpy.sin(self.angle(z))
+
+    # Curvature:
+    def K(self, x, z):
+        a1 = self.outer**2
+        a2 = x**2
+        a3 = self.freq*z + self.phase
+        a4 = numpy.cos(a3)
+        a5 = self.inner**2
+        a6 = a4**2
+        a7 = self.freq**2
+        a8 = numpy.sin(a3)**2
+        a9 = self.outer**3
+        a10 = self.inner**3
+        return -((2*a7*a10*self.outer*x*a4 + 2*a7*a5*a1*a2)*a8 + (2*a7*self.inner*a9*x**3 + 2*a7*a10*self.outer*x)*a4 - 4*a7*a5*a1*a2) / ((a7*a5*a2*a8 + a5*a6 - 2*self.inner*self.outer*x*a4 + a1*a2) * numpy.sqrt(4*a7*a5*a1*a2*a8 + 4*a5*a1*a6 - 8*self.inner*a9*x*a4 + 4*a2*self.outer**4))
+
+    def walk_curve(self, x0, z0, eps, thresh = 1e-3, gd_thresh = 1e-7):
+        x, z = x0, z0
+        eps2 = eps*eps
+        verts = []
+        iters = 0
+        # Until we return to the same point at which we started...
+        while True:
+            iters += 1
+            verts.append([x, 0, z])
+            # ...walk around the curve by stepping perpendicular to the
+            # gradient by 'eps'.  So, first find the gradient:
+            dx = self.F_x(x, z)
+            dz = self.F_z(x, z)
+            # Normalize it:
+            f = 1/numpy.sqrt(dx*dx + dz*dz)
+            nx, nz = dx*f, dz*f
+            # Find curvature at this point because it tells us a little
+            # about how far we can safely move:
+            K_val = abs(self.K(x, z))
+            eps_corr = 2 * numpy.sqrt(2*eps/K_val - eps*eps)
+            # Scale by 'eps' and use (-dz, dx) as perpendicular:
+            px, pz = -nz*eps_corr, nx*eps_corr
+            # Walk in that direction:
+            x += px
+            z += pz
+            # Moving in that direction is only good locally, and we may
+            # have deviated off the curve slightly.  The implicit function
+            # tells us (sort of) how far away we are, and the gradient
+            # tells us how to minimize that:
+            #print("W: x={} z={} dx={} dz={} px={} pz={} K={} eps_corr={}".format(
+            #    x, z, dx, dz, px, pz, K_val, eps_corr))
+            F_val = self.F(x, z)
+            count = 0
+            while abs(F_val) > gd_thresh:
+                count += 1
+                dx = self.F_x(x, z)
+                dz = self.F_z(x, z)
+                f = 1/numpy.sqrt(dx*dx + dz*dz)
+                nx, nz = dx*f, dz*f
+                # If F is negative, we want to increase it (thus, follow
+                # gradient).  If F is positive, we want to decrease it
+                # (thus, opposite of gradient).
+                F_val = self.F(x, z)
+                x += -F_val*nx
+                z += -F_val*nz
+                # Yes, this is inefficient gradient-descent...
+            diff = numpy.sqrt((x-x0)**2 + (z-z0)**2)
+            #print("{} gradient-descent iters. diff = {}".format(count, diff))
+            if iters > 100 and diff < thresh:
+                #print("diff < eps, quitting")
+                #verts.append([x, 0, z])
+                break
+        data = numpy.array(verts)
+        return data
+
+    def x_cross(self, z, sign):
+        # Single cross-section point in XZ for y=0.  Set sign for positive
+        # or negative solution.
+        n1 = numpy.sqrt(self.radical(z))
+        n2 = self.inner * numpy.cos(self.angle(z))
+        if sign > 0:
+            return (n2-n1) / self.outer
+        else:
+            return (n2+n1) / self.outer
+
+    def turn(self, points, dz):
+        # Note one full revolution is dz = 2*pi/freq
+        # How far to turn in radians (determined by dz):
+        rad = self.angle(dz)
+        c, s = numpy.cos(rad), numpy.sin(rad)
+        mtx = numpy.array([
+            [ c,  s,  0],
+            [-s,  c,  0],
+            [ 0,  0,  1],
+        ])
+        return points.dot(mtx) + [0, 0, dz]
+
+    def screw_360(self, z0_period_start, x_init, z_init, eps, dz, thresh, endcaps=False):
+        #z0 = -10 * 2*numpy.pi/freq / 2
+        z0 = z0_period_start * 2*numpy.pi/self.freq / 2
+        z1 = z0 + 2*numpy.pi/self.freq
+        #z1 = 5 * 2*numpy.pi/freq / 2
+        #z0 = 0
+        #z1 = 2*numpy.pi/freq
+        init_xsec = self.walk_curve(x_init, z_init, eps, thresh)
+        num_xsec_steps = init_xsec.shape[0]
+        zs = numpy.append(numpy.arange(z0, z1, dz), z1)
+        num_screw_steps = len(zs)
+        vecs = num_xsec_steps * num_screw_steps * 2
+        offset = 0
+        if endcaps:
+            offset = num_xsec_steps
+            vecs += 2*num_xsec_steps
+        print("Generating {} vertices...".format(vecs))
+        data = numpy.zeros(vecs, dtype=stl.mesh.Mesh.dtype)
+        v = data["vectors"]
+        # First endcap:
+        if endcaps:
+            center = init_xsec.mean(0)
+            for i in range(num_xsec_steps):
+                v[i][0,:] = init_xsec[(i + 1) % num_xsec_steps,:]
+                v[i][1,:] = init_xsec[i,:]
+                v[i][2,:] = center
+        # Body:
+        verts = init_xsec
+        for i,z in enumerate(zs):
+            verts_last = verts
+            verts = self.turn(init_xsec, z-z0)
+            if i > 0:
+                for j in range(num_xsec_steps):
+                    # Vertex index:
+                    vi = offset + (i-1)*num_xsec_steps*2 + j*2
+                    v[vi][0,:] = verts[(j + 1) % num_xsec_steps,:]
+                    v[vi][1,:] = verts[j,:]
+                    v[vi][2,:] = verts_last[j,:]
+                    #print("Write vertex {}".format(vi))
+                    v[vi+1][0,:] = verts_last[(j + 1) % num_xsec_steps,:]
+                    v[vi+1][1,:] = verts[(j + 1) % num_xsec_steps,:]
+                    v[vi+1][2,:] = verts_last[j,:]
+                    #print("Write vertex {} (2nd half)".format(vi+1))
+        # Second endcap:
+        if endcaps:
+            center = verts.mean(0)
+            for i in range(num_xsec_steps):
+                vi = num_xsec_steps * num_screw_steps * 2 + num_xsec_steps + i
+                v[vi][0,:] = center
+                v[vi][1,:] = verts[i,:]
+                v[vi][2,:] = verts[(i + 1) % num_xsec_steps,:]
+        v[:, :, 2] += z0 + self.ext_phase / self.freq
+        v[:, :, :] /= self.scale
+        mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
+        print("Beginning z: {}".format(z0/self.scale))
+        print("Ending z: {}".format(z1/self.scale))
+        print("Period: {}".format((z1-z0)/self.scale))
+        return mesh
+
+surf1 = ExplicitSurfaceThing(
+    freq = 20,
+    phase = 0,
+    scale = 1/16, # from libfive
+    inner = 0.4 * 1/16,
+    outer = 2.0 * 1/16,
+    rad = 0.3 * 1/16,
+    ext_phase = 0)
+
+z_init = 0
+x_init = surf1.x_cross(z_init, 1)
+mesh1 = surf1.screw_360(-10, x_init, z_init, 0.000002, 0.001, 5e-4)
+fname = "spiral_inner0_one_period.stl"
+print("Writing {}...".format(fname))
+mesh1.save(fname)
+
+surf2 = ExplicitSurfaceThing(
+    freq = 10,
+    phase = 0,
+    scale = 1/16, # from libfive
+    inner = 0.9 * 1/16,
+    outer = 2.0 * 1/16,
+    rad = 0.3 * 1/16,
+    ext_phase = numpy.pi/2)
+
+z_init = 0
+x_init = surf2.x_cross(z_init, 1)
+mesh2 = surf2.screw_360(-5, x_init, z_init, 0.000002, 0.001, 5e-4)
+fname = "spiral_outer90_one_period.stl"
+print("Writing {}...".format(fname))
+mesh2.save(fname)

python_extrude_meshgen/README.md

@@ -0,0 +1,84 @@
+# To-do items, wanted features, bugs:
+
+## Cool
+- More complicated: Examples of *merging*. I'm not sure on the theory
+  behind this.
+  
+## Annoying/boring
+- https://en.wikipedia.org/wiki/Polygon_triangulation - do this to
+  fix my wave example!
+  - http://www.polygontriangulation.com/2018/07/triangulation-algorithm.html
+- Clean up `examples.ram_horn_branch()`. The way I clean it up might
+  help inform some cleaner designs.
+- I really need to standardize some of the behavior of fundamental
+  operations (with regard to things like sizes they generate). This
+  is behavior that, if it changes, will change a lot of things that I'm
+  trying to keep consistent so that my examples still work.
+- Winding order.  It is consistent through seemingly everything,
+  except for reflection and `close_boundary_simple`.  (When there are
+  two parallel boundaries joined with something like
+  `join_boundary_simple`, traversing these boundaries in their actual
+  order to generate triangles - like in `close_boundary_simple` - will
+  produce opposite winding order on each. Imagine a transparent clock:
+  seen from the front, it moves clockwise, but seen from the back, it
+  moves counter-clockwise.)
+- File that bug that I've seen in trimesh/three.js
+  (see trimesh_fail.ipynb)
+- Why do I get the weird zig-zag pattern on the triangles,
+  despite larger numbers of them? Is it something in how I
+  twist the frames?
+  - How can I compute the *torsion* on a quad? I think it
+    comes down to this: torsion applied across the quad I'm
+    triangulating leading to neither diagonal being a
+    particularly good choice.  Subdividing the boundary seems
+    to help, but other triangulation methods (e.g. turning a
+    quad to 4 triangles by adding the centroid) could be good
+    too.
+  - Facets/edges are just oriented the wrong way...
+  - Picking at random the diagonal on the quad to triangulate with
+    does seem to turn 'error' just to noise, and in its own way this
+    is preferable.
+- Integrate parallel_transport work and reuse what I can
+- `/mnt/dev/graphics_misc/isosurfaces_2018_2019` - perhaps include my
+  spiral isosurface stuff from here
+
+## Abstractions
+
+- This has a lot of functions parametrized over a lot
+  of functions.  Need to work with this somehow.  (e.g. should
+  it subdivide this boundary? should it merge opening/closing
+  boundaries?)
+- Some generators produce boundaries that can be directly merged
+  and produce sensible geometry.  Some generators produce
+  boundaries that are only usable when they are further
+  transformed (and would produce degenerate geometry).  What sort
+  of nomenclature captures this?
+
+- How can I capture the idea of a group of parameters which, if
+  they are all scaled in the correct way (some linearly, others
+  inversely perhaps), generated geometry that is more or less
+  identical except that it is higher-resolution?
+- Use mixins to extend 3D transformations to things (matrices,
+  cages, meshes, existing transformations)
+- I can transform a Cage.  Why not a CageGen?
+  
+## ????
+- Embed this in Blender?
+  
+## Future thoughts
+
+- What if I had a function that could generate a Cage as if
+  from a parametric formula and smoothly vary its orientation?
+  My existing tools could easily turn this to a mesh. If I could vary
+  the detail of the Cage itself (if needed), then I could also
+  generate a mesh at an arbitrary level of detail simply by sampling at
+  finer and finer points on the parameter space.  (This might also tie
+  into the Parallel Transport work.)
+- What are the limitations of using Cages?
+- Current system is very "generative".  Could I do basically L-system
+  if I have rules for how a much is *refined*?  What about IFS?
+- Do this in Rust once I understand WTF I am doing
+
+## Other thoughts
+
+- Why do I never use the term "extruding" to describe what I'm doing?

python_isosurfaces_2018_2019/build_osx.sh

@@ -0,0 +1,10 @@
+#!/bin/sh
+
+g++ -L/usr/local/Cellar/cgal/4.12/lib -I/usr/local/Cellar/cgal/4.12/include \
+    -L/usr/local/Cellar/gmp/6.1.2_2/lib \
+    -L/usr/local/Cellar/mpfr/4.0.1/lib \
+    -L/usr/local/Cellar/boost/1.67.0_1/lib \
+    -DCGAL_CONCURRENT_MESH_3 \
+    -lCGAL -lgmp -lmpfr -lboost_thread-mt \
+    ./mesh_an_implicit_function.cpp \
+    -o mesh_an_implicit_function.o

python_isosurfaces_2018_2019/cgal_dabbling.org

@@ -0,0 +1,173 @@
+#+TITLE: CGAL dabbling
+#+DATE: <2018-08-06 Mon>
+#+AUTHOR: Hodapp
+#+EMAIL: hodapp87@gmail.com
+#+OPTIONS: ':nil *:t -:t ::t <:t H:3 \n:nil ^:t arch:headline
+#+OPTIONS: author:t c:nil creator:comment d:(not "LOGBOOK") date:t
+#+OPTIONS: e:t email:nil f:t inline:t num:t p:nil pri:nil stat:t
+#+OPTIONS: tags:t tasks:t tex:t timestamp:t toc:t todo:t |:t
+#+DESCRIPTION:
+#+EXCLUDE_TAGS: noexport
+#+KEYWORDS:
+#+LANGUAGE: en
+#+SELECT_TAGS: export
+
+# By default I do not want that source code blocks are evaluated on export. Usually
+# I want to evaluate them interactively and retain the original results.
+#+PROPERTY: header-args :eval never-export :export both
+
+- CGAL is one of the most insanely cryptic and impenetrable libraries
+  I have found.
+- Where I am stuck now:
+  - I can use 1D features in
+    [[https://doc.cgal.org/latest/Mesh_3/Mesh_3_2mesh_two_implicit_spheres_with_balls_8cpp-example.html][Mesh_3/mesh_two_implicit_spheres_with_balls.cpp]] but this is the
+    wrong sort of data (it's a 3D mesh, yes, but not a surface mesh)
+    and "medit" and "vtu" are all it can write.  How do I extract a
+    surface mesh to get a Polyhedron_3 so that I can write an OBJ?
+  - Or: Is there any way to use 1D features with a surface mesh?
+  - Even if I manage to put this Mesh_3 into the right form I have no
+    guarantees that it is a good surface mesh.  It uses a different
+    algorithm than for the surface mesh - I'm simply trying to take
+    the surface of the result.  I should also expect this algorithm to
+    be something more like O(N^3) with mesh resolution because it must
+    fill the volume with tetrahedrons, and I then just throw away all
+    of these.
+- This is based around the source for [[https://doc.cgal.org/latest/Mesh_3/Mesh_3_2mesh_implicit_sphere_8cpp-example.html][mesh_implicit_sphere.cpp]] from
+  [[https://doc.cgal.org/latest/Mesh_3/index.html][CGAL: 3D Mesh Generation]].
+- Why am I not using [[https://doc.cgal.org/latest/Mesh_3/group__PkgMesh__3Functions.html#ga68ca38989087644fb6b893eb566be9ea][facets_in_complex_3_to_triangle_mesh()]]?
+- [[https://doc.cgal.org/latest/Mesh_3/Mesh_3_2mesh_two_implicit_spheres_with_balls_8cpp-example.html][Mesh_3/mesh_two_implicit_spheres_with_balls.cpp]] shows the use of
+  [[https://doc.cgal.org/latest/Mesh_3/classCGAL_1_1Mesh__domain__with__polyline__features__3.html][Mesh_domain_with_polyline_features_3]] to explicitly give features to
+  preserve
+  - This uses the ~make_mesh_3~ call while I'm using
+    ~make_surface_mesh~...
+  - I can only call ~CGAL::print_polyhedron_wavefront~ if I have a
+    ~Polyhedron~ and I cannot figure out how to get one from a ~C3t3~.
+    So, perhaps I am stuck with this "medit" mesh format, or else VTU.
+  - I think the terminology is trying to tell me: this isn't a surface
+    mesh, it's a 3D mesh made of tetrahedra.  It sounds like the type
+    of data from a Delaunay triangulation is inherently rather
+    different from a surface mesh.
+  - https://doc.cgal.org/latest/Mesh_3/index.html#title24 does the
+    opposite direction of what I need
+- There are [[https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__detect__features__grp.html][Feature Detection Functions]], but ~detect_features~ is in
+  [[https://doc.cgal.org/latest/Mesh_3/classCGAL_1_1Polyhedral__complex__mesh__domain__3.html][Polyhedral_complex_mesh_domain_3]] and
+  [[https://doc.cgal.org/latest/Mesh_3/classCGAL_1_1Polyhedral__mesh__domain__with__features__3.html#a5a868ac7b8673436766d28b7a80d2826][Polyhedral_mesh_domain_with_features_3]]
+- Domains:
+  - [[https://doc.cgal.org/latest/Mesh_3/classMeshDomain__3.html][MeshDomain_3]] concept
+  - [[https://doc.cgal.org/latest/Mesh_3/classMeshDomainWithFeatures__3.html][MeshDomainWithFeatures_3]] concept
+- Why make_surface_mesh, Implicit_surface_3,
+  Complex_2_in_triangulation_3 vs. make_mesh_3,
+  Implicit_mesh_domain_3, Mesh_complex_3_in_triangulation_3?
+- [[https://doc.cgal.org/latest/Mesh_3/group__PkgMesh3Parameters.html#ga0a990b28d55157c62d4bfd2624d408af][parameters::features()]] - can do 1-dimensional features
+  - now how do I use it?
+- [[https://doc.cgal.org/latest/Surface_mesher/index.html][make_surface_mesh]] args:
+  - ~SurfaceMeshC2T3& c2t3~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurfaceMeshComplex__2InTriangulation__3.html][SurfaceMeshComplex_2InTriangulation_3]]
+  - ~SurfaceMeshTraits::Surface_3 surface~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurface__3.html][Surface_3]]
+  - ~SurfaceMeshTraits traits~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurfaceMeshTraits__3.html][SurfaceMeshTraits_3]]
+  - ~FacetsCriteria criteria~ - must be a model of concept [[https://doc.cgal.org/latest/Surface_mesher/classSurfaceMeshFacetsCriteria__3.html][SurfaceMeshFacetsCriteria_3]]
+  - ~Tag~
+- "This algorithm of CGAL::make_surface_mesh is designed for smooth
+  implicit surfaces. If your implicit surface is not smooth, then the
+  sharp features of the surface will not be meshed correctly."
+
+#+BEGIN_SRC elisp
+  (setq org-confirm-babel-evaluate nil)
+  (setq org-src-fontify-natively t)
+  (setq org-src-tab-acts-natively t)
+  (org-version)
+#+END_SRC
+
+#+RESULTS:
+: 9.1.9
+
+#+BEGIN_SRC sh :results verbatim
+  gcc --version
+#+END_SRC
+
+#+RESULTS:
+: Apple LLVM version 9.1.0 (clang-902.0.39.2)
+: Target: x86_64-apple-darwin17.7.0
+: Thread model: posix
+: InstalledDir: /Library/Developer/CommandLineTools/usr/bin
+
+#+NAME: includes
+#+BEGIN_SRC C++
+  #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+  #include <CGAL/Mesh_triangulation_3.h>
+  #include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+  #include <CGAL/Mesh_criteria_3.h>
+
+  #include <CGAL/Implicit_mesh_domain_3.h>
+  #include <CGAL/make_mesh_3.h>
+#+END_SRC
+
+#+NAME: typesDomain
+#+BEGIN_SRC C++
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+  typedef K::FT FT;
+  typedef K::Point_3 Point;
+  typedef FT (Function)(const Point&);
+  typedef CGAL::Implicit_mesh_domain_3<Function,K> Mesh_domain;
+
+  #ifdef CGAL_CONCURRENT_MESH_3
+  typedef CGAL::Parallel_tag Concurrency_tag;
+  #else
+  typedef CGAL::Sequential_tag Concurrency_tag;
+  #endif
+#+END_SRC
+
+#+NAME: typesTriangulation
+#+BEGIN_SRC C++
+  typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
+  typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
+#+END_SRC
+
+#+NAME: typesCriteria
+#+BEGIN_SRC C++
+  typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+#+END_SRC
+
+#+NAME: sphereFunction
+#+BEGIN_SRC C++
+  FT sphere_function (const Point& p) {
+      return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-1;
+  }
+#+END_SRC
+
+#+BEGIN_SRC C++ :noweb yes :flags  -L/usr/local/Cellar/cgal/4.12/lib -I/usr/local/Cellar/cgal/4.12/include  -L/usr/local/Cellar/gmp/6.1.2_2/lib  -L/usr/local/Cellar/mpfr/4.0.1/lib  -L/usr/local/Cellar/boost/1.67.0_1/lib  -DCGAL_CONCURRENT_MESH_3  -lCGAL -lgmp -lmpfr -lboost_thread-mt
+  <<includes>>
+
+  using namespace CGAL::parameters;
+
+  <<typesDomain>>
+  <<typesTriangulation>>
+  <<typesCriteria>>
+  <<sphereFunction>>
+
+  int main()
+  {
+      // Domain (Warning: Sphere_3 constructor uses squared radius !)
+      Mesh_domain domain(sphere_function,
+                         K::Sphere_3(CGAL::ORIGIN, 2.));
+      // Mesh criteria
+      Mesh_criteria criteria(facet_angle=30, facet_size=0.1, facet_distance=0.025,
+                             cell_radius_edge_ratio=2, cell_size=0.1);
+
+      std::cout << "Generating..." << std::endl;
+      // Mesh generation
+      C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+
+      // Output
+      std::ofstream medit_file("out.mesh");
+      c3t3.output_to_medit(medit_file);
+
+      std::cout << "Done" << std::endl;
+
+      return 0;
+  }
+#+END_SRC
+
+#+RESULTS:
+| Generating... |
+| Done          |

python_isosurfaces_2018_2019/mesh_an_implicit_function.cpp

@@ -0,0 +1,171 @@
+// Taken from:
+// https://doc.cgal.org/latest/Surface_mesher/Surface_mesher_2mesh_an_implicit_function_8cpp-example.html
+// https://doc.cgal.org/latest/Surface_mesher/index.html
+// https://doc.cgal.org/latest/Mesh_3/index.html
+
+#include <CGAL/Surface_mesh_default_triangulation_3.h>
+#include <CGAL/Complex_2_in_triangulation_3.h>
+#include <CGAL/IO/Complex_2_in_triangulation_3_file_writer.h>
+#ifdef CGAL_FACETS_IN_COMPLEX_2_TO_TRIANGLE_MESH_H
+  #include <CGAL/IO/facets_in_complex_2_to_triangle_mesh.h>
+#else
+  // NixOS currently has CGAL 4.11, not 4.12.  I guess 4.12 is needed
+  // for this.  The #ifdef is unnecessary, but the header and call for
+  // below are deprecated.
+  #include <CGAL/IO/output_surface_facets_to_polyhedron.h>
+#endif
+#include <CGAL/make_surface_mesh.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Implicit_surface_3.h>
+#include <CGAL/IO/print_wavefront.h>
+#include <CGAL/Polyhedron_3.h>
+
+#include <iostream>
+#include <fstream>
+#include <limits>
+#include <algorithm>
+
+// Triangulation
+typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
+typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
+
+// Domain?
+typedef Tr::Geom_traits GT;
+typedef GT::Sphere_3 Sphere_3;
+typedef GT::Point_3 Point_3;
+typedef GT::Vector_3 Vector_3;
+typedef GT::FT FT;
+typedef FT (*Function)(Point_3);
+typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;
+// how does this differ from CGAL::Implicit_mesh_domain_3<Function,K>?
+
+typedef CGAL::Polyhedron_3<GT> Polyhedron;
+
+FT sphere_function(Point_3 p) {
+    Point_3 p2(p.x() + 0.1 * cos(p.x() * 20),
+               p.y(),
+               p.z() + 0.1 * sin(p.z() * 4));
+    const FT x2=p2.x()*p2.x(), y2=p2.y()*p2.y(), z2=p2.z()*p2.z();
+    return x2+y2+z2-1;
+}
+
+Vector_3 sphere_gradient(Point_3 p) {
+	float A = 0.1;
+	float B = 0.1;
+	float F1 = 20;
+	float F2 = 4;
+	return Vector_3(2*(A*cos(p.x()*F1) + p.x())*(1 - A*F1*sin(p.x()*F1)),
+					2*p.y(),
+					2*(B*sin(p.z()*F2) + p.z())*(1 + B*F2*cos(p.z()*F2)));
+}
+
+FT spiral_function(Point_3 p) {
+    float outer = 2.0;
+    float inner = 0.4; // 0.9
+    float freq = 20; // 10
+    float phase = M_PI; // 3 * M_PI / 2;
+    float thresh = 0.3;
+    const FT d1 = p.y()*outer - inner * sin(p.x()*freq + phase);
+    const FT d2 = p.z()*outer - inner * cos(p.x()*freq + phase);
+    return std::max(sqrt(d1*d1 + d2*d2) - thresh,
+                    p.x()*p.x() + p.y()*p.y() + p.z()*p.z() - 1.9*1.9);
+}
+
+Vector_3 spiral_gradient(Point_3 p) {
+    float outer = 2.0;
+    float inner = 0.4;
+    float freq = 20;
+    float phase = M_PI;
+    float thresh = 0.3;
+	// "block([%1,%2,%3,%4,%5,%6],%1:P+x*F,%2:cos(%1),%3:z*O-I*%2,%4:sin(%1),%5:y*O-I*%4,%6:1/sqrt(%5^2+%3^2),[((2*F*I*%3*%4-2*F*I*%2*%5)*%6)/2,O*%5*%6,O*%3*%6])"
+	float v1 = phase + p.x() * freq;
+	float v2 = cos(v1);
+	float v3 = p.z() * outer - inner * v2;
+	float v4 = sin(v1);
+	float v5 = p.y() * outer - inner * v4;
+	float v6 = 1.0 / sqrt(v5*v5 + v3*v3);
+	return Vector_3(((2*freq*inner*v3*v4-2*freq*inner*v2*v5)*v6)/2,
+					outer * v5 * v6,
+					outer * v3 * v6);
+}
+
+int main() {
+    Tr tr;            // 3D-Delaunay triangulation
+    C2t3 c2t3 (tr);   // 2D-complex in 3D-Delaunay triangulation
+
+    FT bounding_sphere_rad = 2.0;
+    
+    // defining the surface
+    Surface_3 surface(spiral_function,             // pointer to function
+                      Sphere_3(CGAL::ORIGIN,
+                               bounding_sphere_rad*bounding_sphere_rad)); // bounding sphere
+
+    std::string fname("spiral_thing4.obj");
+    float angular_bound = 30;
+    float radius_bound = 0.02;
+    float distance_bound = 0.02;
+    
+    // Note that "2." above is the *squared* radius of the bounding sphere!
+    // defining meshing criteria
+    CGAL::Surface_mesh_default_criteria_3<Tr> criteria(
+        angular_bound, radius_bound, distance_bound);
+
+    std::cout << "angular bound = " << angular_bound << ", "
+              << "radius bound = " << radius_bound << ", "
+              << "distance bound = " << distance_bound << std::endl;
+
+    std::cout << "Making surface mesh..." << std::endl;
+    // meshing surface
+    CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Manifold_tag());
+    std::cout << "Vertices: " << tr.number_of_vertices() << std::endl;
+
+    // This didn't work on some calls instead of 'poly':
+    //CGAL::Surface_mesh<Point_3> sm;
+    Polyhedron poly;                                        
+    std::cout << "Turning facets to triangle mesh..." << std::endl;
+    
+#ifdef CGAL_FACETS_IN_COMPLEX_2_TO_TRIANGLE_MESH_H
+    CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, poly);
+#else
+    CGAL::output_surface_facets_to_polyhedron(c2t3, poly);
+#endif
+
+    FT err = 0.0;
+	FT inf = std::numeric_limits<FT>::infinity();
+    for (Polyhedron::Point_iterator it = poly.points_begin();
+         it != poly.points_end();
+         ++it)
+	{
+
+		FT rate = 2e-6;
+		FT f0 = abs(spiral_function(*it));
+		FT f;
+		for (int i = 0; i < 100; ++i) {
+			f = spiral_function(*it);
+			Vector_3 grad(spiral_gradient(*it));
+
+			*it -= grad * rate * (f > 0 ? 1 : -1);
+			/*
+			std::cout << "Iter " << i << ": "
+					  << "F(" << it->x() << "," << it->y() << "," << it->z()
+					  << ")=" << f << ", F'=" << grad << std::endl;
+			*/
+		}
+		//FT diff = (abs(f) - abs(f0)) / f0;
+        /*
+		std::cout << "F(" << it->x() << "," << it->y() << "," << it->z()
+				  << "): " << f0 << " to " << f << std::endl;
+        */
+		
+        err += f * f;
+    }
+    err = sqrt(err);
+    std::cout << "RMS isosurface distance: " << err << std::endl;
+    
+    std::cout << "Mesh vertices: " << poly.size_of_vertices() << ", "
+              << "facets: " << poly.size_of_facets() << std::endl;
+    
+    std::cout << "Writing to " << fname << "..." << std::endl;
+    std::ofstream ofs(fname);
+    CGAL::print_polyhedron_wavefront(ofs, poly);
+}

python_isosurfaces_2018_2019/mesh_implicit_sphere.cpp

@@ -0,0 +1,55 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/Implicit_mesh_domain_3.h>
+#include <CGAL/make_mesh_3.h>
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef K::FT FT;
+typedef K::Point_3 Point;
+typedef FT (Function)(const Point&);
+typedef CGAL::Implicit_mesh_domain_3<Function,K> Mesh_domain;
+
+#ifdef CGAL_CONCURRENT_MESH_3
+typedef CGAL::Parallel_tag Concurrency_tag;
+#else
+typedef CGAL::Sequential_tag Concurrency_tag;
+#endif
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+// Function
+FT sphere_function (const Point& p)
+{ return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-1; }
+
+int main()
+{
+  // Domain (Warning: Sphere_3 constructor uses squared radius !)
+  Mesh_domain domain(sphere_function,
+                     K::Sphere_3(CGAL::ORIGIN, 2.));
+  // Mesh criteria
+  Mesh_criteria criteria(facet_angle=30, facet_size=0.1, facet_distance=0.025,
+                         cell_radius_edge_ratio=2, cell_size=0.1);
+  
+  // Mesh generation
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+  
+  // Output
+  std::ofstream medit_file("out.mesh");
+  c3t3.output_to_medit(medit_file);
+  
+  return 0;
+}
+

python_isosurfaces_2018_2019/mesh_two_implicit_spheres_with_balls.cpp

@@ -0,0 +1,106 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/Labeled_mesh_domain_3.h>
+#include <CGAL/Mesh_domain_with_polyline_features_3.h>
+#include <CGAL/make_mesh_3.h>
+
+#include <CGAL/IO/print_wavefront.h>
+
+// Kernel
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+
+// Domain
+typedef K::FT FT;
+typedef K::Point_3 Point;
+typedef FT (Function)(const Point&);
+typedef CGAL::Mesh_domain_with_polyline_features_3<
+  CGAL::Labeled_mesh_domain_3<K> > Mesh_domain;
+
+// Polyline
+typedef std::vector<Point>        Polyline_3;
+typedef std::list<Polyline_3>       Polylines;
+
+#ifdef CGAL_CONCURRENT_MESH_3
+typedef CGAL::Parallel_tag Concurrency_tag;
+#else
+typedef CGAL::Sequential_tag Concurrency_tag;
+#endif
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
+
+typedef CGAL::Mesh_complex_3_in_triangulation_3<
+  Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+// Function
+FT sphere_function1 (const Point& p)
+{ return CGAL::squared_distance(p, Point(CGAL::ORIGIN))-2; }
+
+FT sphere_function2 (const Point& p)
+{ return CGAL::squared_distance(p, Point(2, 0, 0))-2; }
+
+FT sphere_function (const Point& p)
+{
+  if(sphere_function1(p) < 0 || sphere_function2(p) < 0)
+    return -1;
+  else
+    return 1;
+}
+
+#include <cmath>
+
+int main()
+{
+  // Domain (Warning: Sphere_3 constructor uses squared radius !)
+  Mesh_domain domain =
+    Mesh_domain::create_implicit_mesh_domain(sphere_function,
+                                             K::Sphere_3(Point(1, 0, 0), 6.));
+
+  // Mesh criteria
+  Mesh_criteria criteria(edge_size = 0.15,
+                         facet_angle = 25, facet_size = 0.15,
+                         cell_radius_edge_ratio = 2, cell_size = 0.15);
+
+  // Create edge that we want to preserve
+  Polylines polylines (1);
+  Polyline_3& polyline = polylines.front();
+
+  for(int i = 0; i < 360; ++i)
+  {
+    Point p (1, std::cos(i*CGAL_PI/180), std::sin(i*CGAL_PI/180));
+    polyline.push_back(p);
+  }
+  polyline.push_back(polyline.front()); // close the line
+
+  // Insert edge in domain
+  domain.add_features(polylines.begin(), polylines.end());
+
+  // Mesh generation without feature preservation
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
+                                      CGAL::parameters::no_features());
+
+  std::ofstream medit_file("out-no-protection.vtu");
+  c3t3.output_to_medit(medit_file);
+  medit_file.close();
+  c3t3.clear();
+
+  // Mesh generation with feature preservation
+  c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+
+  // Output
+  medit_file.open("out-with-protection.vtu");
+  c3t3.output_to_medit(medit_file);
+  medit_file.close();
+
+  return 0;
+}

python_isosurfaces_2018_2019/nix-shell.sh

@@ -0,0 +1,3 @@
+#!/bin/sh
+
+nix-shell -p python37Packages.numpy-stl

python_isosurfaces_2018_2019/scratch_ply.py

@@ -0,0 +1,21 @@
+#!/usr/bin/env python
+
+import math
+
+def f1(T, I, O, P, F):
+    return lambda x: (x, (T + I * math.sin(P + x*F)) / O, (I * math.cos(P + x*F)) / O)
+
+def f2(T, I, O, P, F):
+    return lambda x: (x, -(T - I * math.sin(P + x*F)) / O, (I * math.cos(P + x*F)) / O)
+
+f = f2(O=2.0, I=0.4, F=20, P=math.pi, T=0.3)
+print("ply")
+print("format ascii 1.0")
+r = range(-400, 400)
+print("element vertex %d" % (len(r)))
+print("property float32 x")
+print("property float32 y")
+print("property float32 z")
+print("end_header")
+for xi in r:
+    print("%f %f %f" % f(float(xi) / 200))

python_isosurfaces_2018_2019/shell.nix

@@ -0,0 +1,9 @@
+{ pkgs ? import <nixpkgs> {} }:
+
+let python_with_deps = pkgs.python3.withPackages
+      (ps: [ps.numpy ps.numpy-stl]);
+in pkgs.stdenv.mkDerivation rec {
+  name = "gfx_scratch";
+
+  buildInputs = with pkgs; [ python_with_deps ];
+}

python_isosurfaces_2018_2019/shell_cgal.nix

@@ -0,0 +1,10 @@
+{ pkgs ? import <nixpkgs> {} }:
+
+let stdenv = pkgs.stdenv;
+in stdenv.mkDerivation rec {
+  name = "cgal_scratch";
+
+  buildInputs = with pkgs; [ cgal boost gmp mpfr ];
+}
+
+# g++ -lCGAL -lmpfr -lgmp mesh_an_implicit_function.cpp -o mesh_an_implicit_function.o

python_isosurfaces_2018_2019/spiral_parametric.py

@@ -0,0 +1,109 @@
+#!/usr/bin/env python3
+
+import sys
+import numpy
+import stl.mesh
+
+# TODO:
+# - This is a very naive triangulation strategy.  It needs fixing - the
+# way it handles 'flatter' areas isn't optimal at all, even if the
+# sharper areas are much better than from CGAL or libfive.
+# - Generate just part of the mesh and then copy.  It is rotationally
+# symmetric, as well as translationally symmetric at its period.
+
+fname = "spiral_outer0.stl"
+freq = 20
+phase = 0
+scale = 1/16 # from libfive
+inner = 0.4 * scale
+outer = 2.0 * scale
+rad = 0.3 * scale
+
+angle = lambda z: freq*z + phase
+
+# z starting & ending point:
+z0 = -20*scale
+z1 = 20*scale
+# Number of z divisions:
+m = 1600
+# Number of circle points:
+n = 1000
+
+dz = (z1 - z0) / (m-1)
+
+data = numpy.zeros((m-1)*n*2 + 2*n, dtype=stl.mesh.Mesh.dtype)
+# Vertex count:
+# From z0 to z0+dz is n circle points joined with 2 triangles to next -> n*2
+# z0+dz to z0+dz*2 is likewise... up through (m-1) of these -> (m-1)*n*2
+# Two endcaps each have circle points & center point -> 2*n
+# Thus: (m-1)*n*2 + 2*n
+v = data["vectors"]
+print("Vertex count: {}".format(m*n*2 + 2*n))
+
+verts = numpy.zeros((n, 3), dtype=numpy.float32)
+
+# For every z cross-section...
+for z_idx in range(m):
+    #sys.stdout.write(".")
+    # z value:
+    z = z0 + dz*z_idx
+    # Angle of center point of circle (radians):
+    # (we don't actually need to normalize this)
+    rad = angle(z)
+    c,s = numpy.cos(rad), numpy.sin(rad)
+    # Center point of circle:
+    cx, cy = (inner + outer)*numpy.cos(rad), (inner + outer)*numpy.sin(rad)
+    # For every division of the circular cross-section...
+    if z_idx == 0:
+        # Begin with z0 endcap as a special case:
+        verts_last = numpy.zeros((n, 3), dtype=numpy.float32)
+        verts_last[:, 0] = cx
+        verts_last[:, 1] = cy
+        verts_last[:, 2] = z
+    else:
+        verts_last = verts
+    verts = numpy.zeros((n, 3), dtype=numpy.float32)
+    for ang_idx in range(n):
+        # Step around starting angle (the 'far' intersection of the
+        # line at angle 'rad' and this circle):
+        rad2 = rad + 2*ang_idx*numpy.pi/n
+        # ...and generate points on the circle:
+        xi = cx + outer*numpy.cos(rad2)
+        yi = cy + outer*numpy.sin(rad2)
+        verts[ang_idx, :] = [xi, yi, z]
+        #print("i={}, z={}, rad={}, cx={}, cy={}, rad2={}, xi={}, yi={}".format(i,z,rad,cx,cy, rad2, xi, yi))
+    if z_idx == 0:
+        for i in range(n):
+            v[i][0,:] = verts[(i + 1) % n,:]
+            v[i][1,:] = verts[i,:]
+            v[i][2,:] = verts_last[i,:]
+            #print("Write vertex {}".format(i))
+    else:
+        for i in range(n):
+            # Vertex index:
+            vi = z_idx*n*2 + i*2 - n
+            v[vi][0,:] = verts[(i + 1) % n,:]
+            v[vi][1,:] = verts[i,:]
+            v[vi][2,:] = verts_last[i,:]
+            #print("Write vertex {}".format(vi))
+            v[vi+1][0,:] = verts_last[(i + 1) % n,:]
+            v[vi+1][1,:] = verts[(i + 1) % n,:]
+            v[vi+1][2,:] = verts_last[i,:]
+            #print("Write vertex {} (2nd half)".format(vi+1))
+
+# then handle z1 endcap:
+for i in range(n):
+    # See vi definition above.  z_idx ends at m-1, i ends at n-1, and
+    # so evaluate vi+1 (final index it wrote), add 1 for the next, and
+    # then use 'i' to step one at a time:
+    vi = (m-1)*n*2 + (n-1)*2 - n + 2 + i
+    v[vi][0,:] = verts[i,:]
+    v[vi][1,:] = verts[(i + 1) % n,:]
+    v[vi][2,:] = [cx, cy, z]
+    # Note winding order (1 & 2 flipped from other endcap)
+    #print("Write vertex {} (endcap)".format(vi))
+
+print("Writing {}...".format(fname))
+mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
+mesh.save(fname)
+print("Done.")

python_isosurfaces_2018_2019/spiral_parametric2.py

@@ -0,0 +1,201 @@
+#!/usr/bin/env python3
+
+import sys
+import numpy
+import stl.mesh
+
+# TODO:
+# - Fix correction around high curvatures.  It has boundary issues
+# between the functions.
+# - Check every vertex point against the actual isosurface.
+# - Check rotation direction
+# - Fix phase, which only works if 0!
+
+fname = "spiral_inner0_one_period.stl"
+freq = 20
+phase = 0
+scale = 1/16 # from libfive
+inner = 0.4 * scale
+outer = 2.0 * scale
+rad = 0.3 * scale
+ext_phase = 0
+
+"""
+fname = "spiral_outer90_one_period.stl"
+freq = 10
+#phase = numpy.pi/2
+phase = 0
+scale = 1/16 # from libfive
+inner = 0.9 * scale
+outer = 2.0 * scale
+rad = 0.3 * scale
+ext_phase = numpy.pi/2
+"""
+
+def angle(z):
+    return freq*z + phase
+
+def max_z():
+    # This value is the largest |z| for which 'radical' >= 0
+    # (thus, for x_cross to have a valid solution)
+    return (numpy.arcsin(rad / inner) - phase) / freq
+
+def radical(z):
+    return rad*rad - inner*inner * (numpy.sin(angle(z)))**2
+
+def x_cross(z, sign):
+    # Single cross-section point in XZ for y=0.  Set sign for positive
+    # or negative solution.
+    n1 = numpy.sqrt(radical(z))
+    n2 = inner * numpy.cos(angle(z))
+    if sign > 0:
+        return (n2-n1) / outer
+    else:
+        return (n2+n1) / outer
+
+def curvature_cross(z, sign):
+    # Curvature at a given cross-section point.  This is fugly because
+    # it was produced from Maxima's optimized expression.
+    a1 = 1/outer
+    a2 = freq**2
+    a3 = phase + z*freq
+    a4 = numpy.cos(a3)
+    a5 = a4**2
+    a6 = numpy.sin(a3)
+    a7 = a6**2
+    a8 = inner**2
+    a9 = numpy.sqrt(rad**2 - a8*a7)
+    a10 = -a2*(inner**4)*a5*a7 / (a9**3)
+    a11 = 1 / a9
+    a12 = -a2*a8*a5*a11
+    a13 = a2*a8*a7*a11
+    a14 = 1/(outer**2)
+    a15 = -freq*a8*a4*a6*a11
+    if sign > 0:
+        return -a1*(a13+a12+a10+a2*inner*a4) / ((a14*(a15+freq*inner*a6)**2 + 1)**(3/2))
+    else:
+        return a1*(a13+a13+a10-a2*inner*a4) / ((a14*(a15-freq*inner*a6)**2 + 1)**(3/2))
+
+def cross_section_xz(eps):
+    # Generate points for a cross-section in XZ.  'eps' is the maximum
+    # distance in either axis.
+    verts = []
+    signs = [-1, 1]
+    z_start = [0, max_z()]
+    z_end = [max_z(), 0]
+    # Yes, this is clunky and numerical:
+    for sign, z0, z1 in zip(signs, z_start, z_end):
+        print("sign={} z0={} z1={}".format(sign, z0, z1))
+        z = z0
+        x = x_cross(z, sign)
+        while (sign*z) >= (sign*z1):
+            verts.append([x, 0, z])
+            x_last = x
+            dz = -sign*min(eps, abs(z - z1))
+            if abs(dz) < 1e-8:
+                break
+            x = x_cross(z + dz, sign)
+            #curvature = max(abs(curvature_cross(z, sign)), abs(curvature_cross(z + dz, sign)))
+            curvature = abs(curvature_cross((z + dz)/2, sign))
+            dx = (x - x_last) * curvature
+            print("start x={} dx={} z={} dz={} curvature={}".format(x, dx, z, dz, curvature))
+            while abs(dx) > eps:
+                dz *= 0.8
+                x = x_cross(z + dz, sign)
+                curvature = abs(curvature_cross((z + dz)/2, sign))
+                #curvature = max(abs(curvature_cross(z, sign)), abs(curvature_cross(z + dz, sign)))
+                dx = (x - x_last) * curvature
+                print("iter x={} dx={} z={} dz={} curvature={}".format(x, dx, z, dz, curvature))
+            z = z + dz
+            print("finish x={} z={} curvature={}".format(x, z, curvature))
+    n = len(verts)
+    data = numpy.zeros((n*2, 3))
+    data[:n, :] = verts
+    data[n:, :] = verts[::-1]
+    data[n:, 2] = -data[n:, 2]
+    return data
+
+def turn(points, dz):
+    # Note one full revolution is dz = 2*pi/freq
+    # How far to turn in radians (determined by dz):
+    rad = angle(dz)
+    c, s = numpy.cos(rad), numpy.sin(rad)
+    mtx = numpy.array([
+        [ c,  s,  0],
+        [-s,  c,  0],
+        [ 0,  0,  1],
+    ])
+    return points.dot(mtx) + [0, 0, dz]
+
+def screw_360(eps, dz):
+    #z0 = -10 * 2*numpy.pi/freq / 2
+    z0 = -5 * 2*numpy.pi/freq / 2
+    z1 = z0 + 2*numpy.pi/freq
+    #z1 = 5 * 2*numpy.pi/freq / 2
+    #z0 = 0
+    #z1 = 2*numpy.pi/freq
+    init_xsec = cross_section_xz(eps)
+    num_xsec_steps = init_xsec.shape[0]
+    zs = numpy.arange(z0, z1, dz)
+    num_screw_steps = len(zs)
+    vecs = num_xsec_steps * num_screw_steps * 2 + 2*num_xsec_steps
+    print("Generating {} vertices...".format(vecs))
+    data = numpy.zeros(vecs, dtype=stl.mesh.Mesh.dtype)
+    v = data["vectors"]
+    # First endcap:
+    center = init_xsec.mean(0)
+    for i in range(num_xsec_steps):
+        v[i][0,:] = init_xsec[(i + 1) % num_xsec_steps,:]
+        v[i][1,:] = init_xsec[i,:]
+        v[i][2,:] = center
+    # Body:
+    verts = init_xsec
+    for i,z in enumerate(zs):
+        verts_last = verts
+        verts = turn(init_xsec, z-z0)
+        if i > 0:
+            for j in range(num_xsec_steps):
+                # Vertex index:
+                vi = num_xsec_steps + (i-1)*num_xsec_steps*2 + j*2
+                v[vi][0,:] = verts[(j + 1) % num_xsec_steps,:]
+                v[vi][1,:] = verts[j,:]
+                v[vi][2,:] = verts_last[j,:]
+                #print("Write vertex {}".format(vi))
+                v[vi+1][0,:] = verts_last[(j + 1) % num_xsec_steps,:]
+                v[vi+1][1,:] = verts[(j + 1) % num_xsec_steps,:]
+                v[vi+1][2,:] = verts_last[j,:]
+                #print("Write vertex {} (2nd half)".format(vi+1))
+    # Second endcap:
+    center = verts.mean(0)
+    for i in range(num_xsec_steps):
+        vi = num_xsec_steps * num_screw_steps * 2 + num_xsec_steps + i
+        v[vi][0,:] = center
+        v[vi][1,:] = verts[i,:]
+        v[vi][2,:] = verts[(i + 1) % num_xsec_steps,:]
+    v[:, :, 2] += z0 + ext_phase / freq
+    v[:, :, :] /= scale
+    mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
+    print("Beginning z: {}".format(z0/scale))
+    print("Ending z: {}".format(z1/scale))
+    print("Period: {}".format((z1-z0)/scale))
+    return mesh
+
+#print("Writing {}...".format(fname))
+#mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
+#mesh.save(fname)
+#print("Done.")
+
+# What's up with this?  Note the jump from z=0.0424 to z=0.037.
+"""
+finish x=0.13228756555322954 z=0.042403103949074046 curvature=2.451108140319032
+sign=1 z0=0.042403103949074046 z1=0
+__main__:75: RuntimeWarning: invalid value encountered in double_scalars
+start x=0.0834189730812818 dx=nan z=0.042403103949074046 dz=-0.005 curvature=nan
+finish x=0.0834189730812818 z=0.03740310394907405 curvature=nan
+"""
+# Is it because curvature is undefined there - thus the starting step
+# size of 0.005 is fine?
+
+m = screw_360(0.01, 0.001)
+print("Writing {}...".format(fname))
+m.save(fname)

python_isosurfaces_2018_2019/spiral_parametric3.py

@@ -0,0 +1,234 @@
+#!/usr/bin/env python3
+
+import sys
+import numpy
+import stl.mesh
+
+# TODO:
+# - This still has some strange errors around high curvature.
+# They are plainly visible in the cross-section.
+# (errr... which errors were those? I can't see in the render)
+# - Check rotation direction
+# - Fix phase, which only works if 0 (due to how I work with y)
+# Things don't seem to line up right.
+# - Why is there still a gap when using Array modifier?
+# Check beginning and ending vertices maybe
+# - Organize this so that it generates both meshes when run
+# - Use SymPy instead of all this hard-coded stuff?
+
+# This is all rather tightly-coupled.  Almost everything is specific
+# to the isosurface I was trying to generate. walk_curve may be able
+# to generalize to some other shapes.
+class ExplicitSurfaceThing(object):
+    def __init__(self, freq, phase, scale, inner, outer, rad, ext_phase):
+        self.freq = freq
+        self.phase = phase
+        self.scale = scale
+        self.inner = inner
+        self.outer = outer
+        self.rad = rad
+        self.ext_phase = ext_phase
+        
+    def angle(self, z):
+        return self.freq*z + self.phase
+
+    def max_z(self):
+        # This value is the largest |z| for which 'radical' >= 0
+        # (thus, for x_cross to have a valid solution)
+        return (numpy.arcsin(self.rad / self.inner) - self.phase) / self.freq
+
+    def radical(self, z):
+        return self.rad*self.rad - self.inner*self.inner * (numpy.sin(self.angle(z)))**2
+
+    # Implicit curve function
+    def F(self, x, z):
+        return (self.outer*x - self.inner*numpy.cos(self.angle(z)))**2 + (self.inner*numpy.sin(self.angle(z)))**2 - self.rad**2
+
+    # Partial 1st derivatives of F:
+    def F_x(self, x, z):
+        return 2 * self.outer * self.outer * x - 2 * self.inner * self.outer * numpy.cos(self.angle(z))
+
+    def F_z(self, x, z):
+        return 2 * self.freq * self.inner * self.outer * numpy.sin(self.angle(z))
+
+    # Curvature:
+    def K(self, x, z):
+        a1 = self.outer**2
+        a2 = x**2
+        a3 = self.freq*z + self.phase
+        a4 = numpy.cos(a3)
+        a5 = self.inner**2
+        a6 = a4**2
+        a7 = self.freq**2
+        a8 = numpy.sin(a3)**2
+        a9 = self.outer**3
+        a10 = self.inner**3
+        return -((2*a7*a10*self.outer*x*a4 + 2*a7*a5*a1*a2)*a8 + (2*a7*self.inner*a9*x**3 + 2*a7*a10*self.outer*x)*a4 - 4*a7*a5*a1*a2) / ((a7*a5*a2*a8 + a5*a6 - 2*self.inner*self.outer*x*a4 + a1*a2) * numpy.sqrt(4*a7*a5*a1*a2*a8 + 4*a5*a1*a6 - 8*self.inner*a9*x*a4 + 4*a2*self.outer**4))
+
+    def walk_curve(self, x0, z0, eps, thresh = 1e-3, gd_thresh = 1e-7):
+        x, z = x0, z0
+        eps2 = eps*eps
+        verts = []
+        iters = 0
+        # Until we return to the same point at which we started...
+        while True:
+            iters += 1
+            verts.append([x, 0, z])
+            # ...walk around the curve by stepping perpendicular to the
+            # gradient by 'eps'.  So, first find the gradient:
+            dx = self.F_x(x, z)
+            dz = self.F_z(x, z)
+            # Normalize it:
+            f = 1/numpy.sqrt(dx*dx + dz*dz)
+            nx, nz = dx*f, dz*f
+            # Find curvature at this point because it tells us a little
+            # about how far we can safely move:
+            K_val = abs(self.K(x, z))
+            eps_corr = 2 * numpy.sqrt(2*eps/K_val - eps*eps)
+            # Scale by 'eps' and use (-dz, dx) as perpendicular:
+            px, pz = -nz*eps_corr, nx*eps_corr
+            # Walk in that direction:
+            x += px
+            z += pz
+            # Moving in that direction is only good locally, and we may
+            # have deviated off the curve slightly.  The implicit function
+            # tells us (sort of) how far away we are, and the gradient
+            # tells us how to minimize that:
+            #print("W: x={} z={} dx={} dz={} px={} pz={} K={} eps_corr={}".format(
+            #    x, z, dx, dz, px, pz, K_val, eps_corr))
+            F_val = self.F(x, z)
+            count = 0
+            while abs(F_val) > gd_thresh:
+                count += 1
+                dx = self.F_x(x, z)
+                dz = self.F_z(x, z)
+                f = 1/numpy.sqrt(dx*dx + dz*dz)
+                nx, nz = dx*f, dz*f
+                # If F is negative, we want to increase it (thus, follow
+                # gradient).  If F is positive, we want to decrease it
+                # (thus, opposite of gradient).
+                F_val = self.F(x, z)
+                x += -F_val*nx
+                z += -F_val*nz
+                # Yes, this is inefficient gradient-descent...
+            diff = numpy.sqrt((x-x0)**2 + (z-z0)**2)
+            #print("{} gradient-descent iters. diff = {}".format(count, diff))
+            if iters > 100 and diff < thresh:
+                #print("diff < eps, quitting")
+                #verts.append([x, 0, z])
+                break
+        data = numpy.array(verts)
+        return data
+
+    def x_cross(self, z, sign):
+        # Single cross-section point in XZ for y=0.  Set sign for positive
+        # or negative solution.
+        n1 = numpy.sqrt(self.radical(z))
+        n2 = self.inner * numpy.cos(self.angle(z))
+        if sign > 0:
+            return (n2-n1) / self.outer
+        else:
+            return (n2+n1) / self.outer
+
+    def turn(self, points, dz):
+        # Note one full revolution is dz = 2*pi/freq
+        # How far to turn in radians (determined by dz):
+        rad = self.angle(dz)
+        c, s = numpy.cos(rad), numpy.sin(rad)
+        mtx = numpy.array([
+            [ c,  s,  0],
+            [-s,  c,  0],
+            [ 0,  0,  1],
+        ])
+        return points.dot(mtx) + [0, 0, dz]
+
+    def screw_360(self, z0_period_start, x_init, z_init, eps, dz, thresh, endcaps=False):
+        #z0 = -10 * 2*numpy.pi/freq / 2
+        z0 = z0_period_start * 2*numpy.pi/self.freq / 2
+        z1 = z0 + 2*numpy.pi/self.freq
+        #z1 = 5 * 2*numpy.pi/freq / 2
+        #z0 = 0
+        #z1 = 2*numpy.pi/freq
+        init_xsec = self.walk_curve(x_init, z_init, eps, thresh)
+        num_xsec_steps = init_xsec.shape[0]
+        zs = numpy.append(numpy.arange(z0, z1, dz), z1)
+        num_screw_steps = len(zs)
+        vecs = num_xsec_steps * num_screw_steps * 2
+        offset = 0
+        if endcaps:
+            offset = num_xsec_steps
+            vecs += 2*num_xsec_steps
+        print("Generating {} vertices...".format(vecs))
+        data = numpy.zeros(vecs, dtype=stl.mesh.Mesh.dtype)
+        v = data["vectors"]
+        # First endcap:
+        if endcaps:
+            center = init_xsec.mean(0)
+            for i in range(num_xsec_steps):
+                v[i][0,:] = init_xsec[(i + 1) % num_xsec_steps,:]
+                v[i][1,:] = init_xsec[i,:]
+                v[i][2,:] = center
+        # Body:
+        verts = init_xsec
+        for i,z in enumerate(zs):
+            verts_last = verts
+            verts = self.turn(init_xsec, z-z0)
+            if i > 0:
+                for j in range(num_xsec_steps):
+                    # Vertex index:
+                    vi = offset + (i-1)*num_xsec_steps*2 + j*2
+                    v[vi][0,:] = verts[(j + 1) % num_xsec_steps,:]
+                    v[vi][1,:] = verts[j,:]
+                    v[vi][2,:] = verts_last[j,:]
+                    #print("Write vertex {}".format(vi))
+                    v[vi+1][0,:] = verts_last[(j + 1) % num_xsec_steps,:]
+                    v[vi+1][1,:] = verts[(j + 1) % num_xsec_steps,:]
+                    v[vi+1][2,:] = verts_last[j,:]
+                    #print("Write vertex {} (2nd half)".format(vi+1))
+        # Second endcap:
+        if endcaps:
+            center = verts.mean(0)
+            for i in range(num_xsec_steps):
+                vi = num_xsec_steps * num_screw_steps * 2 + num_xsec_steps + i
+                v[vi][0,:] = center
+                v[vi][1,:] = verts[i,:]
+                v[vi][2,:] = verts[(i + 1) % num_xsec_steps,:]
+        v[:, :, 2] += z0 + self.ext_phase / self.freq
+        v[:, :, :] /= self.scale
+        mesh = stl.mesh.Mesh(data, remove_empty_areas=False)
+        print("Beginning z: {}".format(z0/self.scale))
+        print("Ending z: {}".format(z1/self.scale))
+        print("Period: {}".format((z1-z0)/self.scale))
+        return mesh
+
+surf1 = ExplicitSurfaceThing(
+    freq = 20,
+    phase = 0,
+    scale = 1/16, # from libfive
+    inner = 0.4 * 1/16,
+    outer = 2.0 * 1/16,
+    rad = 0.3 * 1/16,
+    ext_phase = 0)
+
+z_init = 0
+x_init = surf1.x_cross(z_init, 1)
+mesh1 = surf1.screw_360(-10, x_init, z_init, 0.000002, 0.001, 5e-4)
+fname = "spiral_inner0_one_period.stl"
+print("Writing {}...".format(fname))
+mesh1.save(fname)
+
+surf2 = ExplicitSurfaceThing(
+    freq = 10,
+    phase = 0,
+    scale = 1/16, # from libfive
+    inner = 0.9 * 1/16,
+    outer = 2.0 * 1/16,
+    rad = 0.3 * 1/16,
+    ext_phase = numpy.pi/2)
+
+z_init = 0
+x_init = surf2.x_cross(z_init, 1)
+mesh2 = surf2.screw_360(-5, x_init, z_init, 0.000002, 0.001, 5e-4)
+fname = "spiral_outer90_one_period.stl"
+print("Writing {}...".format(fname))
+mesh2.save(fname)