A comparison of some deconvolution priors

Continuing with experiments based on ADMM I implemented and experimented with a few different deconvolution priors.  I chose several $\ell_2$ and $\ell_1$ priors based on sparse gradients, sparse curvature and simply the norm of the solution vector as well as an $\ell_{0.5}$ HyperLaplacian prior.

For all results the ground-truth signal is shown in black, the blurred and corrupted signal in red and the deblurred signal in blue.  As a ground-truth signal I'm using a simple piecewise linear signal with a slanted region, some constant regions and discontinuous jumps.  These features help to show the issues with different regularizers.  The input signal is 128 samples long, blurred by a Gaussian filter with $\sigma=4$ and then corrupted by Gaussian noise with $\sigma=2.5%$ and $\mu=0$. 

For the $\ell_2$ results I'm solving for the optimal solution directly in the Fourier domain, while for the $\ell_1$ and $\ell_0.5$ I'm using ADMM, solving the data-subproblem in the Fourier domain and the splitting variable subproblem with shrinkage/proximal operators.  For all problems I've played with the prior weights a bit to try to get good results, but don't claim they are optimal.

$\ell_2$ norm of solution vector

First up is the $\ell_2$ norm of the solution.  This does sharpen up the signal somewhat, but has bad ringing artifacts.  These are expected for the $\ell_2$ solution, particularly given the high noise levels in the corrupted signal.  The relatively large blur also leaves lots of freedom for low-frequency ringing.

$\ell_2$ norm of solution first derivative

Results for the $\ell_2$-norm of the signal derivatives do a bit better, particularly for the slanted segment, but still show low-frequency ringing.  The same is true for the $\ell_2$-norm of the second derivative, although the ringing is a bit more pronounced in the flat regions:

$\ell_2$ norm of solution second derivative

Moving on to the $\ell_1$ priors, the total-variation prior is simply the $1$-norm of the signal first derivatives. It favours piecewise constant solutions which do well on the large flat regions but introduce staircase artifacts for the slanted region:

$\ell_1$ norm of solution first derivative

The total variation prior produces a much sharper signal than the $\ell_2$ priors, however for the slanted region it also introduces spurious features in the slanted region that are hard to distinguish from (correct) features elsewhere in the signal.  This occurs because the assumption of the total variation prior is that the input is piecewise constant, which is not true for this signal.  A better assumption is that the solution is piecewise linear, leading to the sparse-curvature prior.  The sparse-curvature assumes that the second derivative of the signal is non-zero in only a few locations:

$\ell_1$ norm of solution second derivative

The sparse-curvature prior does much better in the slanted region than the total variation prior and still improves the signal, but the edges are not as sharp at discontinuities.  Often people combine these priors to try to achieve a compromise between the two behaviours.

The final prior is the Hyper-Laplacian prior, which is simply the $\ell_{0.5}$-norm of the signal first derivatives:

$\ell_{0.5}$ norm of solution first derivative

The Hyper-Laplacian prior is non-convex, making it difficult to optimize for.  Here it appears that the optimization got stuck: some edges are extremely sharp, however the trailing edge of the high-amplitude constant region is lost entirely, possible because of the slanted region adjacent to it.  I played with the parameters quite a bit but did not get a result better than this. 

Matlab example code generating these results is available at my GitHub repository:

https://github.com/jamesgregson/Matlab-deblurring-example

Matlab/MEX utility for loading VTK Rectilinear Grid files (.vtr files)

The title says it all, I've written some basic code for loading VTK rectilinear grid files into Matlab.  The code supports uniformly spaced meshes in up to four dimensions for both point and cell data.

You can download the code from my Github repository:
https://github.com/jamesgregson/matlab_vtr

Index transformation between bounding boxes and uniform grids

I end up rewriting this code all the time. It transforms from points in space to grid coordinates world_to_grid() and back grid_to_world(), such that p = grid_to_world( aabb, dim, world_to_grid( aabb, dim, p ) ).  It's simple, but bugs here can mess up lots of things in ways that are hard to detect.  No range checking is performed in order to make it easy to apply custom boundary conditions.

    template< typename real, typename index, typename real3 >
    inline real3 _world_to_grid( const real *aabb, const index *dim, const real3 &pos ){
        return real3( real(dim[0]-1)*(pos[0]-aabb[0])/(aabb[1]-aabb[0]), real(dim[1]-1)*(pos[1]-aabb[2])/(aabb[3]-aabb[2]), real(dim[2]-1)*(pos[2]-aabb[4])/(aabb[5]-aabb[4]) );
    }
    
    template< typename real, typename index, typename real3 >
    inline real3 _grid_to_world( const real *aabb, const index *dim, const real3 &pos ){
        return real3( aabb[0]+real(pos[0])*(aabb[1]-aabb[0])/real(dim[0]-1), aabb[2]+real(pos[1])*(aabb[3]-aabb[2])/real(dim[1]-1), aabb[4]+real(pos[2])*(aabb[5]-aabb[4])/real(dim[2]-1) );
    }  

The code assumes that the dim variable holds the number of points in each direction (e.g. int dim[] = { nx, ny, nz };) and that the first point is coincident with the minimum value of the axis-aligned bounding box with the last point coincident with the maximum value of the bounding box. The aabb parameter is the axis-aligned bounding box defining the grid in world space, e.g. int aabb[] = { xmin, xmax, ymin, ymax, zmin, zmax };

A follow-up to fluid simulation on non-uniform grids

In my last post, I discussed preliminary results for fluid simulation on non-uniform Cartesian grids.  In that post I showed some preliminary results, but there were some bugs that added disturbing artifacts.

I have fixed a number of bugs and now have a solver based on BFECC advection using min-max limited cubic interpolation for non-uniform and often highly anisotropic meshes for the velocity and pressure solve, with high-resolution uniform density fields for the density field.  The results are fairly impressive:

This image shows a volume rendering (in Paraview) of a simulation computed using the grid in the background.  Near the area of interest, the grid is uniform, but it grows very quickly (geometrically, with growth rate ~1.5) outside this region.  The velocity/pressure grid is 133x127x134, but covers nearly a 10x6x10 m cubic volume, with 2cm cells in the fine region.  The density field is 1x6x1 m with 1cm uniform resolution.

