syntax-highlighter

Showing posts with label random. Show all posts
Showing posts with label random. Show all posts

Wednesday, June 18, 2014

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

Sunday, December 8, 2013

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

Tuesday, October 22, 2013

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

Thursday, June 13, 2013

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!

Friday, May 31, 2013

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

Sunday, May 26, 2013

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.

Thursday, May 23, 2013

Smooth Feedrate Envelopes for Motion Control

When running motors it is desirable to run them as smoothly as possible to minimize vibrations and possible missed steps. This is why controllers for 3D printers and CNC machines typically incorporate some notion of acceleration rather than instantly switching from one feedrate to another.  Often this is done with a simple ramp of the feedrate, i.e. using a constant acceleration profile.

As an example, consider a machine starting at feedrate $f_0$ and then performing a very long linear move at a constant feedrate $f_1$. For this example, If $f(t)$ is the machine feedrate over time and $a$ is a constant acceleration, this profile would be defined mathematically as follows:

\begin{equation}
f(t) = \begin{cases}
f_0      & \mbox{if } t < t_0 \\
f_0 + a (t - t_0) \hspace{0.5cm} & \mbox{if } 0 \leq t - t_0 \leq \frac{f_1-f_0}{a} \\
f_1 & \mbox{otherwise}
\end{cases}
\end{equation}

Graphically, here's a plot of the feedrate over time for a move starting at rest at $t=0$ and accelerating up to a feedrate of 2 over one unit of time:


The problem with a constant acceleration profile is that there are sharp kinks in the feedrate plotted over time.  These kinks imply instantaneous changes in acceleration, which in turn imply infinite forces for infinitely short periods of time.  Of course, there is no mechanical way to produce these forces, so what actually happens is the machine overshoots very slightly and averages the forces out over a short time. For low feedrates with a light machine this actually works okay, but for a heavy machine at high-feedrates, the overshoot can be more than a motor step which causes the machine to lose position.  In an open-loop design, once the machine loses position, it never recovers and in all likelihood, the part is ruined.

There are a few ways to address this problem:
  • Use lower feedrates
  • Use higher torque motors
  • Use a closed-loop control scheme, e.g. with encoders on the motors
  • Make the acceleration smooth
The first is clearly not an option because it wastes time and feedrates may be chosen specifically for valid reasons such as minimizing local part heating or reducing machining time.  In an ideal world we'd do the remaining three items, but options two and three are expensive, particularly for hobby gear.  However the fourth option can be tackled in firmware with minimal hardware overhead. 

In order to smoothly transition between accelerations we can simply use a different curve to interpolate the feedrates.  The conditions needed are that the feedrates match the desired rates at the beginning and end of the curve and that the slope of the feedrate curves (i.e. the acceleration) is zero at the endpoints.  In between the endpoints we want the curve to be smooth.

The one of the simplest classes of functions that meet these requirements are cubic polynomials.  These are defined by four coefficients $a$, $b$, $c$ and $d$ using the following equation, where $\tau$ is the fraction of the total time spent accelerating:

\begin{equation}
f(\tau) = a \tau^3 + b \tau^2 + c \tau + d
\end{equation}

We now want to solve for the coefficients needed to reproduce the move.  There are four coefficients so we need four equations.  Two come from the requirement that we match the feedrates at the curve endpoints:

\begin{eqnarray}
f(\tau=0) = a 0^3 + b 0^2 + c 0 + d &=& f_0 \\
f(\tau=1) = a 1^3 + b 1^2 + c 1 + d &=& f_1
\end{eqnarray}

From these, we see that $d=f_0$ and $a+b+c=f_1-f_0$. The remaining two equations can be found using the requirements that the slope of the feedrate curve is zero at the endpoints. To enforce these constraints we need the derivative of the cubic function:

\begin{equation}
f'(\tau) = 3 a \tau^2 + 2 b \tau + c
\end{equation}

The constraints can now be enforced by requiring that:

\begin{eqnarray}
f'(\tau=0) = 3 a 0^2 + 2 b 0 + c &=& 0 \\
f'(\tau=1) = 3 a 1^2 + 2 b 1 + c &=& 0 
\end{eqnarray}

These equations make it clear that $c=0$ and $3 a + 2 b = 0$. Combining these with the previous conditions leaves two equations and two unknowns:

