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 };

Updated C Mathematical Expression Parser

I've updated and posted my recursive descent mathematical expression parsing code.  The new code is available from my Github repository:

https://github.com/jamesgregson/expression_parser

The original post describing the library is at: http://jamesgregson.blogspot.ca/2012/06/mathematical-expression-parser-in-c.html

New features include Boolean operations and a callback interface for custom named variables and functions.  You can optionally disable the Boolean operations (and the 10 levels of parsing they trigger) if not needed.

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!

pyPolyCSG library updated

I have updated the pyPolyCSG python constructive solid geometry library to use the most recent version of the Carve CSG library. The updates to Carve appear to make it significantly more robust as well as much faster.

The update to pyPolyCSG:

• Improves the robustness, at least on my existing scripts
• Improves the speed of performing Boolean operations
• Fixes a bug where vertex normals and texture coordinates were treated as vertices (thanks to Ryan Rix for finding this)
• Adds the ability to transform polyhedra by a 3x3 or 4x4 matrix (again thanks to Ryan Rix!)

You can get the updated version from github: https://github.com/jamesgregson/pyPolyCSG

Smooth Feedrate Envelopes for Motion Control, Part II

In the previous post, I derived equations for smooth feedrate control of stepper motors, claiming that by using a smoother feedrate envelope that the motors could be driven faster with less chance of skipping steps and losing position.

In this post, I demonstrate this with a real stepper motor and show that it actually does work: using the envelopes does actually prevent the motors from losing steps.  My test setup is a single NEMA 17 stepper, driven by one of my A4988 driver breakouts, which is controlled by an Arduino sketch running on the Arduino Due.  I'm using half-stepping on the motors, driven by a 200KHz timer interrupt step callback which decides whether or not to step based on the interpolated supplied delays for the start and end of each move.  The move itself approximates a square wave, first accelerating from a slow feedrate, then performing a constant speed portion, then decelerating back to the initial feedrate.

The video below shows the stepper being driven using a constant acceleration profile, which causes the kinks in the feedrate graph in my previous post. You can clearly see it moving around on the table and stalling frequently before it reaches the top speed.


In contrast, here is the result using the third-order cubic feedrate envelope for the same set of moves. The stepper is easily able to handle the top speed and jerks around considerably less on the table.  Of course this comes at a price, a higher pulse-frequency must be used to resolve the acceleration profile.


You can get the code I used for this from following link: https://sites.google.com/site/jamesgregson/tmp/linear_move.zip, it includes a multi-axis DDA implementation suitable for use with timer-interrupts as well as the code for evaluating the feedrate envelopes.

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' );

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.

Setting up the WiFly RN-XV with a Teensy 3.0

A recent order from Sparkfun arrived, including a 3.3V Serial LCD and a Roving Networks RN-XV WiFly module.  The RN-XV module is intended to be a drop-in replacement for an XBee, except that it operates over WiFi.  At about $35, it is just about the cheapest way to make your project wireless enabled.

The module is 3.3V, meaning some form of level shifting is needed with a 5V system like an Arduino.  You can use this module with an Arduino via an XBee shield pretty easily.  However it is even easier to use with the 3.3V Teensy 3.0 ARM board, provided you have a breakout for small-pitch XBee module footprint.  The Teensy is also nice for this application because it has multiple serial ports, so you don't need to use the SoftwareSerial library, or program the board, then disconnect to use the wireless.

Setting everything up was pretty easy once I knew what to do, but this post summarizes the process should I ever need to do it again.

The setup that I am using is shown below:

Only four connections are needed once you're set up, 3.3V, GND and two data connections.  DOUT from the Teensy3.0 Serial3 connects to DIN of the RN-XV and DIN from the Teensy to DOUT to the RN-XV.  This makes the module operate as just a serial port, making it pretty easy to interface with. The remaining orange wire connects the Serial LCD display, more on this later.

To get started, I found it was easiest to set the RN-XV in ad-hoc mode. This can be done by connecting pin 8 to 3.3V and will cause the module to create its own wireless network.  When this happens you will see the status LEDs blinking green, orange and red; they're doing it, but you can't really see in this picture.  Note the additional green wire to 3.3V connected to the 8th pin.

You can then look for the network. On a Mac it's pretty easy, it just shows up in the list of networks in the status bar:

The WiFly shows up towards the bottom as WiFly-GSX-a8 or something similar.  If you connect to this network, you can then telnet to the module using the IP address: 169.254.1.1, port 2000.  The module should then respond with a *HELLO* string, at which point you type $$$ to enter command mode.  Command mode allows you to set up the module for your network.

When the module is ready, it will respond with the CMD message to indicate that you're in command mode.  To set up your network you can issue the commands:

set wlan phrase (password);
set lan ssid (your network name);
save
reboot