Being able to run different resolutions and gradings for the velocity and density fields is extremely helpful: fine fluid details distract from a poor velocity solution, and high-resolution densities help avoid diffusion in the physics.  The image above shows the density as resolved by the fluid grid. It is terrible.  However the density as resolved by the density grid is -way- better:

It's still not perfect, but given the cell size and anisotropy, I think it does extremely well.  Although there are definite artifacts, the payoff is in the memory usage and runtime. The whole simulation is 5 seconds in real time and takes approximately 22 seconds per output frame, meaning my 5 second simulation completes in under an hour.  These simulations used to take on the order of 4-5 hours.

The results look pretty good.  I think the grading is too steep to get really nice results, but it's an excellent proof-of-concept.

Fluid Simulation for Graphics on Non-Uniform Structured Grids

I've been playing around with fluid simulation on non-uniform structured grids.  They have some charms in that it is easy to have very large domains with isolated regions of interest; e.g. near the camera or salient fluid features.  One of the big advantages is that it makes far-field boundaries easy, you simply extend the mesh far away from the domain.

My solver pretty run-of-the-mill; only the pressure solver was updated in order to handle anisotropic cells.  I've used a finite-volume formulation to derive the pressure-correction equation for this case, but it suffices to say that only the weights used in forming the Poisson equation change.

Here is an example, a grid that is roughly 10x10x5 meters, with 2cm fine cells for roughly 100x100x100 cells:

The aspect ratios get quite high, 100:1 is not uncommon.   You can see a buoyant flow simulation that I'm running as the contour values.  The simulation takes about 20 seconds per frame, of which 7 seconds is writing the VTK output files that I use to post-process.  In spite of that, the visual detail in the region of interest is excellent (given the resolution):

I find that it helps to advect a uniform resolution density field that is roughly 2X the fine cell resolution.  Writing the code with general-purpose fields that can be point-sampled arbitrarily makes this a trivial feature to implement, simply allocate different grids for density and velocity.

Finally here is a rendering of the final simulation in Paraview, total simulation time ~35 minutes. There are still some bugs to track down, but the results are pretty promising:

Updated 4x4 transformation matching OpenGL

I've updated the code from my previous post: 4x4 transformation matching OpenGL to remove the requirement for a separate linear algebra library by introducing some 4x4 inverse code from Rodolphe Vaillant (originally at http://www.irit.fr/~Rodolphe.Vaillant/?e=7) along with some other minor improvements.  You can still use Eigen or GMM to perform matrix inverses, but it is no longer a requirement.

As before, the code re-implements the main OpenGL matrix functions and includes code for testing the output against the system OpenGL.  I use the code for a number of my imaging projects where it is desirable (but optional) to have an OpenGL preview.

The code is now also hosted at Github which should simply future updates:

https://github.com/jamesgregson/transformation

C++/Mex Image Deblurring using ADMM

I've posted some sample code on Github for performing image deblurring in Matlab using Mex.  I wrote it as a way to play around with the ADMM algorithm for sparse signal reconstruction, as described in Stephen Boyd's ADMM paper, as well as to get some experience using C++ code from Matlab.

The code uses matrix-free linear operators to evaluate the forward and reverse measurement (blurring) process, solving the data-term subproblem with gradient descent.  While you can do this much faster using FFTs for deconvolution, the code should generalize quite well to tomography with non-orthographic projections.  I've tried to comment it well so that I can figure out what I did if I ever need to pick it up again.

In the image above, the left plot is the ground-truth image and the right is blurred by a point-spread function and corrupted by Gaussian noise.  After running the demo, the following results are obtained:

The reconstructed image is shown on the right, while a slice through the center scanline shows the reconstruction (blue dots), ground-truth (solid black line) and input (red line).  From the results it's clear that the method does well at preserving discontinuities even for fairly large blurs and noise levels.

The demo is pretty quick to run and is helpful for choosing parameters since you can see how the behavior of the method changes with the regularizer weight and the penalty term weight.

You can download the code from https://github.com/jamesgregson/MexADMMDeblurDemo

FLIP fluid simulation with sample code

I've been continuing to play around with incompressible fluid simulation and have implemented the FLIP method (FLuid Implicit Particle). The FLIP algorithm represents the fluid concentration and velocity field using a large set of particles.  At every iteration these particles are transferred to an auxiliary grid using kernel density estimation which defines a grid-based velocity field which is made incompressible by pressure projection and then used to advect the particles. Finally the velocity of each particle is corrected using the difference between the projected velocity field and the initial velocity field obtained by kernel density estimation.

The primary benefit of this approach is that it has extremely low numerical viscosity, making simulating fluids such as smoke and water possible at sane grid resolutions.  The downside is considerably more complexity than a purely grid-based solver; you need to track particles, perform the mapping between particles and grid (and vice versa) and have some sensible scheme to reseed particles when either there are too few or too many in a region.

To play with the method I implemented it in C++, largely following Robert Bridson's Fluid Simulation for Computer Graphics.  It is a canned example demonstrating variable density flow under the influence of gravity using an body force proportional to the fluid 'concentration', a sort of poor-man's Boussinesq approximation.  The simulation below was performed on a 100x100 auxiliary grid with 16 particles per grid-cell. This is more than the 2x2 typically used, but helps get around reseeding.  As this is not free-surface flow, there are particles everywhere which allows causes the smoke to form nice instabilities and finally turbulent mixing (although I did not run it particularly far).

In order to have an unmangled version of the code available to myself, I am posting it here.  It is uncommented and unstructured, but (hopefully) fairly understandable once you understand the basic algorithm.  The code is bundled with the GMM++ sparse matrix library for the pressure solves and required VTK and CMake for output and as a build-system respectively.

You can download the code here: FLIP2.zip

Latex formulas as images using Python

I find it endlessly frustrating to incorporate math into vector-graphics documents and Powerpoint presentations. While Powerpoint does allow importing equations from the equation editor, they usually look terrible.  However there is an online Latex equation editor run by codecogs.com that allows you to generate images of expressions.  This works great, but it's inconvenient to save and regenerate images after using it.

