## Pencil milling of images

Friday, December 30th, 2011 at 2:27 pm

The big idea is to use this 3D printer makerbot as a 2D plotter (by the application of a rubber band and a felt tip pen) to draw a representation of someone’s face on a post-it note.

We have the printer driver technology, we have a camera. All that is in between is software. So it should be easy, right?

Computer vision is a field that has been worked over for decades. Unfortunately it hasn’t resulted in a straightforward bitmap to vector drawing function I can call from Python. What I do have, however, are a lot of specialized CAM algorithms that I can apply to the problem. Specifically, the Pencil Milling algorithm. (That’s the one where you run a blunt object along all the surface to surface corners as if in the gutter of a kerb.)

First off, we use the FIND_EDGES filter which converts a bitmap into another bitmap with the sudden changes in brightness highlighted.

(The code which follows is as much for my reference — once I forget exactly how I did this — as much as it could be for anyone else’s benefit.)

```import Image, PIL.ImageFilter, PIL.ImageOps
i = Image.open("a_dscf1150.jpg")
i = PIL.ImageOps.grayscale(i)
j = i.filter(PIL.ImageFilter.FIND_EDGES)
j = j.filter(PIL.ImageFilter.GaussianBlur)
writestl("test0.stl", j)
j.save("test0.jpg")
```

This is the edge detection before it was blurred.

And here’s what it looks like after it is blurred:

The blurring is so I can produce a 3D triangulated surface without too much noise. The function writestl() above looks like this:

```import struct, math

def pushtriangle(fout, x0, y0, z0, x1, y1, z1, x2, y2, z2):
v1x,v1y,v1z = x1-x0, y1-y0, z1-z0
v2x,v2y,v2z = x2-x0, y2-y0, z2-z0
nx,ny,nz = (v1y*v2z - v1z*v2y), -(v1x*v2z - v1z*v2x), (v1x*v2y - v1y*v2x)
nv = math.sqrt(nx*nx+ny*ny+nz*nz) or 1
n0,n1,n2 = nx/nv, ny/nv, nz/nv
fout.write(struct.pack("<12f2c", n0, n1, n2, x0, y0, z0, x1, y1, z1, x2, y2, z2, " ", " "))

def writestl(fname, j):
w, h = j.size
fout = open(fname, "wb")
fout.write(struct.pack("<i", 2*(h-1)*(w-1)))
zfac = -40.0/255.0
hfac = 1.0
for y in xrange(1, h):
for x in xrange(1, w):
x00, y00, z00 = (x-1)*hfac, (y-1)*hfac, pix[x-1, y-1]*zfac
x01, y01, z01 = (x-1)*hfac, (y)*hfac, pix[x-1, y]*zfac
x10, y10, z10 = (x)*hfac, (y-1)*hfac, pix[x, y-1]*zfac
x11, y11, z11 = (x)*hfac, (y)*hfac, pix[x, y]*zfac
pushtriangle(fout, x00, y00, z00, x01, y01, z01, x11, y11, z11)
pushtriangle(fout, x00, y00, z00, x11, y11, z11, x10, y10, z10)
fout.close()
```

This creates a regular grid from the pixels of the image and outputs the triangulated surface (as a binary STL file) which we can use as input for the CAM pencil algorithm.

Here’s how that surface looks. It has been reflected in the Y-axis because bitmaps put their origin at the top left hand corner rather than the bottom left.

There are many ways to explain the pencil milling toolpath. In this case you can imagine dropping a ball-bearing onto the surface and rolling it around until it falls into one of the grooves, and then tilting the surface so it rolls along the groove.

Accessing this code is reasonably complicated, but it builds on the earlier work of running these algorithms from minimal Python code. There is a sophisticated system in there struggling to get out. But as with all developing systems, the design can get messy when only one application ever uses it. A direct unit testing process would be a second application.

First step is to load the triangulated surface (see that earlier work) and define a “Toolsurface”, which is just a triangulated surface together with a particular cutter shape. The reason for the existence of this class is to hold the hefty pre-calculations you can have on the triangles once you know a specific tool shape (generally offsetting them in space along their planar normal). The toolshape definition also looks weird, putting the ball radius into the function 5 times like that. There’s probably a good reason, which I don’t have time to go into right now.

```fssurf = LoadSurface("test0.stl")
fstoolsurf = FsToolSurface.New()
r = 4
fstoolsurf.SetLiftCutterShape(-20, 0, r, r, 0, r, r, 100, r)
```

The next step is to make the machining boundary. I know there is a function for making a rectangular boundary in one call, but I like to do it from first principles. A boundary is a set of paths (because there can be islands). The Build function is the same as the one on FsSurface — it puts things into boxes so we can do particular function calls very fast (eg measure the distance to the boundary from a given point if it is less than 1mm away).

```ulo, uhi, vlo, vhi = -50, 550, -50, 550
fsboundary = FsBoundaries.New()
fspath = FsPath2X.New(0.0)
fspath.Start(ulo, vlo)
fsboundary = FsBoundaries.New()
fsboundary.Build(5.0)
```

