Add note on parallel_transport

.gitignore

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

README.md

@@ -3,6 +3,14 @@
 This is repo has a few projects that are related in terms of
 high-level goal, but almost completely unrelated in their descent.
 
+- `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
@@ -24,10 +32,21 @@ high-level goal, but almost completely unrelated in their descent.
   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.
+  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.
+
+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.
+
+## Projects not covered here
+
+- curl-noise work (both in Clojure and in Python/vispy)
+- parallel transport
+- prosha, of course

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?)

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

@@ -8,20 +8,20 @@
 - 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
+- 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
+- Winding order.  It is consistent through seemingly everything,
-  everything, except for reflection and close_boundary_simple.
+  except for reflection and `close_boundary_simple`.  (When there are
-  (When there are two parallel boundaries joined with something like
+  two parallel boundaries joined with something like
-  join_boundary_simple, traversing these boundaries in their actual order
+  `join_boundary_simple`, traversing these boundaries in their actual
-  to generate triangles - like in close_boundary_simple - will produce
+  order to generate triangles - like in `close_boundary_simple` - will
-  opposite winding order on each. Imagine a transparent clock: seen from the
+  produce opposite winding order on each. Imagine a transparent clock:
-  front, it moves clockwise, but seen from the back, it moves
+  seen from the front, it moves clockwise, but seen from the back, it
-  counter-clockwise.)
+  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,
@@ -39,7 +39,7 @@
     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
+- `/mnt/dev/graphics_misc/isosurfaces_2018_2019` - perhaps include my
   spiral isosurface stuff from here
 
 ## Abstractions

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)