To that end, I've written a python function that will generate the URL for the Latex image-generation script, get the image data using HTTP and write it to disk as a file of your choosing.  It uses the python requests module (available via pip) and ImageMagick:

import os, requests

def formula_as_file( formula, file, negate=False ):
    tfile = file
    if negate:
        tfile = 'tmp.png'
    r = requests.get( 'http://latex.codecogs.com/png.latex?\dpi{300} \huge %s' % formula )
    f = open( tfile, 'wb' )
    f.write( r.content )
    f.close()
    if negate:
        os.system( 'convert tmp.png -channel RGB -negate -colorspace rgb %s' %file )

The code is appallingly simple for how much time it has saved me; it allows me to regenerate formulas for conference posters by simply running a python script and since most vector-graphics packages simply link to bitmap images allows the resulting layouts to update automatically.  Even better is that CodeCogs image generation script produces transparent PNGs so you can drop them on top of lightly colored backgrounds.  By setting the negate flag you will get a white on transparent PNG that can be used over dark backgrounds.

Here are some examples:

formula_as_file( r'\Gamma_{Levin}(x) = \| \nabla p(x) \|_2^{0.8} + \sum_i |\frac{\partial^2 p(x)}{\partial x_i^2}|^{0.8}', 'reg_levin.png', True )

formula_as_file( r'\Gamma_{MTV}(x) = \| \sum_{i=1}^3 \left( \| \nabla p^{(i)}(x)\|_2^2 \right) \|_2', 'reg_MTV.png' )

Finally no more ugly power-point presentations!

Efficiently Sample from Normal (Gaussian) Distribution

The following code snippet samples from a standard normal distribution using the Box-Muller transform, which avoids nastiness like cropping to a fixed number of standard-deviations or rejection sampling.  I have shamelessly cribbed this from Wikipedia (the previous link), but use it all the time, so I'm putting it up as a snippet.

void utils_sample_normal( float &x, float &y ){
    float u=drand48(), v=drand48();
    float lnu = sqrt( -2.0*log(u) );
    x = lnu*cos(2.0*M_PI*v);
    y = lnu*sin(2.0*M_PI*v);
}

Matlab Signal Deblurring & Denoising Example

To date my research has been largely focused on inverse problem such as tomography or image deblurring.  These problems are often highly under-determined and so must include strong priors to obtain good solutions and finding efficient solvers for these priors is challenging.  A labmate recently pointed me towards the ADMM method, which splits the full problem into coupled sub-problems.  It has a number of advantages:

• Flexibility - ADMM handles a number of non-trivial problems in a common framework
• Efficiency - Often the subproblems have efficient, often embarassingly-parallel, solvers

An excellent introduction to the method can be found in the paper: Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers by Boyd et al.. It's a pretty gentle and practical introduction, with numerous example problems.

I've quickly implemented ADMM for combined deblurring and denoising of 1D input signals using the total-variation regularization in a Generalized-Lasso problem definition.  Unlike the Boyd paper, I've chosen to use Landweber iterations to solve the data subproblem as these are commonly used in large-scale deblurring and tomography.  My well-commented sample implementation (see the end of this post) allows all the method parameters to be tweaked to see the effect of using inexact solves, different penalty parameters and different regularization weights.

As an example, the image below shows a 100 sample square-wave signal, blurred by a Gaussian with standard deviation of 3 and corrupted by Gaussian noise with a sigma of 5%.  Only two inner Landweber iterations were used for each outer iteration.  The red line shows the original blurred and noisy input, while the blue line shows the reconstructed result after 100 iterations.

Overall I find that the method works as advertised. It's fairly easy to implement, either in Matlab or in C/C++, can handle 1-norms of general linear regularizers and converges quickly.  It also allows matrix-free solvers to be used for the data-subproblem, while the solve for the prior is embarrassingly-parallel, so it should scale quite well too.

You can get the implementation that produced this plot below. Note that you may get slightly different results since the noise is generated randomly.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                                                     %%%
%%% Sample signal denoising & deblurring using ADMM                 %%%
%%% code written by James Gregson (james.gregson@gmail.com), 2013       %%%
%%% Use the code for whatever you'd like                                %%%
%%%                                                                     %%%
%%% Demonstrates performing Total-Variation (TV) denoising and          %%%
%%% deblurring of a 1D input signal using the ADMM iterative scheme     %%%
%%% applied to a Generalized-Lasso problem.  More information on the    %%%
%%% approach can be found in [1].  The approach from [1] is modified    %%%
%%% slightly by using gradient-descent (Landweber iterations) to solve  %%%
%%% the first subproblem in place of a direct solver or iterative       %%%
%%% method such as conjugate gradient.  Matrix-free Landweber           %%%
%%% iterations are commonly used for large-scale linear inverse         %%%
%%% problems and this sample code allows the accuracy of the            %%%
%%% subproblem solves to be adjusted to see the effect on the final     %%%
%%% reconstructions.                                                    %%%
%%%                                                                     %%%
%%% [1] Boyd et al., Distributed Optimization and Statistical Learning  %%%
%%%     via the Alternating Direction Method of Multipliers, 2010       %%%
%%%                                                                     %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Problem Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N           = 100;     % Signal sample count
lambda      = 0.1;     % Total-variation weight
sigma       = 3.0;     % PSF sigma for generating blurred input
noise       = 0.05;    % Gaussian-noise sigma

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ADMM Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rho         = 1.0;     % ADMM constraint weight 0 < rho <= 2
outer_iters = 100;     % Number of ADMM iterations
inner_iters = 2;       % Number of Landweber steps per ADMM iteration
relax       = 1.9;     % Under-relaxation factor for Landweber steps [0,2]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Construct System PSF & Image Formation Model %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

off = -ceil(sigma*3):ceil(sigma*3);  % psf pixel offsets
psf = exp( -off.^2/sigma^2 );        % psf values for offsets
psf = psf/sum(psf);                  % normalize to 1.0

% generate psf matrix by setting diagonals based on 1D psf above
M = zeros( N, N );                   
for i=1:numel(psf),
    M = M + diag( psf(i)*ones(N-abs(off(i)),1), off(i) );
end
M = sparse( M );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Generate noisy and blurred synthetic input %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