You can also issue commands to assign a static IP address to the module, but I didn't do this.  For more information, see this excellent introduction http://www.tinkerfailure.com/2012/02/setting-up-the-wifly-rn-xv/

I found that sometimes the module would respond with a confirmation and sometimes would not. I repeated the process a few times in the hopes that some combination would stick.  After this process, remove the power and and connection from pin 8 to 3.3V.  This will cause the device to try to connect to your wireless network.

You should now be able to telnet to the device, but this time with your computer and it connected to your normal WiFi network rather than the ad-hoc network that the device creates.  However first you need to find the IP address of the module.  To do this, I went into my router configuration page:

Conveniently the WiFly module had an entry: 192.168.1.106. Depending on your router, you should be able to set up a specific IP address for the router to assign to the module based on the MAC address.  However my POS router does not allow this.

I could then telnet to the module's IP address, again using port 2000.  This module responds with the same *HELLO* prompt, indicating that everything was successful and the module is on the network and communicating.

With the connections above the Teensy should now see the module as just another serial port.  To test this, I attached the Serial LCD and uploaded the following code to the Teensy:

#include<stdio.h>

void setup(){
  Serial.begin(9600);
  Serial2.begin(9600);
  Serial3.begin(9600);
}

void write_lines( const char *L0, const char *L1 ){
  
  Serial2.write( 0xFE );
  Serial2.write( 0x01 );
  delay(10);
  Serial2.write( 0xFE );
  Serial2.write( 128 );
  delay(10);
  Serial2.print( L0 );
  Serial2.write( 0xFE );
  Serial2.write( 192 );
  delay(10);
  Serial2.print( L1 );
}


void loop(){
  if( Serial3.available() ){
    char L0[17];
    char L1[17];
    int pos = 0;

    L0[0] = '\0';
    L1[0] = '\0';

    while( Serial3.available() ){
      char c = Serial3.read();
      if( c == '\n' ){
        pos = 0;
        Serial.print('\n');
      } else if( c == '\r' ){
        
      } else {
        if( pos < 16 ){
          L0[pos] = c;
          pos++;
          L0[pos] = '\0';
        } else if( pos < 32 ){
          L1[pos-16] = c; 
          pos++;  
          L1[pos-16] = '\0'; 
        }
        Serial.print( (char)c );
      }
    }
    write_lines( L0, L1 );
  }
  delay(100);
}

My LCD is a 2x16 character display.  The code above just polls for available data on the third serial port and, when a newline is encountered, prints it out onto the display.  Lo and behold, after the following session:

Jamess-MacBook-Pro:~ jgregson$ telnet 192.168.1.106 2000
Trying 192.168.1.106...
Connected to 192.168.1.106.
Escape character is '^]'.
*HELLO*
This is James

The result on the display is below:

Hooray! An utterly useless internet thingy!

Putting CMake binaries in a specific directory

I've started using CMake for pretty much all my development work and have been quite happy except for a few minor issues.  One of those is that the build-products are put in Debug and Release directories by default, instead of single output directory. You can add an install target to do this, but XCode won't apply it through the UI. Ditto goes for a post-build copy step.

It took a while to find, but it's actually simple to get binaries placed into a directory of your choosing. In the example project, code and the CMakeLists.txt file is in the code/ subdirectory, I build using CMake/XCode in the build/ subdirectory and I want the build products placed in the bin/ subdirectory, regardless of whether they are built as debug or release.

Here's the CMakeLists.txt file:

cmake_minimum_required( VERSION 2.6 )

project( entropy )

add_executable( entropy main.cpp )
target_link_libraries( entropy ${LIBS} )
set_target_properties( entropy PROPERTIES
  RUNTIME_OUTPUT_DIRECTORY_DEBUG   ${CMAKE_SOURCE_DIR}/../bin 
  RUNTIME_OUTPUT_DIRECTORY_RELEASE ${CMAKE_SOURCE_DIR}/../bin
)

The set_target_properties() command allows you to set the runtime output directory. You'll see this a fix similar to this online in many places, but it will put the built executables in bin/Debug or bin/Release, instead of just bin/. What was not immediately obvious is that it can be set per-build-configuration, causing all build products to be put in bin/.

Start of CNC Firmware Running on the Teensy 3.0

I've been intermittently working on some CNC firmware with the goal of getting it to run on everything from an Arduino to a standard PC.  It's a fun project because it's very resource constrained but involves a number of different parts, like GCode parsing, asynchcronous programming, portability and motion control for non-Cartesian machines.

In the interests of portability I've been writing the various components in plain ANSI C (except the comments, I like my double slashes!).  By doing this I expect the code to be portable across a wide variety of platforms, from AVRs to the Propoeller child to ARM to PCs, with minimal effort.