\begin{eqnarray}
a + b &=& f_1 - f_0 \\
3 a + 2b &=& 0
\end{eqnarray}

So $a = -\frac{2 b}{3}$ which means that $b = 3 (f_1-f_0)$ and $a = -2 (f_1 - f_0)$. This gives the following equation for the interpolating curve:

\begin{equation}
f(\tau) = -2(f_1-f_0)\tau^3 + 3(f_1-f_0)\tau^2 + f_0
\end{equation}

The only remaining thing is to define $\tau$ in terms of $t$.  This is a simple linear interpolation from the start of the acceleration $t_0$ to the end of the acceleration $t_1=\frac{f_1-f_0}{a}$:

\begin{equation}
\tau = \frac{t-t_0}{t_1-t_0} = \frac{t-t_0}{\frac{f_1-f_0}{a}-t_0}
\end{equation}

Plotting this for the same parameters as before gives a smooth, kink-free curve that considerably reduces the time-rate-of-change of acceleration:


In the post-to-come I will demonstrate applying this to a real stepper motor being driven aggressively.  Although seemingly complicated, for a cost of only a few operations per step, it is possible to switch from the linear acceleration profile to the cubic one derived here and get considerably smoother operation.

Sunday, March 24, 2013

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/.

Friday, March 22, 2013

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.

Monday, September 10, 2012

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.

Tuesday, September 4, 2012

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.

Sunday, June 10, 2012

3-Axis A4988 Stepper Driver Carrier

More on the perpetually in-progress CNC (see the mechanical stuff in parts I, II, III, IV)

I've recently built a prototype board for the Pololu A4988 stepper drivers carriers.  These little drivers are inexpensive (about $12) and fairly gutsy (2A per phase), but can blow fairly easily.  Using male headers they can easily be used as drop in modules for a larger CNC controller board.  My board breaks out the microstepping pins to DIP switches and adds screw terminal connections for the high-power side, with female headers for the TTL control inputs and low-power supply.


The two-pin set of female headers on the left is the low-power supply, the middle six-pin set of headers is the step/direction controls for each of the axes and the bottom six-pin header is for upper and lower limit switches for each axis.  The top set of two-pin screw-terminals is the motor power-supply and the remaining 3x6 pin connectors are the motor winding connections with the bottom two terminals for upper and lower limit inputs.  Unfortunately I ran out of space on the board to provide a 5V supply to the limit switches, so these will have to be wired externally to 5V.  I may also add some pull-down resistors to the backside of the board, since the switches are currently floating, although this does not seem to be a problem in practice for some reason.

There's a surprising number of solder joints needed for this simple board, largely due to using point-to-point wiring, but it seems to be electrically sound:


This board really cleans up my testing rig for the CNC, just a few jumper wires are all that's needed to connect my Arduino to the two mostly-finished translation stages.  Using screw-terminals instead of soldered on connectors is also much more convenient when rewiring the steppers as bipolar parallel/serial during testing:


The final machine will probably not end up using this board, I'm considering investing in a proper 3-axis board, since they can be found quite cheaply on Ebay.  That said, it's a nice tidy little package that allows the stepper drivers to be replaced as needed as well as preventing the inevitable wiring mistakes that happen when developing on breadboards.

Thursday, May 31, 2012

Boost::python + Distutils Python Extension Module Example

In getting my ongoing python CSG library functioning, I needed a way to compile and build python extension modules.  Typically I use qmake from the QT library to build projects, as it's relatively cross-platform (Xcode support is broken as of version 4+) and has relatively little magic sauce making it function (I like transparency in a build system).

However, this isn't so hot for python modules where it's nice to have a simple setup.py script that compiles and installs your module.   Distutils provides just this functionality, while being pretty easy to pick up.  So I decided to do the building with distutils.

Then you have to choose how to actually generate the wrapper code.  Doing this by hand using the python API is manageable but cumbersome and results in a lot more code to maintain.  To improve this, there are a number of other wrapper code generators, including SWIG and boost::python.  I prefer boost::python because it i) is very explicit, you define the exposed interface directly and ii) works with all compilers, without custom build steps or interface file compilation.  It's a bit more typing than SWIG, but I (again) like the transparency.