input = zeros( N, 1 );                  % start with zeros
input(floor(N/4):ceil(3*N/4),:) = 1.0;  % make a center region at 1.0
blur = M * input + noise*(randn(N,1));  % blur the input by the psf


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Construct the difference matrix that computes image gradients %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A = sparse( -diag( ones(N,1), 0 ) + diag( ones(N-1,1), 1 ) );
A(N,N-1)=1; % use a backward difference for the final point

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup ADMM and perform iterations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute the eigenvalues of the first subproblem
% system matrix and use the largest to define a
% Landweber step size, for large system a a power-
% iteration should be used instead.
tmp = eigs(M'*M + rho*A'*A);
step = relax/tmp(1);

% initialize the ADMM variables
x = blur;             % intrinsic (sharp) signal
z = zeros( N, 1 );    % splitting variable
u = zeros( N, 1 );    % scaled Lagrange multipliers

% define an anonymous shrinkage operator to implement
% the second sub-problem solve
shrink = @(kappa,x) max( abs(x)-kappa, 0 ).*sign(x);

% perform outer_iters outer iterations, plot the
% solver progress as the iterations proceed
for k=1:outer_iters,
    fprintf( 1, 'iteration %d of %d\n', k, outer_iters );
    plot( x, 'b-+' );
    drawnow;
 
    % define an anonymous function returning the gradient
    % of the first subproblem w.r.t. x, holding z and u
    % fixed, then perform gradient-descent (Landweber 
    % iterations).
    gradF = @(x) M'*(M*x) - M'*blur + rho*A'*( A*x - z + u );
    x = gradient_descent( gradF, x, step, inner_iters );
    
    % solve the second sub-problem using the anonymous
    % shrinkage operator
    z = shrink( lambda/rho, A*x + u );
    
    % update the scaled Lagrange multipliers
    u = u + A*x - z;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup ADMM and perform iterations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hold on;
plot( blur, 'r-o' );
plot( x, 'b-o' );

Light-in-Flight Paper accepted to Siggraph 2013!

A recent project focused on reducing the cost of Transient Imaging.  Transient images are effectively videos that depict light-transport at the timescales of light-transport; they allow you to see light 'in-flight' as it bounces around the scene. 

This has been done in the past with extremely expensive optical setups costing hundreds of thousands of dollars that involve very short laser pulses and so-called streak-cameras.  However, it is possible to accomplish much the same effect with modifications to off-the-shelf hardware for on the order of only a $1000.  The paper describing this has been accepted to SIGGRAPH 2013 and there will be an accompanying technical demo in the Emerging Technologies area.

You can get an idea what the project does from the following video. If you happen to be attending SIGGRAPH this year, please stop by the booth to find out more!

Stochastic Deconvolution Paper Accepted to CVPR2013

My most recent submission on image deblurring has been accepted as a poster at CVPR2013, to be held in Portland Oregon.  The approach is an adaptation of the stochastic random walk used in my earlier Stochastic Tomography paper to image deblurring and handles saturated pixels and boundary conditions naturally.  The submission video gives an overview of the method:

More info is available at the Stochastic Deconvolution website.  The site also includes sample C++ code illustrating the basic method.

Symmetric Discretization of Linear Boundary Conditions

When solving large-scale PDEs, it's generally preferable (where possible) to have a symmetric discretization so that the conjugate gradient method can be used.  Conjugate gradient uses relatively little memory compared to more general purpose methods (like GMRES) and can also take advantage of symmetric preconditioners (like incomplete LDL^T) that are faster and use lower memory than preconditioners such as incomplete LU.

However applying boundary conditions to the problem can be an issue, since it is easy to end up with a non-symmetric matrix unless the constrained variables are eliminated from the system matrix.  This means that solvers like GMRES have to be used, which is costly in terms of computation and memory.

The following Matlab script gives an example for how to apply linear boundary conditions via an auxiliary boundary condition matrix suggested by a labmate of mine.  It shows both pin (Dirichlet) constraints and gradient (Neumann) boundary conditions at the domain boundary and in the interior.  The only downside is that Neumann constraints are a bit peculiar to specify since one of the variables must be eliminated from the boundary condition matrix.  However this can usually be specified fairly easily, and certainly more easily than eliminating the variable from the system matrix.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Example for showing symmetric discretization  %%
%% of general linear boundary conditions.        %%
%% (c) James Gregson 2013                        %%
%% Please send questions or comments to:         %%
%% james.gregson@gmai.com                        %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% solve a 1D Laplace equaiton with 8 mesh points
% subject to the constraint that the leftmost 
% grid point has value 1.0 and the rightmost 
% grid point has value 2.0 and the fifth is 
% constrained to the value of the fourth minus 0.2
%
% then solves the same problem, but with a biharmonic
% operator, to show that the result is different 
% from just pinning the values; i.e. the biharmonic
% energy is applied to the constrained grid points. 
% If this weren't the case, both problems would have 
% the same solution.

% the approach is to use a standard symmetric 
% discretization of the Laplace operator A for all
% grid points (including those constrained by
% boundary conditions) and a second boundary
% condition matrix that enforces the constraints
% specifically, the boundary conditions are
% represented as:
%
% u = B v + c
%
% where u has the same dimensions as the Laplace
% operator and B has one column for every
% unconstrained grid-point.  c is a vector of
% constants. this allows the problem: 
%
% A u = f
%
% to be written as:
%
% A*B*v = f - A*c
%
% however this is system is rectangular and 
% non-symmetric.  This can be remedied by 
% multiplying through by B', giving:
%
% B'*A*B*v = B'*(f-A*c)
% 
% which is square and symmetric.

% define the objective function, in this case a 
% 6-variable 1D Laplace equation. The Laplacian
% is discretized with a standard 3-point stencil
A = [  1, -1,  0,  0,  0,  0,  0,  0;
      -1,  2, -1,  0,  0,  0,  0,  0;
       0, -1,  2, -1,  0,  0,  0,  0;
       0,  0, -1,  2, -1,  0,  0,  0;
       0,  0,  0, -1,  2, -1,  0,  0;
       0,  0,  0,  0, -1,  2, -1,  0;
       0,  0,  0,  0,  0, -1,  2, -1;
       0,  0,  0,  0,  0,  0, -1,  1 ];

