.gitignore

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

README.md

@@ -3,14 +3,6 @@
 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
@@ -32,21 +24,10 @@ 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 - but
+  later, particularly for the integration with Blender.
-  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

@@ -1,22 +0,0 @@
-* 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

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

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

@@ -1,232 +0,0 @@
-#!/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 everything,
+- Winding order.  It is consistent through seemingly
-  except for reflection and `close_boundary_simple`.  (When there are
+  everything, except for reflection and close_boundary_simple.
-  two parallel boundaries joined with something like
+  (When there are two parallel boundaries joined with something like
-  `join_boundary_simple`, traversing these boundaries in their actual
+  join_boundary_simple, traversing these boundaries in their actual order
-  order to generate triangles - like in `close_boundary_simple` - will
+  to generate triangles - like in close_boundary_simple - will produce
-  produce opposite winding order on each. Imagine a transparent clock:
+  opposite winding order on each. Imagine a transparent clock: seen from the
-  seen from the front, it moves clockwise, but seen from the back, it
+  front, it moves clockwise, but seen from the back, it moves
-  moves counter-clockwise.)
+  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

@@ -1,10 +0,0 @@
-#!/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

@@ -1,173 +0,0 @@
-#+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

@@ -1,171 +0,0 @@
-// 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

@@ -1,55 +0,0 @@
-#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

@@ -1,106 +0,0 @@
-#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

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

python_isosurfaces_2018_2019/scratch_ply.py

@@ -1,21 +0,0 @@
-#!/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

@@ -1,9 +0,0 @@
-{ 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

@@ -1,10 +0,0 @@
-{ 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

@@ -1,109 +0,0 @@
-#!/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

@@ -1,201 +0,0 @@
-#!/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

@@ -1,234 +0,0 @@
-#!/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)