Of course the downside to boost::python is bjam, the boost build system that appears to summon demons from the ether to build your project through sorcery.  Most tutorials assume you will use bjam but I have yet to get it to work on an even modestly complicated project.  So I set about building python modules wrapped by boost::python using Distutils.  This is actually pretty easy, but I had to dig around through forums and documentation to find out how. Consequently I decided to post a full example from start to finish.

I'm going to assume that we're wrapping two source files functions_to_wrap.cpp and classes_to_wrap.cpp, located in root_dir/src, which have headers functions_to_wrap.h and classes_to_wrap.h located in root_dir/include.  The two headers are:

functions_to_wrap.h
#ifndef FUNCTIONS_TO_WRAP_H
#define FUNCTIONS_TO_WRAP_H

// returns a random string
const char *get_string();

// returns true if values are equal
bool are_values_equal( int a, int b );

// returns the number of supplied arguments to demonstrate 
// boost::python's default argument overloading features
int num_arguments( bool arg0, bool arg1=false, bool arg2=false, bool arg3=false );

#endif

classes_to_wrap.h
#ifndef CLASSES_TO_WRAP_H
#define CLASSES_TO_WRAP_H

class wrapped_class {
private:
 int m_value;
public:
 wrapped_class();
 wrapped_class( const int value );
 
 void set_value( const int value );
 int  get_value( ) const;
};
#endif

There are also the implementation files, but these don't actually matter as far as the build process or wrapper definitions are concerned.  To wrap these, we write an additional c++ file in root_dir/src called boost_python_wrapper.cpp that will define the python interface.  The contents of this file are:


boost_python_wrapper.cpp
#include<boost/python.hpp>

// include the headers for the wrapped functions and classes
#include"functions_to_wrap.h"
#include"classes_to_wrap.h"

using namespace boost::python;

// generate overloads for default arguments, this example is a function
// that takes a minimum of 1 argument and a maximum of four. the 
// overloads are generated in num_arguments_overloads, which is 
// then supplied to boost::python when defining the wrapper
BOOST_PYTHON_FUNCTION_OVERLOADS( num_arguments_overloads, num_arguments, 1, 4 );

// generate wrapper code for both the functions and classes
BOOST_PYTHON_MODULE(ExtensionExample){
 // declare the functions to wrap. note inclusion of the
 // num_arguments_overloads() to include the overloads 
 // that handle default function arguments
 def( "get_string",       get_string );
 def( "are_values_equal", are_values_equal );
 def( "num_arguments",    num_arguments, num_arguments_overloads() );
 
 // declare the classes to wrap. both the default
 // constructor and the constructor taking an 
 // integer argument are exposed
 class_<wrapped_class>("wrapped_class",init<>())
 .def(init<int>())
 .def("get_value", &wrapped_class::get_value )
 .def("set_value", &wrapped_class::set_value )
 ;
}

This file calls a bunch of boost::python macros which instantiate the wrapper code, including handling default arguments for functions and overloaded constructors for classes.  It is also possible to do function overloading, although I've left this out for simplicity.  Note that there's nothing specific about the filenames I chose, as long as the functions and classes to be wrapped are included in the wrapper source code file, you can use any structure/names you want.

When boost_python_wrapper.cpp, functions_to_wrap.cpp and classes_to_wrap.cpp are compiled as a shared library, boost::python will generate wrapper automatically generate wrapper code inside of classes_to_wrap.cpp by macro-expansion, so no additional source files are needed.  Just link the results with the python and boost::python libraries and you get a python module. 

Distutils provides a clean way of doing this and installing the module automatically.  To make a build/install script, we create a script setup.py in the root directory. The contents of this script are little more than compiler and linker flags.

Here's the setup.py script:
from distutils.core import setup, Extension

# define the name of the extension to use
extension_name    = 'ExtensionExample'
extension_version = '1.0'

# define the directories to search for include files
# to get this to work, you may need to include the path
# to your boost installation. Mine was in 
# '/usr/local/include', hence the corresponding entry.
include_dirs = [ '/usr/local/include', 'include' ]

# define the library directories to include any extra
# libraries that may be needed.  The boost::python
# library for me was located in '/usr/local/lib'
library_dirs = [ '/usr/local/lib' ]

# define the libraries to link with the boost python library
libraries = [ 'boost_python-mt' ]