% we are solving the Laplace equation, so the
% right hand side is zero
f = [ 0, 0, 0, 0, 0, 0, 0, 0 ]';

% the boundary condition matrix.  
B = [  0,  0,  0,  0,  0; % first grid point, set to c[1]
       1,  0,  0,  0,  0;
       0,  1,  0,  0,  0;
       0,  0,  1,  0,  0;
       0,  0,  1,  0,  0; % fifth grid point, set to c[4]-0.2
       0,  0,  0,  1,  0;
       0,  0,  0,  0,  1;
       0,  0,  0,  0,  0 ]; % last grid point, set to c[8]

% the boundary condition constant matrix
c = [ 1, 0, 0, 0, -0.2, 0, 0, 2 ]';

% form the system with boundary conditions
M = B'*A*B
rhs = B'*( f - A*c );

% solve for the primary variables
x = M\rhs;

% use the primary variables to solve for the
% grid values
u1 = B * x + c;

% now solve a biharmonic problem with the same constraints
% in the same manner, note the A' factor in the  system
% matrix and right-hand-side.
M = B'*A'*A*B;
rhs = B'*( f - A'*A*c );
x = M\rhs;
u2 = B * x + c;

% plot the results
figure
hold on
plot( u1, '-ro' );
plot( u2, '-b+' );
legend( 'Laplace', 'Biharmonic' );  

The plot produced by this script is shown below:

Both solves reproduce the boundary conditions correctly, but the biharmonic solve finds a solution that is smoother than the Laplace solution around the interior constraints.

Playing with Some Fluid Simulation

A recent project has me playing around with some fluid simulation.  I got some nice looking results and thought I would post a few.

The first is a buoyancy-driven flow simulated using a Marker-And-Cell (MAC) grid with the Boussinesq approximation.  This is simply a source term applied to the vertical velocity based on the fluid density, making it trivial to add to a basic solver.  Although I have no explicitly simulated viscosity, numerical viscosity is enough to cause some interesting unsteady flow effects

The second is a simulation of smoke in box, using a 3D 'Stable Fluids' solver with the first order advection replaced by the second-order 'BFECC' scheme.  This gives -way- better results than standard first order semi-Lagrangian discretizations:

The results compared to the first order scheme (shown below) are pretty impressive:

That's it..

Parametrically Designed Gearboxes

Using my Python CSG and Gear libraries, I've been able to start parametrically designing parts.  This can sort of be done with OpenSCAD, but the lack of proper variables and functions makes it difficult.  Recently I tried designing a full gearbox.  The script is shown below: it's a bit messy, but you get the idea:

import gears
import pyPolyCSG as csg

def make_clamp_hub( B ):
    thickness = 9
    hub = csg.cylinder( B/2.0+6, thickness, True )
    hub = hub + csg.cylinder( 4.0, B+10, True ).rotate( 90.0, 0.0, 0.0 ).translate( B/2+3, 0, 0 )
    hub = hub - csg.cylinder( 2.0, 100, True ).rotate( 90.0, 0.0, 0.0 ).translate( B/2+3, 0, 0 )
    hub = hub - csg.box( B/2+6, thickness, 2, True ).translate( B/2+3.5, 0, 0 )
    return hub.rotate(90,0,0).translate( 0, 0, thickness/2-0.01 )

def make_gear( pressure_angle, pitch, teeth, thickness, bore ):
    px, py = gears.gears_make_gear( pressure_angle, teeth, pitch )
    coords = []
    for i in range( 0, len(px) ):
        coords.append( ( px[i], py[i] ) )
    gear = csg.extrusion( coords, thickness ) + make_clamp_hub( bore ).translate( 0, 0, thickness )
    gear = gear - csg.cylinder( bore/2.0, thickness*100, True ).rotate( 90.0, 0.0, 0.0 )
    return gear

pressure_angle = 20.0
pitch          = 0.8
N1             = 12
B1             = 8.0+0.7
T1             = 10.0

N2             = 36
B2             = 5.0+0.7
T2             = 5.0

backlash       = 1.0

dp1 = gears.gears_pitch_diameter( pressure_angle, N1, pitch )
dp2 = gears.gears_pitch_diameter( pressure_angle, N2, pitch )
do1 = gears.gears_outer_diameter( pressure_angle, N1, pitch )
do2 = gears.gears_outer_diameter( pressure_angle, N2, pitch )
dc  = ( dp1 + dp2 )/2.0 + backlash

gear1 = make_gear( pressure_angle, pitch, N1, T1, B1 )
gear1.save_mesh("gear_%gdeg_P%g_%d_tooth.obj" % ( pressure_angle, pitch, N1 ))

gear2 = make_gear( pressure_angle, pitch, N2, T2, B2 )
gear2.save_mesh("gear_%gdeg_P%g_%d_tooth.obj" % ( pressure_angle, pitch, N2 ))


def hole_xy( x, y, radius, height ):
    return csg.cylinder( radius, height, True ).rotate( 90.0, 0.0, 0.0 ).translate( x, y, height/2.0 )

B1s        = 10.0
B1p        = 22
B2p        = 10

box_pad    = 5.0
screw_diam = 4.0
box_thick  = 4.0
nema_offset  = 31.0/2.0

box_height = max( (do1, do2 ) ) + box_pad*2.0
box_width  = max( ( dc + (max( ( B1p, do1 ) )+do2)/2.0 + box_pad*2.0, nema_offset+dc+dp1/2 + screw_diam*2.0 + box_pad*2.0 ) )
screw_x_off = box_width/2.0-screw_diam
screw_y_off = box_height/2.0-screw_diam

g1_coords = ( box_pad+max((do1,B1p))/2.0, box_height/2.0 )
g2_coords = ( g1_coords[0]+dc, g1_coords[1] )

plate = csg.box( box_width, box_height, box_thick )
plate = plate - hole_xy( box_pad, box_pad, screw_diam/2.0, box_thick*2.0 )
plate = plate - hole_xy( box_width-box_pad, box_pad, screw_diam/2.0, box_thick*2.0 )
plate = plate - hole_xy( box_width-box_pad, box_height-box_pad, screw_diam/2.0, box_thick*2.0 )
plate = plate - hole_xy( box_pad, box_height-box_pad, screw_diam/2.0, box_thick*2.0 )