My hopes are that the firmware will conform to the NIST GCode specification as closely as is manageable.  That said, some sacrifices will have to made as the spec calls for around 20k of addressable space for expressions, which exceeds the total RAM of Arduinos and even the Teensy 3.0.  However I hope to get mostly spec-compliant in terms of order of operations, expressions (if not the addressable space) and features.  Where possible I am also allowing the specific capabilities to be determined by preprocessor macros.

I've done a bit of work on the GCode interpreter, getting some test code together that parses GCode expressions (this led to my Mathematical Expression Parser in C) as well as an arbitrary dimension DDA implementation.  To date these have been just simple tests, nothing actually running on real hardware.  I've also done some early tests on doing feedrate optimization and

However today I made some good progress on actually getting some real firmware started.  I got my early GCode interpreter stripped down and abstracted out all the hardware-specific stuff into a hardware layer.  I then wrote the hardware abstraction layer for the Teensy 3.0 version of the Arduino environment and fired up the code, using the following string to test with:

N1229.0(This is a comment)G01x110.0   y330.0 M7 F20.0(MSG: another comment)

And it worked!

The screenshot above shows the GCode interpreter code loaded up in the Arduino environment, the Teensy loader application that programs the Teensy and the serial output of the parser. Note that this is not just reading the commands linearly, it is fully parsing the input, splitting it into distinct commands, sorting the commands by operation precedence (as described in the NIST spec) and finally calling back to the Arduino sketch with each command.

Parsing the string and calling back takes about 250 us (the remaining time taken is up by serial communications) so more than fast enough to fill up a lookahead buffer.  Total RAM used is about 5k, so this won't run on an Uno but does run handily on the Teensy 3.0. Ditto goes for the flash at 33Kb currently. This is pretty hefty compared to other firmware, but I think portability and cleanness makes up for that.

The next steps will be to get the kinematic/motion control stuff fleshed out.  I think I should be able to use by DDA implementation running with the Teensy timers pretty easily.  Then it's just tweaking and adding features...

Periodic Interrupt Timers on the Teensy 3.0 (Freescale MK20DX128)

I recently ordered a Teensy 3.0 and today it finally arrived!  It definitely is teensy.  Anyway, after soldering on headers and popping it in a breadboard, I got it running with the blink example.  It took some time, but my issue was downloading directly from the Teensy loader page (which does not support Teensy3.0) rather than from the PJRC Teensy forums. To be fair, the Teensy loader page does have a notice at the top about this, but I missed it.

After getting the software, I was able to write a standard Arduino sketch to blink the LED connected to pin 13.  The code is below:

void setup(){
   pinMode(13, OUTPUT);
}

void loop(){
   digitalWrite( 13, HIGH );
   delay( 100 );
   digitalRead( 13, LOW );
   delay( 100 )
}

It's pretty nice that the Teensy3 works like a regular Arduino, but with more pins, at a higher clock rate, in 32 bits and with heaps of extra peripherals.  However my end goal is to use the board as a CNC controller for my ongoing firmware project.  This will make extensive use of timed interrupts to control steppers, so I thought I'd try to get timers working.  Of course, the Teensy3 is not a Atmel uC, so everything changes at this point and, with the help of this forum post to get started, I had to dive into the manual (available from here) for the Freescale MK20DX128 that the Teensy3 is based upon.

According to the forum post, the timer to use is one of the (4?) Periodic Interrupt Timers (PITs).  As expected, these have a number of control registers. Registers listed with an [N], e.g. PIT_LDVAL[N] should have an appropriate timer index substituted, like PIT_LDVAL2.  Here are the registers:

• SIM_SCGC6 - Enables/disables clock used by PIT timers, not exactly clear on the details, set to SIM_SCGC6_PIT in the forum post example.
• PIT_MCR - Enables and disables the PIT timers. Writing zero enables the timers and writing 1 disables them.
• PIT_LDVAL[N] - Sets the timer count value.  Apparently the timer runs at 50MHz, so toggling timer 2 every second should set PIT_LDVAL2 to 0x2fa080 (hex for 50,000,000).  Visually, this appears to be around a second.
• PIT_TCTRL[N] - Bit zero (TEN in the manual) enables (set to 1) or disables (set to zero) the timer. Bit one (TIE in the manual) enables (set to 1) or disables (set to zero) interrupts that can be generated by the timer.
• PIT_TFLG[N] - Flag to indicate timer waiting.  Set to one to start timer and at the end of every called interrupt routine, otherwise interrupts will stop. 

Finally, interupts must be enabled. Again I'm not clear on the details, but calling NVIC_ENABLE_IRQ( IRC_PIT_CH[N] ) results in the interupt "void pit[N]_isr(void){}" being called.  Although it seems like the chip should have four timers, I only succeeded in getting timers 0, 1, and 2 working properly with interrupts, testing with index 3 gave a linker error in the Arduino software.

