## Tuesday, September 18, 2012

### Simple printf style formatting of std::string

Nothing special here but a quick snippet that allows an std::string to be formatted similar to printf, which I prefer.  Not a safe function by any means, but great for debugging and log files. I must have looked up the C variable argument list code dozens of times over the years.  May not work under Visual C++ if Microsoft hasn't yet come to terms with vsprintf() actually being a part of the C standard library.

#include<string>
#include<cstdio>
#include<cstdarg>

#define STR_FORMAT_BUFFER_SIZE 2048

std::string str_format( const char *fmt, ... ){
char buffer[STR_FORMAT_BUFFER_SIZE];
va_list arg;
va_start( arg, fmt );
vsprintf( buffer, fmt, arg );
va_end(arg);
return std::string(buffer);
}


## Sunday, September 16, 2012

### Low Cost CNC Part V - A Redesign

I've been working on a homemade CNC now for quite some time.  My goal for the project was to produce something modular, where shaft-mounts, motor mounts and bearings were entirely separate parts using a standardized mounting pattern.  This would allow the machines to be put together like Legos and has a lot of advantages, like the ability to mix and match drive options, e.g. in the image below, one axis is screw-driver while the other is belt-driven.

I still think this concept has merit, however the part designs that I ended up using made for bulky and not particularly stiff machines. In the photo above, the top of the XY table is close to six inches from the base plate and the machine itself has considerable give.  Even worse, I didn't build enough slop into the designs to accommodate tolerances for manually built mounting plates, so it was actually quite difficult to get all the pieces to play nice.

I've since redesigned the machine to use more compact mounts, merged the shaft and motor supports into a single axis-end part and moved to half-inch shafting.  The switch to larger shafting results in a much stiffer machine, but unfortunately does require bushing style bearings.  Using the combined axis ends simplifies alignment, but unfortunately precludes belt drives.  I've also switched to Nema 17 steppers from Nema 23s, which leads to a more compact overall machine.  Surprisingly they don't seem to have much affect on the overall machine speed.  This makes for a much cleaner design:

Like the old parts, the new parts are 3D printed, but this time on my recently acquired RepRap.  The meshes were parametrically generated using my Python contructive solid geometry library, which is becoming quite usable.  This serves for prototypes, but final parts could be machines from plastic or aluminum.  The fixed-end for the leadscrew uses the same axis-end as the motor side, but with some bearing plates and Nylin nuts to take the axial loads:

After switching to acetal linear bushings, the bearings mounts can be made much more compact, reducing the height of the XY table from about six inches to a shade over 3.  These seem to run nicely on plain hardware store shafting, although hardened precision linear shafting would obviously be stiffer and smoother than cold-rolled bar stock.

Overall the machine is looking much cleaner.  I've decided to go with a knee-mill style machine rather than a gantry arrangement.  This means an XY table horizontally with a third axis mounted vertically. I bought a massive piece of aluminum channel to serve as the base of the machine, 10"x3"x20".  This provides a sturdy mount, and the controller and power-supply can be mounted to the underside.  Here's the thing midway through the build process yesterday:

And here's the finished-except-for-limit-switches XY table mounted to the channel.  The green duct tape is just there to stop it from gouging up my coffee table when I move it.  I will also probably replace the top cheese-plate with a thicker piece of flat-bar that's had all the holes drilled and tapped; I can't tell you how much I'm looking forward to doing that.

The mount for the Z-axis will be attached at the near-end of the machine.  I still have to design and build that mount, as well as print out the parts for the third axis.

After assembly, I simply had to try it out, missing limit switches or not.  I hooked it up to my 3-axis controller board with the Pololu A4988 carrier board carriers. The Arduino is flashed with GRBL, and wired to the carrier boards, which are in turn wired to the steppers.  I need to find some extra heatsinks for the X and Y axes, along with some thermal tape or epoxy, but they seem to run cool enough for testing anyway.

Here is it running, the first time I've gotten one of these CNC projects using actual GCode.  The leadscrew for the Y-axis is a bit loose in the fixed-end bearings, causing the knocking sounds, I have to look into the alignment, but otherwise it's working pretty well.  Speed is about one inch per second.

Next up will be designing a mount for the Z-axis.  I'm quite happy with the current set of printed parts and will probably continue to use them.

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