plate = plate + hole_xy( g1_coords[0], g1_coords[1], B1p/2.0+4.0, box_thick ).translate( 0, 0, box_thick/2+3 )
plate = plate - hole_xy( g1_coords[0], g1_coords[1], B1p/2.0, box_thick ).translate( 0, 0, box_thick/2+3 )
plate = plate - hole_xy( g1_coords[0], g1_coords[1], B1s/2.0, box_thick*2.0 )

plate = plate - hole_xy( g2_coords[0], g2_coords[1], B2p/2.0, box_thick*2.0 )
plate = plate - hole_xy( g2_coords[0]-nema_offset, g2_coords[1]-nema_offset, screw_diam/2.0, box_thick*2.0 )
plate = plate - hole_xy( g2_coords[0]+nema_offset, g2_coords[1]-nema_offset, screw_diam/2.0, box_thick*2.0 )
plate = plate - hole_xy( g2_coords[0]+nema_offset, g2_coords[1]+nema_offset, screw_diam/2.0, box_thick*2.0 )
plate = plate - hole_xy( g2_coords[0]-nema_offset, g2_coords[1]+nema_offset, screw_diam/2.0, box_thick*2.0 )
#plate = plate + gear1.translate( g1_coords[0], g1_coords[1], -10.0 )
#plate = plate + gear2.translate( g2_coords[0], g2_coords[1], -10.0 )

plate.save_mesh( "gearbox_plate.obj" )

This generates two gears and two gearbox plates that are mounted together using M3 screws.  The printed gearbox is shown below, although I'm missing some bearings and hardware.  Each plate has a recess into which a 608ZZ bearing can be press-fit.  Depending on the orientation it can either serve as a fixed or floating end.

I hope to use gearboxes like this to increase the speed that I can drive leadscrews at for my CNC project, which currently has pretty limited feedrates.  This would allow the use of cheap hardware store threaded rods, rather than expensive leadscrews, while still being a self-locking drive.

Python CSG library release

I am releasing the source code to my Python CSG library, named pyPolyCSG.  The code is based on the Carve CSG library and provides Python users with a way to load/save meshes in a variety of formats (.obj, .off, .stl, .vtp) as well as apply transformations (translate/rotate/scale) and perform Boolean operations upon them.  I have now been using the library as a replacement for OpenSCAD for quite some time and finding it useful, as you can see below:

However it is very much a work-in-progress, particularly the installation.  I hope to be cleaning this up and adding features as time goes on.  If anyone has suggestions on how to package the code better, I'd love to hear them.

In the meantime, the code is up on GitHub at: https://github.com/jamesgregson/pyPolyCSG, licensed under the MIT license.  Please send your comments for new features, bug fixes or pretty much anything.

Debugging CImg Code Under XCode

I frequently use CImg to perform basic image operations like loading, saving and resizing. However trying to debug code using the CImg library under XCode has been problematic, since i) CImg uses ImageMagick executables to perform format conversions and ii) XCode doesn't check the system $PATH variable when debugging.

To get around this, you can explicitly specify the path to the ImageMagick executables in your code.  This avoids the command-not-found error that otherwise occurs when using CImg through the debugger.  You can dink around with trying to bend XCode to your will, but I've found this is just way easier.

Here's the related code, note you will need to change the path to suit your system:

cimg::imagemagick_path( "/usr/local/bin/convert" );

I put this at the beginning of the code's main(.) function allowing proper debugging.

2x2 and 3x3 Matrix Inverses and Linear Solvers

Solving 2x2 and 3x3 systems comes up a lot and I have probably rewritten this code dozens of times over the last ten years.  So I'm finally posting it. The header file is below:

/**
 @file solve_2x2_3x3.c
 @author James Gregson (james.gregson@gmail.com)
 @brief 2x2 and 3x3 matrix inverses and linear system solvers. Free for all use. Best efforts made at testing although I give no guarantee of correctness. Although not required, please leave this notice intact and pass along any bug-fixes/reports to the email address above.  See the functions test_2x2() and test_3x3() for usage details.
*/

#ifndef SOLVE_2X2_3X3
#define SOLVE_2x2_3X3

#define SOLVE_SUCCESS 0
#define SOLVE_FAILURE 1
#define SOLVE_EPSILON 1e-8

/**
 @brief Compute the inverse of a 2x2 matrix A and store it in iA, returns SOLVE_FAILURE on failure.
 @param[in] A Input matrix, indexed by row then column, i.e. A[row][column]
 @param[out] iA Output matrix that is the inverse of A, provided a is definite
 @return SOLVE_SUCCESS on success, SOLVE_FAILURE on failure
 */
int invert_2x2( const double A[2][2], double iA[2][2] );

/**
 @brief Solves a 2x2 system Ax=b
 @param[in] A input 2x2 matrix
 @param[out] x output solution vector
 @param[in] b input right-hand-side
 @return SOLVE_SUCCESS on success, SOLVE_FAILURE on failure
 */
int solve_2x2( const double A[2][2], double x[2], const double b[2] );

/**
 @brief Computes the squared residual of a solution x for a linear system Ax=b, for testing purposes.
 @param[in] A input matrix
 @param[in] x input solution vector
 @param[in] b input right-hand-side
 @return squared 2-norm of residual vector
 */
double residual_2x2( const double A[2][2], const double x[2], const double b[2] );

/**
 @brief Generates random systems and solves them using the solve_2x2() function, printing the sum of residual 2-norms.
 */
void test_2x2();

/**
 @brief Compute the inverse of a 3x3 matrix A and store it in iA, returns SOLVE_FAILURE on failure.
 @param[in] A Input matrix, indexed by row then column, i.e. A[row][column]
 @param[out] iA Output matrix that is the inverse of A, provided a is definite
 @return SOLVE_SUCCESS on success, SOLVE_FAILURE on failure
 */
int invert_3x3( const double A[3][3], double iA[3][3] );

/**
 @brief Solves a 3x3 system Ax=b
 @param[in] A input 3x3 matrix
 @param[out] x output solution vector
 @param[in] b input right-hand-side
 @return SOLVE_SUCCESS on success, SOLVE_FAILURE on failure
 */
