Treating Image Edges

There are some good examples floating around of creating customized image processing kernels in the python imaging library (PIL). It's very convenient to simply define your own mask, and the library takes care of the convolution.

One drawback is that the edges don't get filtered in this process. There are various ad-hoc things you can do. You can just drop a line of pixels all the way around the picture, so the gradient image of a 800x600 picture would be 798x598 pixels (if you're using a compact stencil), or you could just copy the edge pixels from the old picture to the new picture. So in effect, the edges don't get filtered. For most applications this is probably no big deal, but it might be, and we can do better.

All we have to do is change our difference operator for the edge pixels, use a forward or backward difference on the edges rather than the central differences used on the interior. Then you get gradient info for your edge pixels too!

Customizing the PIL kernel class has been done, so I'm going to use Fortran90 and F2Py for my difference operators. I'll still use the PIL for all of the convenient image reading/writing utilities though. No need to implement that in Fortran!

We can use our little Maxima script to generate the differences we need on the edges. If you've taken a basic numerical methods course you'll probably recognize them as the three point forward and backward difference approximations for the first derivative.



The Python to do just that is here:

import sys, math, Image, numpy
from filters import *


if __name__ == '__main__':
in_img = Image.open(sys.argv[1])
in_img = in_img.convert("L")
out_img = Image.new("L", in_img.size)
pix = in_img.getdata()
# put the pixels into a fortran contiguous array:
pix_array = numpy.array(pix)
pix_array = pix_array.reshape( in_img.size, order='F' )
# calculate the gradient
g = first_order_gradient(pix_array,in_img.size[0],in_img.size[1])
gmag = 255 - numpy.sqrt( g[0,:,:]*g[0,:,:] + g[1,:,:]*g[1,:,:] )
gmag.astype("int")
# pack the gradient magnitude into a new image
out_img.putdata( gmag.reshape( gmag.size, order='F' ) )
out_img.save(sys.argv[2])

Where I've used F2Py to generate the filters module from a Fortran90 file. The module defines the subroutine first_order_gradient(), which implements the ideas discussed above. The gradients on the interior points are calculated using a central difference:

! apply the central operator to the interior pixels
do j = 2, nj-1
do i = 2, ni-1
g(1,i,j) = (in_image(i+1,j) - in_image(i-1,j)) / 2.0
g(2,i,j) = (in_image(i,j+1) - in_image(i,j-1)) / 2.0
end do
end do

The gradients along the edges are calculated using the forward and backward differences (only x edges shown for brevity):

! loop over the y dimension and apply the asymmetric operators to the
! x min and x max edges
do j = 2, nj-1 ! leave out the corners
i = 1 ! minimum edge, go forward
g(1,i,j) = (-3.0*in_image(i,j) + 4.0*in_image(i+1,j) - in_image(i+2,j))/2.0
i = ni ! maximum edge, go backward
g(1,i,j) = (3.0*in_image(i,j) - 4.0*in_image(i-1,j) + in_image(i-2,j))/2.0
! no change for the y derivative
i = 1
g(2,i,j) = (in_image(i,j+1) - in_image(i,j-1)) / 2.0
i = ni
g(2,i,j) = (in_image(i,j+1) - in_image(i,j-1)) / 2.0
end do

Finally the four corners are calculated separately, because each one has a unique combination of stencils required. Putting the loops in Fortran speeds things up considerably compared to the pure Python implementations I did before. Here is the gradient magnitude of the dials from Poor Man's ADC calculated with the above technique:

Sobel, Scharr and Truncation Error

In image processing the discrete difference operators to calculate gradients and second derivatives are called "masks". They are often referred to as "stencils" in other applications of finite differences for solving PDEs, estimating sensitivities, etc.

Some of the kernels (a mask plus particular values) for calculating gradients are Sobel:

[ 1 2 1]
Gy = [ 0 0 0]
[-1 -2 -1]

[ 1 0 -1]
Gx = [ 2 0 -2]
[ 1 0 -1]

or Sharr:

[ 3 10 3]
Gy = [ 0 0 0]
[-3 -10 -3]

[ 3 0 -3]
Gx = [ 10 0 -10]
[ 3 0 -3]