# define the source files for the extension
source_files = [ 'src/boost_python_wrapper.cpp', 'src/functions_to_wrap.cpp', 'src/classes_to_wrap.cpp' ]

# create the extension and add it to the python distribution
setup( name=extension_name, version=extension_version, ext_modules=[Extension( extension_name, source_files, include_dirs=include_dirs, library_dirs=library_dirs, libraries=libraries )] )

Note that the include and library paths are most likely included by default. I have added them simply to demonstrate how you would include (potentially custom) headers and libraries into the build process without overcomplicating the example.

To build and install the module, you can then type:

$ python setup.py install

This should produce a lot of compiler spew.  If there are no errors, you should then be able to use the module.  Here is an example script that does just that:

demo.py
from ExtensionExample import *

print 'calling get_string(), should print out a random string'
print get_string()

print 'calling get_string(), should print out a random string'
print get_string()

print 'calling are_values_equal( 10, 10 ), should print True'
print are_values_equal( 10, 10 )

print 'calling are_values_equal( 10, 1 ), should print False'
print are_values_equal( 10, 1 )

print 'calling num_arguments( True ), should print 1'
print num_arguments( True )

print 'calling num_arguments( True, True ), should print 2'
print num_arguments( True, True )

print 'calling num_arguments( True, True, True ), should print 3'
print num_arguments( True, True, True )

print 'running: a = wrapped_class()'
a = wrapped_class()
print 'running: a.get_value(), should print -1'
print a.get_value()
print 'running: a.set_value( 5 )'
a.set_value( 5 )
print 'running: a.get_value(), should print 5'
print a.get_value()


If you run this script you should see the expected output from the script (evident from the script itself). You can download the entire example as a zip file here which includes all the source/header files as well as the build scripts and a readme.  Hopefully it will be helpful to other people, or just me later on.

Sunday, May 27, 2012

Example code for building CGAL::Polyhedron_3

I've found myself redoing this code relatively frequently, so I thought I would post it.  Nothing special, just a snippet that loads the file input.obj, builds a CGAL::Polyhedron_3 from it, and writes it out as dump.off. The code also includes a rudimentary Wavefront OBJ loader.


#include<fstream>
#include<vector>
#include<string>
#include<algorithm>

#include<CGAL/Simple_cartesian.h>
#include<CGAL/Polyhedron_incremental_builder_3.h>
#include<CGAL/Polyhedron_3.h>
#include<CGAL/IO/Polyhedron_iostream.h>

typedef CGAL::Simple_cartesian<double>     Kernel;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;
typedef Polyhedron::HalfedgeDS             HalfedgeDS;

// A modifier creating a triangle with the incremental builder.
template<class HDS>
class polyhedron_builder : public CGAL::Modifier_base<HDS> {
public:
 std::vector<double> &coords;
 std::vector<int>    &tris;
    polyhedron_builder( std::vector<double> &_coords, std::vector<int> &_tris ) : coords(_coords), tris(_tris) {}
    void operator()( HDS& hds) {
  typedef typename HDS::Vertex   Vertex;
        typedef typename Vertex::Point Point;

  // create a cgal incremental builder
        CGAL::Polyhedron_incremental_builder_3<HDS> B( hds, true);
        B.begin_surface( coords.size()/3, tris.size()/3 );
  
  // add the polyhedron vertices
  for( int i=0; i<(int)coords.size(); i+=3 ){
   B.add_vertex( Point( coords[i+0], coords[i+1], coords[i+2] ) );
  }
  
  // add the polyhedron triangles
  for( int i=0; i<(int)tris.size(); i+=3 ){
   B.begin_facet();
   B.add_vertex_to_facet( tris[i+0] );
   B.add_vertex_to_facet( tris[i+1] );
   B.add_vertex_to_facet( tris[i+2] );
   B.end_facet();
  }
  
  // finish up the surface
        B.end_surface();
    }
};

// reads the first integer from a string in the form
// "334/455/234" by stripping forward slashes and
// scanning the result
int get_first_integer( const char *v ){
 int ival;
 std::string s( v );
 std::replace( s.begin(), s.end(), '/', ' ' );
 sscanf( s.c_str(), "%d", &ival );
 return ival;
}