int solve_3x3( const double A[3][3], double x[3], const double b[3] );

/**
 @brief Computes the squared residual of a solution x for a linear system Ax=b, for testing purposes.
 @param[in] A input matrix
 @param[in] x input solution vector
 @param[in] b input right-hand-side
 @return squared 2-norm of residual vector
 */
double residual_3x3( const double A[3][3], const double x[3], const double b[3] );

/**
 @brief Generates random systems and solves them using the solve_3x3() function, printing the sum of residual 2-norms.
 */
void test_3x3();

#endif

Here is the corresponding source file:

/**
 @file solve_2x2_3x3.c
 @author James Gregson (james.gregson@gmail.com)
 @brief Implementation for 2x2 and 3x3 matrix inverses and linear system solution. See solve_2x2_3x3.h for details.
*/

#include<math.h>
#include<stdlib.h>
#include<stdio.h>
#include"solve_2x2_3x3.h"

int invert_2x2( const double A[2][2], double iA[2][2] ){
 double det;
 det = A[0][0]*A[1][1] - A[0][1]*A[1][0];
 if( fabs(det) < SOLVE_EPSILON )
  return SOLVE_FAILURE;
 
 iA[0][0] =  A[1][1]/det; iA[0][1] = -A[0][1]/det;
 iA[1][0] = -A[1][0]/det; iA[1][1] =  A[0][0]/det;
 
 return SOLVE_SUCCESS;
}

int solve_2x2( const double A[2][2], double x[2], const double b[2] ){
 double iA[2][2];
 
 if( invert_2x2( A, iA ) )
  return SOLVE_FAILURE;
 
 x[0] = iA[0][0]*b[0] + iA[0][1]*b[1];
 x[1] = iA[1][0]*b[0] + iA[1][1]*b[1];
 
 return SOLVE_SUCCESS;
}

double residual_2x2( const double A[2][2], const double x[2], const double b[2] ){
 double r[2];
 r[0] = A[0][0]*x[0] + A[0][1]*x[1] - b[0];
 r[1] = A[1][0]*x[0] + A[1][1]*x[1] - b[1];
 return r[0]*r[0] + r[1]*r[1];
}

void test_2x2(){
 int i, j, k;
 double A[2][2], b[2], x[2], r2=0.0;
 for( k=0; k<10000; k++ ){
  for( i=0; i<2; i++ ){
   b[i] = drand48();
   for( j=0; j<2; j++ ){
    A[i][j] = drand48();
   }
  }
  if( solve_2x2( A, x, b ) == SOLVE_SUCCESS )
   r2 += residual_2x2( A, x, b );
 }
 
 printf( "2x2 squared residual: %0.16E\n", r2 );
 printf( "solution: [%f, %f]^T\n", x[0], x[1] ); 
}

int invert_3x3( const double A[3][3], double iA[3][3] ){
 double det;
 
 det = A[0][0]*(A[2][2]*A[1][1]-A[2][1]*A[1][2])
 - A[1][0]*(A[2][2]*A[0][1]-A[2][1]*A[0][2])
 + A[2][0]*(A[1][2]*A[0][1]-A[1][1]*A[0][2]);
 if( fabs(det) < SOLVE_EPSILON )
  return SOLVE_FAILURE;
 
 iA[0][0] =  (A[2][2]*A[1][1]-A[2][1]*A[1][2])/det;
 iA[0][1] = -(A[2][2]*A[0][1]-A[2][1]*A[0][2])/det;
 iA[0][2] =  (A[1][2]*A[0][1]-A[1][1]*A[0][2])/det;
 
 iA[1][0] = -(A[2][2]*A[1][0]-A[2][0]*A[1][2])/det;
 iA[1][1] =  (A[2][2]*A[0][0]-A[2][0]*A[0][2])/det;
 iA[1][2] = -(A[1][2]*A[0][0]-A[1][0]*A[0][2])/det;
 
 iA[2][0] =  (A[2][1]*A[1][0]-A[2][0]*A[1][1])/det;
 iA[2][1] = -(A[2][1]*A[0][0]-A[2][0]*A[0][1])/det;
 iA[2][2] =  (A[1][1]*A[0][0]-A[1][0]*A[0][1])/det;
 
 return SOLVE_SUCCESS;
}

int solve_3x3( const double A[3][3], double x[3], const double b[3] ){
 double iA[3][3];
 
 if( invert_3x3( A, iA ) )
  return SOLVE_FAILURE;
 
 x[0] = iA[0][0]*b[0] + iA[0][1]*b[1] + iA[0][2]*b[2];
 x[1] = iA[1][0]*b[0] + iA[1][1]*b[1] + iA[1][2]*b[2];
 x[2] = iA[2][0]*b[0] + iA[2][1]*b[1] + iA[2][2]*b[2];
 
 return SOLVE_SUCCESS;
}

double residual_3x3( const double A[3][3], const double x[3], const double b[3] ){ 
 double r[3];
 r[0] = A[0][0]*x[0] + A[0][1]*x[1] + A[0][2]*x[2] - b[0];
 r[1] = A[1][0]*x[0] + A[1][1]*x[1] + A[1][2]*x[2] - b[1];
 r[2] = A[2][0]*x[0] + A[2][1]*x[1] + A[2][2]*x[2] - b[2];
 return r[0]*r[0] + r[1]*r[1] + r[2]*r[2];
} 

void test_3x3(){
 int i, j, k;
 double A[3][3], b[3], x[3], r2=0.0;
 for( k=0; k<10000; k++ ){
  for( i=0; i<3; i++ ){
   b[i] = drand48();
   for( j=0; j<3; j++ ){
    A[i][j] = drand48();
   }
  }
  if( solve_3x3( A, x, b ) == SOLVE_SUCCESS )
   r2 += residual_3x3( A, x, b );
 }
 
 printf( "3x3 squared residual: %0.16E\n", r2 );
 printf( "solution: [%f, %f, %f]^T\n", x[0], x[1], x[2] );
}