Now we can make the implicit area, which is the object that defines how the pencil milling algorithm will operate, and is based on the standard boundary detection algorithm that is used for slices, except it is a very special case where there is no inside because the boundaries are actually doubled up and down the channels. There’s no way I can explain it very well without a whole lot more diagrams, but it is geometrically natural. Once again there are numerous hidden parameters controlling how well the contours are followed and the tolerance. For example, maxHdiffContour is the largest segment length as we track the contour, deltaHdiffContour is the shortest segment length, and minCNdotContour is the cosine of the largest angle change allowed between two segments (as long as they are both greater than deltaHdiffContour in length). In terms of the pencil milling, minCNdotFibrePencil is the cosine of the bitangency angle.

```fsia = FsImplicitAreaSurface.New(FsImplicitAreaSurface.IAT_PENCIL)
fsia.SetZLow(-100)
fsia.SetZHigh(100)
fsia.SetMachiningBoundaries(fsboundary, 0)
#SetPencilConditions(deltaHdiffFibrePencil, maxHdiffFibrePencil, maxZdiffFibrePencil,
minCNdotFibrePencil, marginZCdiffFibrePencil, marginCNZbothFibrePencil, False)
fsia.SetPencilConditions(0.002, 2, 2, 0.93, 1, 1, False)
#SetContourConditions(minCNdotContour, maxZdiffContour, deltaHdiffContour,
maxHdiffContour, maxCPcuspContour, minBNdotContour)
fsia.SetContourConditions(0.9, -1.0, 0.002, 2, -1.0, 0.9)
```

I can see a good deal more documentation and rationalization needs to go on, as well as some interesting systems for isolation testing to get an idea of how much things change when the settings are changed on experimental parts. Often the tolerance requested by the user implies one sample rate, but the input geometry requires operating at a smaller sample rate so it doesn’t accidentally skip over a feature, which means that they get a tighter tolerance than they wanted. (Tighter tolerances are not desirable as they take longer to calculate. I wonder if every algorithm should return the tolerance it actually thinks it got?)

Now we make the weave, which is the boundary model derived from the implicit area.

```fsweave = FsWeave.New()
fsweave.SetShape(ulo-10, uhi+10, vlo-10, vhi+10, 2.1)
fsia.GenWeavePencilMainFibres(fsweave)
```

Remember I said that pencil contours were special because they are doubled up? It requires special functions to pull them out:

```npp = fsia.PencilPathContours(fsweave, 5.0, 4.0, 0.0)
pths = [ ]
for ipp in range(npp):
fspath = FsPath2X.New(0)
RecordPencilContour(fspath, fsweave, ipp, 0.02)
print "path", ipp, "with points", fspath.GetNpts()
pths.append([ (fspath.GetD(i, 0), fspath.GetD(i, 1), fspath.GetD(i, 2))
for i in range(fspath.GetNpts()) ])
```

I can’t quite run everything under the same version of Python (the PyGL stuff isn’t installed in the main CAM version), so I transferred this path geometry back to the user interface using pickle.dump() as a quick and dirty method (to avoid inventing yet another geometric file format). Here is the code to draw the paths into a display list:

```genList = GL.glGenLists(1)
GL.glNewList(genList, GL.GL_COMPILE)
for pth in pths:
GL.glBegin(GL.GL_LINE_STRIP)
for px, py, pz in pth:
GL.glVertex3d(px, py, pz)
GL.glEnd()
GL.glEndList()
```

And, after all that, here is the result, which has a sort of skin cell network membrane network effect.

Not very good.

I went back to the original greyscale image and simply converted that into a height relief, which looks kind of weird and distorted when you see it from the side.

But what happens with the pencil milling passes hunting out the grooves in this?

Really quite a lot better. No idea what happened to the coffee cup though.

What the pencil milling algorithm is really about is detecting and tracking discontinuities in a vector field. In the CAM case the vectors are the contact angles between a ball nosed cutter and the surface at a given XY axis location. But if there’s a filter which derives a vector directly from the image, there’s no reason the algorithm shouldn’t be applied (aside from having to make new interfaces and so forth).

Or maybe I’ll find the bitmap to contours function I need in python.

Having gone through this all, what I know I definitely need now is a method for serializing a FsWeave into a binary object. All the other entities can be created by calling simple functions which set their parameters, but the centrality and opaqueness of the FsWeave object is what is getting in the way of producing a unit test situation.

• 1. anders replies at 3rd January 2012, 9:40 am :

Hi Julian,

Does the pencil-milling path generation rely on the low-level drop-cutter function?

If so, for a single drop-cutter call can you actually make contact with two triangles at once? And detect the cusps where we want pencil-paths this way?

OR, do you call drop cutter at position x and get one contact-normal, then call it again at a slightly moved position x+dx and get another contact-normal? If these two normals differ from each other sufficiently then we have a discontinuity and should generate a pencil-path somewhere between x and x+dx ?
How do you differentiate between convex “hills” and concave “cusps” with this strategy?

thanks,
Anders

• 2. JT replies at 3rd January 2012, 10:42 am :

The latter case (sampling at x and x+dx) is the way to detect the double contact points.

Actually what it does is looks for large changes in the contact normal between evenly spaced sample points, and then subdivides the sample rate until it has two samples 0.002mm apart where there is still a contact normal change, and marks that as a discontinuity.

It is only able to detect concave cusps. However there you can pick up convex (upward) sharp edges (which you don’t want) if you are doing this with a flat bottomed cutter.