// barebones .OFF file reader, throws away texture coordinates, normals, etc.
// stores results in input coords array, packed [x0,y0,z0,x1,y1,z1,...] and
// tris array packed [T0a,T0b,T0c,T1a,T1b,T1c,...]
void load_obj( const char *filename, std::vector<double> &coords, std::vector<int> &tris ){
 double x, y, z;
 char line[1024], v0[1024], v1[1024], v2[1024];
 
 // open the file, return if open fails
 FILE *fp = fopen(filename, "r" );
 if( !fp ) return;
 
 // read lines from the file, if the first character of the
 // line is 'v', we are reading a vertex, otherwise, if the
 // first character is a 'f' we are reading a facet
 while( fgets( line, 1024, fp ) ){
  if( line[0] == 'v' ){
   sscanf( line, "%*s%lf%lf%lf", &x, &y, &z );
   coords.push_back( x );
   coords.push_back( y );
   coords.push_back( z );
  } else if( line[0] == 'f' ){
   sscanf( line, "%*s%s%s%s", v0, v1, v2 );
   tris.push_back( get_first_integer( v0 )-1 );
   tris.push_back( get_first_integer( v1 )-1 );
   tris.push_back( get_first_integer( v2 )-1 );
  }
 }
 fclose(fp); 
}

int main() {
 // two vectors to hold point coordinates and
 // triangle vertex indices
 std::vector<double> coords;
 std::vector<int>    tris;
 
 // load the input file
 load_obj( "input.obj", coords, tris );
 if( coords.size() == 0 )
  return 1;
 
 // build a polyhedron from the loaded arrays
 Polyhedron P;
 polyhedron_builder<HalfedgeDS> builder( coords, tris );
 P.delegate( builder );
 
 // write the polyhedron out as a .OFF file
 std::ofstream os("dump.off");
 os << P;
 os.close();
 
    return 0;
}

Thursday, April 19, 2012

Simple Command-Line Argument Handling

Much of the code I write is run from the command-line and often has a large number of input parameters.  I quickly got tired of handling command-line arguments and of having to remember the ordering of parameters.

The following code makes this easier by allowing parameters to be specified as key-value pairs.  You then call command_line_get_argument(.), with the desired argument name and type.  It searches through the list of arguments supplied to the program and attempts to read them into the correct type.  If the argument doesn't exist or doesn't convert to the correct type, the function returns false.

e.g.
// some variables declared here...
bool ok = true;
ok &= command_line_get_argument( argc, argv, "-input_file",      ARGUMENT_STRING, input_filename );
ok &= command_line_get_argument( argc, argv, "-output_basename", ARGUMENT_STRING, output_basename );
ok &= command_line_get_argument( argc, argv, "-num_iterations",  ARGUMENT_INT,    &num_iterations );
ok &= command_line_get_argument( argc, argv, "-prior_weight",    ARGUMENT_FLOAT,  &prior_weight );
if( !ok ){
 std::cout << "Usage:" << argv[0] << " [arguments]" << std::endl;
 std::cout << "arguments:" << std::endl;
 std::cout << "\t-input_file       [filename]         REQUIRED - input ground-truth image" << std::endl;
 std::cout << "\t-output_basename  [partial filename] REQUIRED - base name for output files, e.g. 'test01/image0_'" << std::endl;
 std::cout << "\t-num_iterations   [integer]          REQUIRED - number of iterations to perform, each is width*height mutations" << std::endl;
 std::cout << "\t-prior_weight     [real]             REQUIRED - weight for regularizer" << std::endl;
 return 1;
}

The 20 or so lines above replace what was initially closer to one hundred with error checking and ends up being much more readable. The programs are also nicer to work with, since they don't rely on a specific argument order.

The code for this can be downloaded from https://sites.google.com/site/jamesgregson/tmp/command_line_arguments.h. It can be used for whatever you'd like, provided it's attributed as mine.

Sunday, November 27, 2011

Got my face did!

I recently got to visit a local studio and while I was there they asked if I wanted my face 3D scanned.  Of course I did. Who wouldn't? Well at least one person.  But I did.  And here's the result, after runing it through the quadric-collapse decimator of Meshlab.



It's pretty damn good quality wise, even on challenging areas like my beard.  Now I just have to merge the scans and create a watertight mesh, ready for 3D printing!