The code is based on these formulas for 2x2 and 3x3 matrix inverses using Cramer's rule.  The formulas are just translated to C, which I have found in the past to be an error-prone process. They are not optimized in any way, but are tested using the test_2x2() and test_3x3() functions above. Note that Visual C++ users will need a functioning implementation of drand48(), and will need to add the header defining it to the list of included headers in the source file. OS-X and Linux users will already have this function defined in stdlib.h.

Python Constructive Solid Geometry Update

In a earlier posts I've alluded to a Python Constructive Solid Geometry (CSG) library that I was working on to allow parametric design.  You can do this with OpenSCAD, which is great software, but in my opinion the language leaves a bit to be desired. I wanted a solution that worked with existing languages, specifically C, C++ and Python, so that the results could be integrated easily with other software such as remeshers or FEA packages.

Of course writing a robust CSG library is a daunting undertaking.  Fortunately there are existing libraries such as CGAL and Carve that handle this.  In my opinion CGAL is the more robust of the two, however it currently has compilation issues under OS-X and is substantially slower than Carve.

Regardless, neither have the interface that I'm looking for, like the ability to directly load meshes, affine transformations and minimal-code ways to perform boolean operations on meshes.  So I started work on a C++ wrapper for Carve that would give me the interface I wanted, with a wrapper for Python.

I'm pleased to say that it's coming along quite well and able to produce parts that are non-trivial.  The interface is considerably cleaned up from before and I'm now starting to use it for projects.  Here's two examples from (another) CNC project:

 

The code that generated these models is here:

from pyCSG import *

def inch_to_mm( inches ):
    return inches*25.4

def mm_to_inch( mm ):
    return mm/25.4

def hole_compensation( diameter ):
    return diameter+1.0

mounting_hole_radius = 0.5*hole_compensation( inch_to_mm( 5.0/16.0 ) )


def axis_end():
    obj = box( inch_to_mm( 4.5 ), inch_to_mm( 1.75 ), inch_to_mm( 0.75 ), True )
    
    screw_hole = cylinder( mounting_hole_radius, inch_to_mm( 3.0 ), True, 20 )
    
    shaft_hole = cylinder( 0.5*hole_compensation( inch_to_mm( 0.5 ) ), inch_to_mm(1.0), True, 20 ).rotate( 90.0, 0.0, 0.0 )

    center_hole = cylinder( 0.5*hole_compensation( inch_to_mm( 1.0 ) ), inch_to_mm(1.0), True, 20 ).rotate( 90.0, 0.0, 0.0 )
    mount_hole = cylinder( 0.5*hole_compensation( 4.0), inch_to_mm(1.0), True, 10 ).rotate( 90.0, 0.0, 0.0 )
    
    notch = box( inch_to_mm( 1.5 ), 2.0, inch_to_mm( 1.0 ), True )

    
    obj = obj - ( shaft_hole.translate( inch_to_mm( 1.5 ), 0.0, 0.0 ) + shaft_hole.translate( inch_to_mm( -1.5 ), 0.0, 0.0 ) )
    obj = obj - ( notch.translate( inch_to_mm( 2.25 ), 0.0, 0.0 ) + notch.translate( inch_to_mm( -2.25 ), 0.0, 0.0 ) )
    obj = obj - ( center_hole + mount_hole.translate( -15.5, -15.5, 0.0 ) + mount_hole.translate(  15.5, -15.5, 0.0 ) + mount_hole.translate( 15.5, 15.5, 0.0 ) + mount_hole.translate( -15.5, 15.5, 0.0 ) ) 
    
    obj = obj - ( screw_hole.translate( inch_to_mm(1.0), 0.0, 0.0 ) + screw_hole.translate( inch_to_mm(-1.0), 0.0, 0.0 ) )
    obj = obj - ( screw_hole.translate( inch_to_mm(2.0), 0.0, 0.0 ) + screw_hole.translate( inch_to_mm(-2.0), 0.0, 0.0 ) )
    
    return obj

def carriage():
    obj = box( inch_to_mm( 5 ), inch_to_mm( 5 ), inch_to_mm( 1.0 ), True )
    shaft_hole = cylinder( inch_to_mm( 0.75 )/2.0, inch_to_mm( 5.5 ), True )
    screw_hole = cylinder( inch_to_mm( 0.5 )/2.0, inch_to_mm( 5.5 ), True )
    
    leadnut_hole = cylinder( inch_to_mm(0.25)*0.5, inch_to_mm( 1.0 ), True );
    leadnut_access = box( inch_to_mm( 1.5 ), inch_to_mm( 3.0/8.0 ), inch_to_mm( 1.0 ), True )
    
    
    mhole = cylinder( mounting_hole_radius, inch_to_mm( 2.0 ), True ).rotate( 90.0, 0.0, 0.0 )
    
    obj = obj - ( shaft_hole.translate( inch_to_mm( 1.5 ), 0.0, 0.0 ) + shaft_hole.translate( inch_to_mm( -1.5 ), 0.0, 0.0 ) + screw_hole )
    obj = obj - ( leadnut_hole.translate( inch_to_mm( 0.5 ), inch_to_mm( -2.5 ), 0.0 ) + leadnut_hole.translate( inch_to_mm( -0.5 ), inch_to_mm( -2.5 ), 0.0 ) + leadnut_access.translate( 0.0, inch_to_mm( -2.0 ), inch_to_mm( 0.2 ) ) )
    
    for i in range( -2, 3 ):
        for j in range( -2, 3 ):
            if i != 0 and j != 0:
                obj = obj - ( mhole.translate( inch_to_mm( 1.0*i ), inch_to_mm( 1.0*j ), 0.0 ) )
    return obj
              

axis_end().save_mesh("axis_end.obj" )
carriage().save_mesh("carriage.obj" )

As you can see, this approach gives lots of flexibility in terms of manipulating and patterning objects using custom code.  The examples above are not great examples of parametric design, but I'm sure you can imagine the sort of stuff that can be done.

I still have to perform a bit of cleanup outside the library to get printable models.  I just run each model through MeshLab's planar edge-flipping optimizer. This is a pretty simple step and I plan to integrate it into the library shortly, along with the ability to extrude custom profiles and build surfaces of revolution.  When these features are finished I plan to release the code for the library and Python wrapper.