or the simpler stencil we used in the previous posts, which only used neighbouring pixels in one direction to estimate the gradient. What if we want to design our own masks? A useful tool is the multidimensional Taylor series:

In two dimensions we have:

which includes up to third order terms. We can't actually include the two pure third order terms because we only have three points in either the x or y direction, and we need a minimum of four. We can include the mixed third order terms though, so that gives us 8 constraints on our coefficients. Instead of including all 9 points in the stencil we have to drop one, so in the Maxima code below you'll notice that the central point is commented out.

Lets examine the truncation error for the two masks above. First we define a vector of delta x and delta y corresponding to each point in the stencil:

/* define the points in the stencil or mask: */
steps : matrix( /* dx, dy */
[-1, 1],
[0, 1],
[1, 1],
[-1, 0],
/* [0,0], */
[1, 0],
[-1, -1],
[0,-1],
[1, -1] );

Then vectors describing the coefficients in the two masks (we'll just use the y-component for this example):

/* coefficients in sobel mask */
sobel_y : matrix(
[1],
[2],
[1],
[0],
/*[0],*/
[0],
[-1],
[-2],
[-1]
);
/* coefficients in sharr mask */
sharr_y : matrix(
[3],
[10],
[3],
[0],
/*[0],*/
[0],
[-3],
[-10],
[-3]
);

To evaluate the truncation error we also need to create an operator that is based on the multidimensional Taylor series shown above. This will let us simply multiply those vectors by the operator to see what the masks are actually estimating. Here's the operator:

n : matrix_size( steps );
A : zeromatrix(n[1],n[1]);

for j : 1 thru n[1] do ( A[1,j] : 1 );/* f */
for j : 1 thru n[1] do ( A[2,j] : steps[j,1] );/* df/dx */
for j : 1 thru n[1] do ( A[3,j] : steps[j,2] );/* df/dy */
for j : 1 thru n[1] do ( A[4,j] : steps[j,1]^2/2! );/* d2f/dx2 */
for j : 1 thru n[1] do ( A[5,j] : steps[j,2]^2/2! );/* d2f/dy2 */
for j : 1 thru n[1] do ( A[6,j] : 2*steps[j,1]*steps[j,2]/2! );/* d2f/dxdy */
for j : 1 thru n[1] do ( A[7,j] : 3*steps[j,1]*steps[j,2]^2/3! );/* d3f/dxdy2 */
for j : 1 thru n[1] do ( A[8,j] : 3*steps[j,1]^2*steps[j,2]/3! );/* d3f/dx2dy */

Now the truncation error terms are shown by a couple of matrix vector multiplies:

sobel_te : A.sobel_y;
sharr_te : A.sharr_y;

Rather than just getting a single non-zero entry in the third place of the resulting vector we get a non-zero entry in the last place as well. Looking back at our operator, this is the row for the third order mixed derivative term, (d^3 f)/(dx^2 dy). In other words both the Sobel and Sharr operators introduce an error in their estimation of the gradient that is proportional to this mixed derivative.

How to avoid this error? Well, we can choose a right hand side vector that only has a 1 in the third position and zeros everywhere else, and then invert our operator to solve for what the values in our mask should be.

/* define which derivative we want to estimate, constrain the unwanted
ones to be zero: */
derivs : matrix(
[0], /* f */
[0], /* df/dx */
[1], /* df/dy */
[0], /* d2f/dx2 */
[0], /* d2f/dy2 */
[0], /* d2f/dxdy */
[0], /* d3f/dx2dy */
[0] /* d3f/dxdy2 */
);
inv_A : invert(A);
coeffs : inv_A.derivs;

This gives a result with only non-zero coefficients for the points at [0,1] and [0,-1], in other words, the simple estimate we've already been using! There's really no need to include all those other points for the gradient calculation.

Open Source Tools

I just wanted to put in another plug for Maxima, the open source computer algebra system that I use regularly. The great thing about using an open source tool, other than customizability and the ability to "look under the hood", is the responsiveness of the user and developer community.

I found a little bug in the f90() function that caused it to use the old fixed-format FORTRAN77 style line continuation characters when outputting a long matrix expression (bug report) rather than the correct Fortran90 free-form style. I posted a bug report to the list on the 18th of October, I had a reply from another user who understood the problem that same day. A work-around was posted (here by me), and the fix was committed to the code on the 24th of November by one of the project's developers. I had the fixed version of the code shortly after that without any effort on my part because I do automatic updates with yum.

Talk about sweet, if I had sent in a bug report like this to a commercial software company I certainly wouldn't have seen such a quick response, a little nuisance like this with a fairly easy work-around wouldn't rate very high on anyone's priority list (but mine), and they probably would have made me pay for a new version to get the fix when it finally did happen! There's a lot of value in such responsive support that makes up for a lack of some of the fancier (but usually unnecessary) features available in some of the commercial computer algebra packages.

Implementing Numerical Methods: an IDE (sort of)

What would an integrated development environment for numerical algorithms look like? To me it looks a lot like Fedora. That’s right, Fedora Linux is your integrated development environment. Oh, separate programs in an operating system isn't integrated enough for you? They ought to be, Python or Bash are right there waiting to glue it all together.

An implementer or developer of numerical methods ought to have a good computer algebra system(CAS) available. This is where the high-level algorithm design needs to take place. All the rest of the low level coding needs to be automated as much as possible. You as the human in the loop need to spend time thinking about the bigger picture and get Maxima to code Fortran for you.

Numerical software needs testing. Lots of rigorous testing. You can use your (CAS) to generate unit tests by evaluating the expressions you converted to Fortran for specific cases, and then converting these answers to Fortran too. Then you can easily wrap those routines in a unit testing framework using F2Py.

This approach lets you automate the generation of the low-level functions that will be in your inner loops (which need to be in Fortran for performance sake). Then use Python's powerful standard library and additional packages for things like parsing input decks, file reading and creation, setting up parallel jobs and generating graphics. All things which can be a real pain in Fortran. The real benefit is that the human effort is kept focused where it is most profitably employed: at the higher levels of abstraction of the mathematical concepts rather than the low level coding.

There it is, a complete "environment" in three building blocks: Maxima, Fortran and Python. Almost completely automated from governing equations in the CAS to low-level number crunching in Fortran. Some human intervention required to bring it all together into a working program, but that’ll be in Python, which is easy.

This is the basic method used in the second half of the magnetron problem. I'll post some more simple examples of the whole work-flow in use after the holidays.

Stupid Human Tricks in Peer Review

I felt like peer review is under attack after reading the second Slashdot article in as many days (1 and 2) about embarrassing failures in peer review of scientific journals/conference articles. The Sokal affair is a pretty well publicized hoax along these lines. The perpetrator wanted to prove a point that articles on post-modern cultural studies were indistinguishable from gobbledygook, surprising no-one outside the field, he was right!

It seems like the above two recent occurrences demonstrate that peer review is in need of an overhaul (or at least a minor upgrade). The second story is about a paper generated by a context free grammar being accepted to a conference. It seems like therein lies a clue to a possible solution.

We need some way to score the "goodness" of a reviewer's work. It's intractable to focus on creating an algorithm for finding "good" original research papers, that basically amounts to implementing some sort of beyond the state-of-the-art artificial intelligence.

What about focusing on "badness"? What if a certain percentage of the articles sent to a reviewer are generated by a context free grammar trained on the archives of relevant journals. These are now known "bad" articles. If the reviewer passes these papers on to publication or presentation then his reputation/score as a reviewer drops. It would need to be a double-blind process. Neither the reviewer nor the institution sending articles for review should know which are true papers and which are generated.

This would be a simple way to check and see if reviewers were paying attention, though it would increase the work-load on a group of folks who are not reimbursed for their time reviewing articles. It would probably improve the results of peer review though, especially if the fact that someone passed along a bogus paper were publicized in the journal itself, right after that page in the front that lists all of the editors and whatnot.

More on Speeding Up Octave

In a previous post we found that speeding up Octave amounted mostly to avoiding for loops. The simplest way to do that is to operate directly on vectors using the built-in operators and functions (which are fast because they are compiled) or cast your problem as a sparse matrix solve. This second option is especially helpful when subsequent calculations depend on previous iterations (the code can't really be vectorized). This option is a very general way of avoiding loops in those cases.

The simple example given previously was a numerical integration of x-squared:

p1(1) = 0;
for( i=2:N )
t = t + dt;
p1(i) = p1( i-1 ) + dt*2*t;
endfor

The way to recast the problem as a sparse matrix solve is to think of p1 as the vector of unknowns, and each iteration of the loop as an equation in the system we want to solve. Writing down the system gives us:


The important detail to remember is to use the functions in Octave to allocate the sparse matrix, or you could easily find yourself writing more really slow for-loops just to create the sparse matrix which you hoped would save lots of time by avoiding for-loops. Talk about the long, slow way around!

Two very useful functions are speye() and spdiag(). The first returns a sparse identity matrix, which is often a good initial building block for many useful operators. The second allows you to place vectors (allocated quickly with the usual vector suspects such as ones() and zeros() ) along the diagonals of the sparse matrix.

# create the main diagonal
A = speye( N );
# alternatively could use spdiag:
# A = spdiag( ones(N,1), 0 );
A = A + spdiag( -ones(N-1,1), -1 ); # add the first sub-diagonal

If your operator isn't banded then you'll need to use spconvert(), which takes as its argument a three (or four if you need a complex result) column matrix. Each row of the argument defines a non-zero entry in the sparse matrix. The first column is the row index, the second column is the column index and the third column is the entry value (fourth column being the imaginary part of the value, if necessary).

Using those three functions should allow you to allocate a sparse matrix in Octave without resorting to for loops (which was why we embarked on this journey to begin with).

This method is counter-intuitive to folks who come from a Fortran (or other compiled language) background, because writing down the loop is the simple, efficient way to solve the problem. It also seems like 'wasting' memory to store all of those redundant coefficients. The timing results speak for themselves though, if you want to stay completely in Octave (or Matlab or Python) sometimes you have to do weird things to get reasonable performance. Of course, if speed really becomes a problem then the inner loops of your calculations need to move to a compiled language that can then be called from Octave, this is a bit more complicated than our little sparse matrix method.

More Edge Detection

I can't really decide on the best way to treat the different channels (red, green, blue) in a color image when doing edge detection. Is it better to just average the magnitude of the gradient in each channel? Treat the channels separately and find red edges, green edges and blue edges? My pragmatic (if ad hoc) choice is to just add them all up, and use a cludge to keep from getting negative gray-scale values.

Here's the Python for gradient-based edges. The magnitude of each channel's gradient is subtracted from 255. Since gray-scale values range from 0 to 255 we use a max() function to set any resulting negative values to 0. The edges detected using this scheme are shown below, on the right.

inpix = img.load()
s = 1.0;
# calculate the gradient in all three channels, sum of squares for mag
for y in xrange(2,h-1):
for x in xrange(2,w-1):
rdx = ( inpix[x+1,y][0] - inpix[x-1,y][0] )/2
gdx = ( inpix[x+1,y][1] - inpix[x-1,y][1] )/2
bdx = ( inpix[x+1,y][2] - inpix[x-1,y][2] )/2
rdy = ( inpix[x,y+1][0] - inpix[x,y-1][0] )/2
gdy = ( inpix[x,y+1][1] - inpix[x,y-1][1] )/2
bdy = ( inpix[x,y+1][2] - inpix[x,y-1][2] )/2
outpix[x,y] = max( 0, 255 - s*( math.sqrt( rdx**2 + rdy**2 ) + math.sqrt( gdx**2 + gdy**2 ) + math.sqrt( bdx**2 + bdy**2 ) ) )

outimg.save(sys.argv[2],"JPEG")

Taking a careful look at the code, you may notice the parameter s floating around in there. This lets you "darken" the edges by choosing something larger than 1. Here's the same image with s=3.0 (on the left):

Like in the previous image processing posts, these examples were made using the Python Imaging Library for image input/output and easy access to the pixels.

My Own Little Pareto Distribution

I signed up for Google Analytics, which collects bunches of neat data. One of the interesting reports is on the number of views each page on your website gets. I plotted them (in Octave of course) with pages across the x-axis, number of views on the y-axis, it looks something like a power law distribution.

x = [11,6,4,2,1,1,1,1,1,1];
bar(x);
print("-dpng","pageviews.png");

Judging by the numbers of pageviews I fear my website is in the trivial many, but that's ok a few people found the octave posts useful.

Asymmetric Complex Step

To solve Laplace's equation with a complex step finite difference we needed to develop asymmetric difference relations to solve on the +/- imaginary planes.

On the plane in the negative imaginary direction we difference centrally in the real direction and forward two steps in the imaginary direction:

And on the plane in the positive imaginary direction we difference centrally in the real direction and backward two steps in the imaginary direction:

We still want the real part of the answer as we did before, but now we have imaginary coefficients for our function values in the real direction. While we have a compact stencil on the real plane (only nearest neighbors are used), we lose that property on either imaginary plane of the solution. This means our system matrix bandwidth is larger than it would be if we had a purely compact stencil for all parts of the solution.

It's fairly easy to calculate arbitrary derivative approximation formulas using Maxima.

/* calculate the coefficients for a finite difference approximation of
a derivative, can include complex steps */

steps : matrix([1,2,%i,2*%i]); /* multiple of the step size h */

/* sum df d2f d3f d4f */
derivs : matrix([0], [1], [0], [0], [0]); /* the rhs for our system */

n : matrix_size( steps );

A : zeromatrix(n[2]+1,n[2]+1);

for j : 1 thru n[2]+1 do
( A[1,j] : 1,
for i : 2 thru n[2]+1 do
(
if j = 1 then A[i,j] : 0 else A[i,j] : steps[1,j-1]^(i-1)/(i-1)!
)
);

invA : invert(A);

coeffs : expand( invA.derivs );
expand( 20*coeffs );

Where the vector steps tells which points will be included in the stencil, and the vector derivs tells which derivative we want to approximate. This example is calculating a real and complex forward difference approximation for the first derivative.

You can download this Maxima batch file here. This example can be run from the command line:

$ maxima -b fd.mac

or using the batch command from an interactive Maxima session:

(%i1) batch("fd.mac")

Latex for Blogging

I've been using htlatex command that comes with the tex4ht package to generate images of nicely typeset equations for use in my blog posts. There are several alternatives for turning latex into html. This article has a decent overview of a few of the methods you might use if you want to have a single tex source file for publishing on the web and on paper. I'm not that serious about my content management. I don't mind a little manual cutting and pasting.

I tried using latex2html, but that is really more suited for processing a large structured document and turning it into a large, structured set of web-pages. All I need is equations, so htlatex is sufficient. If I have a list of equations in 'example.tex', then the command:

$ htlatex example.tex

Produces png image files of all the equations right there in that directory, it even includes the equation labels and alternative text in the generated html, which latex2html does not do. I've noticed that the images generated for a 12pt document look a little better (because they are bigger) when you insert them into a post.

If you're a Fedora user, getting this useful little tool is as easy as

$ yum -y install tetex-tex4ht

Complex Step: Subtractive Cancellation

We showed in an earlier post that approximating first derivatives with complex steps leads to a second order accurate formula that avoids the subtractive cancellation associated with a normal finite difference.

Following the simple convergence study for second derivative approximations we did for cos(). We'll look at the convergence of some first derivative approximation formulas.

The simple complex step approximation is given by

The standard second order finite difference is given by

The 4th order combined complex and real step approximation is

Notice that this formula is not able to avoid subtractive cancellation the way the second order approximation does, but we get a fourth order accurate approximation on a compact stencil.

The convergence of these three formulas for approximating the derivative of cosine is shown in the figure below.

The 4th order scheme runs into subtractive cancellation earlier than the standard 2nd order scheme, but it remains more accurate over the examined range of resolution. The complex step approximation continues to converge as the step size decreases even after the other two formulas diverge.

So for doing sensitivity analysis the second order complex step formula should be used, because only one step needs to be taken, and arbitrarily high precision can be achieved. If you need to approximate a derivative over a region, then the higher order method will pay-off because you can get better accuracy with fewer points.