Anyway, here's the code for my tests:

#define TIE 0x2
#define TEN 0x1

void pit0_isr(void){
  digitalWrite( 13, !digitalRead(13) );
  PIT_TFLG0 = 1; 
}

void setup(){
  pinMode(13,OUTPUT);
  SIM_SCGC6 |= SIM_SCGC6_PIT;
  PIT_MCR = 0x00;
  NVIC_ENABLE_IRQ(IRQ_PIT_CH0);
  PIT_LDVAL0 = 0x2faf080;
  PIT_TCTRL0 = TIE;
  PIT_TCTRL0 |= TEN;
  PIT_TFLG0 |= 1;
}

void loop(){
  delay(2000);
}

Hope this helps someone get up to speed, and perhaps serves as a reference for me later on.

CMake add_subdirectory link errors

This is probably something that everyone already knows, but it took me a few hours to figure out...

I'm currently working on a project that uses CMake to build several libraries and executables all in the same project.  Everything worked fine on OS-X, but when I tried to build on Linux I kept getting linker errors when compiling executable using functions and objects from interdependent libraries that I thought were being linked correctly via the target_link_libraries() command.  Each library was compiled by CMake using the add_subdirectory() command.

I found out that the issue is that you need to specify dependencies between the libraries themselves, not just between executables and libraries.  This was a bit of a surprise, but just adding the target_link_libraries() command to the CMakeLists.txt file for each library and listing the dependencies fixed the problem.

Simplified DDA/Bresenham-like algorithm for line drawing in 2D and 3D

Implementing Bresenham's line drawing algorithm is a pain and has some drawbacks for motion control applications. In particular, it relies on swapping endpoints of the line-segments to achieve specific preconditions and has eight configurations (in 2D alone!) that must be implemented to draw arbitrarily oriented lines.  This has big disadvantages for implementation in high dimension (6D in my case) and even bigger disadvantages for motion control, where you want to move from point A to point B, not arbitrarily (and incorrectly) assume you're already at B and then move backward to A because, hey, order doesn't matter!

I recently ran across Cris Luengo's blog post on a Bresenham-equivalent algorithm in arbitrary dimension based on floating point rounding.  This is likely to be faster on modern machines than the integer only version, simply because i) floating point ops and casting are cheap and ii) branches are expensive, which is effectively the opposite of the case when the original Bresenham algorithm was developed.  He claims the method reproduces the original Bresenham algorithm's plotted pixels, although in my tests I seem to get occasional errant pixels even with double precision.  It may well be an error on my part.

Regardless, this post gave me the idea of just implementing Bresenham's algorithm in a parametric form, accumulating error and advancing pixels in both (all) coordinates in the same manner.  This is considerably more costly in 2D, where you double the operations required, but the relative cost decreases as the dimension increases. Since I'm targeting 6D, the relative cost is only 1/5th, and the simplicity of implementation and lack of endpoint switching more than makes up for it.

I've put up simple C++ code for 2D and 3D integer-only line drawing using a test implementation that is freely available for any use.  It draws 50 random lines with the integer only algorithm described above into the green channel and the same 50 lines with the floating point rounding method (in double precision) into the red channel. The results are below, if the algorithms matched exactly, all the lines would be yellow, with no red or green pixels anywhere.

It looks okay at first, but if you zoom in and look carefully you can see a few red and green pixels around. Regardless, they both do a good job of plotting lines, and if my code is off by a step or two occasionally it's not the end of the world.

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.

An End to Captive Nuts

I've been using 3D printed gears and timing pulleys for a while now, but have been very disappointed with the captive nuts that are used by most scripts available on Thingiverse. Generally I've found that the material being printed isn't stiff enough to allow the set screws to be tightened enough to secure to the shaft without deforming. This causes the pulley or gear to deform out of true.

To get around this I've started using clamping hubs.  They're not too much larger than the captive nut hubs but deform uniformly so the pulley/gear run true.

The picture above shows an example gear.  I've used an M3 screw in the clamping hub.  It ends up fixing very securely to the 8mm shaft and has the additional advantage of not marking the shaft.

I printed the gear above after using my Python Involute Gear Script to generate the involute profile, followed by the Python Constructive Solid Geometry Library to generate the hub and 3D model. The full source of the script that I used is shown below:

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

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

I've found that being able to use Python is much more convenient than OpenSCAD, mostly because its possible to define (real) variables and functions/classes.  As a result I've pretty much switched to using the Python CSG library from OpenSCAD.

Read text file into std::string

Taken (unmodified) from a StackOverflow reply, but I use this frequently and always have to search for it.

std::ifstream in("file.txt");
std::string contents((std::istreambuf_iterator<char>(in)), 
    std::istreambuf_iterator<char>());