Computer Algebra and Graphics
In part 1, we walked through building a computation graph of implicit surfaces in Tensorflow. In this post we’re going to talk about incorporating the computer algebra system Sympy into the picture by generating an algebraic abstract syntax tree that we compile at render time.
Algebraic Abstract Syntax Trees
Abstract syntax trees are a data structure used in compilers to represent the intended computation written by the programmer as a hierarchy of operations. The compiler can then use heuristic or mathematical rewrite rules on that tree to produce a more efficient computation that computes the same result, and then can traverse the tree and translate each node of the tree into an executable piece of (usually lower-level) code. We’re going to apply this basic pipeline to algebraic expressions and see how far we can get!
The first thing to convince ourselves of is that algebraic expressions can be thought of as compose-able trees. As motivation, let’s take one of the more useful transformations of space, that we use to do replication of unit cells:
def saw_wave(amplitude, frequency):
def replicated(x, y, z):
fx, fy, fz = [
amplitude * (coord / frequency) - floor((coord / frequency) + 1/2
for coord in [x, y, z]
]
return (fx, fy, fz)
return replicatedHere, we map repeated cells to a unit cube of length amplitude. To see what this looks like in practice, here’s an example I wrote in GLSL. Note that the size of the replicated unit cube is specified by the frequency input. But what if we wanted to have the replicated unit cell size to vary as a function of x,y,z? Practically, this looks like the pictures in this paper. In category theory, this would be called an Operad, which is roughly a fancy term for “an operation that can be composed”.
We want to treat Operads as nodes in a tree, and when it comes time to evaluate the computation, we will call the top Operad and it will traverse the tree, substituting the values from the input dictionary into the underlying expressions. We’ll use sympy.Expressions as the data structure for the algebraic expression of each Operad.
With this piece, we can have Operads that take other Operads as input, and then recursively substituted the symbols in each expression. Returning to our motivating example, we can have the replicated cell size depend on expressions like sqrt(x**2+y**2), the distance from the z axis, which can be useful if you’re trying to model a microstructure that gets thicker as it approaches the medial axis of your part.
The next part is building up these Operad trees from smaller pieces.
Mutable Closures
In Python, there’s a language-level concept of a closure, which is what we’ve been using so far to contain the information required to render our primitives. Think of closures as having dictionaries with keys that are variable names and values that are the variable values at definition time. For instance:
def func(r):
def expr(x):
return x - r
return exprfunc(2).__closure__[0].cell_contents # <- 2
func(2).__closure__[0].cell_contents = 3 # <- error: not writeable
Closures do have their downsides though, because it’s impossible to modify the closure scope once it’s been initialized. We want a closure-like data structure that we can substitute with values at will, and that is compatible with Sympy. In order to do that, we’ll decorate our primitives with a class that calls the decorated function with sympy.Symbols, and then remembers the argument values the user supplied at call time for later substitution.
Why not just use keyword args to solve the closure problem? On the one hand we could use it to symbolize function inputs:
from sympy import Symbol
def func(r=Symbol('r')):
def expr(X=Symbol('x')):
return X - r
return exprfunc()() # <- x - r
func(2)() # <- x - 2
However, I couldn’t figure out a way to both symbolize the inner expression and retain the information about what the function was originally called with without a wrapper class.
Compiling Sympy Expressions to Tensorflow Graphs
This part is much simpler than it sounds, mostly because of the work done by the Sympy team to provide a blueprint for compilation. When I picked Sympy up, it already supported compilation to Javascript, C++, Fortran, and even Numpy for numerical evaluation. Building on what was already there for Numpy, I filled out a code printer for Tensorflow that simply took each node of the expression graph and translated it to a Tensorflow-compatible code string. Then that string gets evaled into a lambda expression, and we call that lambda with tf.Tensors.
Symbolic Algebraic Geometry
Now that we have built ASTs for each function in our library, we can use Sympy’s algebraic manipulation capabilities to build higher-order expressions. The use case of interest is rendering sweep surfaces, and we’ll start with a special case of sweeps called channel surfaces.
Imagine a sphere sphere(x,y,z,t) := radius**2 — (x — p(t)[0])**2 — (y — p(t)[1])**2 — (z - p(t)[2])**2 whose center is moving along a space curve p(t) := (px, py, pz), creating a channel through some imaginary material. The channel surface related to that curve is the resulting channel through the material. Mathematically, this surface is specified by the set of all points (x,y,z,t) such that sphere(x,y,z,t) = 0 and d(sphere)/dt (x,y,z,t) = 0.
The intuition here is that the second equation is satisfied at a time t by all points (x,y,z)that lie on the plane orthogonal to the tangent of the curve at time t. So if you take the intersection of the sphere at time t with that plane, you get the outer rind of the sphere, and when you union these rinds for every point in time, you get the full channel surface.
Now we have to solve two equations to find the set of points on the surface. In algebraic geometry, this is an example of a variety; succinctly, a set of simultaneous solutions of a system of polynomial equations. We can use symbolic algebra, specifically the theory of resultants, to eliminate the variable t from our system of equations and produce a polynomial in just x,y,z.
This works by first viewing our two polynomials as members of the polynomial ring with coefficients in R[x,y,z] , so the coefficients of the polynomials are polynomials themselves in x,y,z. Then, we build a matrix from multiplication of the polynomials by1, t, t^2, ..., to produce a square matrix with linearly independent columns, called the Sylvester Matrix. Finally, we take the determinant of this matrix, which is the resultant of the two polynomials. The resultant is in general a polynomial in x,y,z, since the matrix entries are all polynomials in x,y,z, and is zero whenever the two original polynomials are simultaneously zero, because the matrix is designed to be invertible whenever one of the polynomials is not zero. Note that Sympy uses a different implementation of the resultant that is more efficient and general, but the basic idea is the same.
I relied on ideas from this paper and this book heavily for this part of the project. In practice, I use these surfaces as pipes that interpolate between points, and I generate the directrix via 3 dimensional Hermite interpolation.
Next Time: Nonlinear Optimization to Render Arbitrary Sweeps
This approach to computing sweeps only works for a very small set of primitives that are polynomials, and doesn’t guarantee anything about the magnitude of the scalar field, which can cause numerical issues when you try to combine channel surfaces with regular distance fields. In order to render arbitrary sweeps and improve on the issues mentioned above by rendering true signed distance channel surfaces, we need to reformulate the problem to allow us to solve a nonlinear function of t at every point in our spatial grid.
Here’s a picture of what we’ll be able to draw:
