Implementing technical paper algorithms in C++ or MATLAB

I am an Electrical Engineering undergrad. I've reading many technical papers about Signal and Image processing algorithms (reconstruction, segmentation, filtering, etc). Most of the algorithms shown in these papers are defined over continuous time and continuous frequency, and often give the solutions in terms of complicated equations. How would you implement a technical paper from scratch in C++ or MATLAB in order to replicate the results obtained in said paper?

More specifically, I was looking at the paper "A general cone-beam reconstruction algorithm" by Wang et al (IEEE Trans Med Imaging. 1993;12(3):486-96), and I was wondering, how do I even start implementing their algorithm? Equation 10 gives you the formula of the reconstructed image at . How would you code that up? Would you have a for-loop going through each voxel and computing the corresponding formula? How would you code functions of functions in that formula? How would you evaluate functions at arbitrary points?

I've read the book "Digital Image Processing" by Gonzalez and Woods but I'm still at a loss. I have also read about the Numerical Recipes book series. Would that be the correct way?

What are your experiences programming algorithms from research papers? Any tips or suggestions?

• I'm gonna take a look at the paper when I have a chance. But I believe this is all about X Y Z points in a given graphic. You define a vertex and then work from there.
– Cassio
Aug 23 '11 at 23:12
• Typically, one discretizes the signals by sampling, and then convert integrals to sums. Aug 23 '11 at 23:54
• So I've read about sampling and converting the integrals to sums, but how do you evaluate the integrand at each sampling point if the functions in the integrand are stored as matrices?
– Damian
Aug 24 '11 at 0:39
• Damian, have you seen how the radon transform is inverted via backprojection? This is a slightly simpler example which I could explain if it would interest you. It is used for tomography using plane waves rather than the conical sampling described in the paper you posted. en.wikipedia.org/wiki/Radon_transform Aug 24 '11 at 0:57
• @mr-crt, is it possible to migrate to dsp.SE instead? Aug 24 '11 at 6:07

Signal processing algorithms defined in continuous time/space/frequency are typically implemented by sampling the signal on a discrete grid and converting integrals into sums (and derivatives into differences). Spatial filters are implemented through convolution with a convolution kernel (i.e. weighted sum of neighbors).

There is a huge body of knowledge concerning filtering sampled time-domain signals; time-domain filters are implemented as either finite impulse response filters, where the current output sample is computed as a weighted sum of the previous N input samples; or infinite impulse response filters, where the current output is a weighted sum of the previous inputs and previous outputs. Formally, the discrete time filters are described using the z-transform, which is the discrete-time analog to the Laplace transform. The bilinear transform maps one to the other (c2d and d2c in Matlab).

How would you evaluate functions at arbitrary points?

When you need the value of a signal at a point that does not lie directly on your sampling grid, you interpolate its value from nearby points. Interpolation can be as simple as choosing the nearest sample, computing a weighted average of the nearest samples, or fitting an arbitrarily complicated analytic function to the sampled data and evaluating this function at the needed coordinates. Interpolating onto a uniform finer grid is upsampling. If your original (continuous) signal does not contain details (i.e. frequencies) finer than half the sampling grid, then the continuous function can be perfectly reconstructed from the sampled version (the Nyquist-Shannon sampling theorem). For an example of how you can interpolate in 2D, see bilinear interpolation.

In Matlab you can use interp1 or interp2 to interpolate 1D or regularly sampled 2D data (respectively), or griddata to interpolate from irregularly sampled 2D data.

Would you have a for-loop going through each voxel and computing the corresponding formula?

Yes, exactly.

Matlab saves you from having to do this via explicit for-loops because it is designed to operate on matrices and vectors (i.e. multidimensional arrays). In Matlab this is called "vectorization". Definite integrals can be approximated with sum, cumsum, trapz, cumtrapz, etc.

I've read the book "Digital Image Processing" by Gonzalez and Woods but I'm still at a loss. I have also read about the Numerical Recipes book series. Would that be the correct way?

Yes, Numerical Recipes would be a great start. It is very practical and covers most of the numerical methods you'll end up needing. (You will find that Matlab already implements everything you need, but Numerical Recipes will provide excellent background.)

I have taken an "algorithms and data structures" class, but I don't see the relation between the material presented there and implementing scientific algorithms.

The material treated in "Algorithms and data structures" courses tends to concentrate on structures like lists, arrays, trees, and graphs containing integers or strings and operations like sorting and selecting: problems for which there is typically a single correct result. When it comes to scientific algorithms, this is only half of the story. The other half concerns methods to estimate real numbers and analytic functions. You'll find this in a course on "Numerical Methods" (or "Numerical Analysis"; like this one - scroll down for the slides): how to estimate special functions, how to estimate integrals and derivatives, etc. Here one of the main tasks is to estimate the accuracy of your result, and one common pattern is to iterate a routine that improves an estimate until it is sufficiently accurate. (You might ask yourself how Matlab knows how to do something as simple as estimate a value of sin(x) for some x.)

As a simple example, here is a short script that computes a radon transform of an image in Matlab. The radon transform takes projections of an image over a set of projection angles. Instead of trying to compute the projection along an arbitrary angle, I instead rotate the entire image using imrotate, so that the projection take is always vertical. Then we can take the projection simply using sum, since the sum of a matrix returns a vector containing the sum over each column.

You could write your own imrotate if you prefer, using interp2.

%# load a density map (image).
A = phantom;

n_pixels = size(A, 1);  %# image width (assume square)

%# At what rotation angles do we want to take projections?
n_thetas = 101;
thetas = linspace(0, 180, n_thetas);

result = zeros(n_thetas, n_pixels);

%# Loop over angles
for ii=1:length(thetas)
theta = thetas(ii);
rotated_image = imrotate(A, theta, 'crop');
result(ii, :) = sum(rotated_image);
end

%# display the result
imagesc(thetas, 1:n_pixels, result.');
xlabel('projection angle [degrees]');

What was once an integral of the density along a ray is now a sum over a column of a discretely sampled image, which was in turn found by interpolating the original image over a transformed coordinate system.

• Wow @nibot, thank you for such detailed answer. I have taken an "algorithms and data structures" class, but I don't see the relation between the material presented there and implementing scientific algorithms. I'll read the links you gave me and start practicing with simpler algorithms (from books instead of papers). Thanks again Aug 24 '11 at 14:14
• Hi Damian, I edited my answer to address your comment. I think you will find what you seek in a course or book on numerical methods / numerical analysis. Aug 24 '11 at 16:36
• Throughout answer! Aug 24 '11 at 17:24
• @nibot: thanks for the edit. I really like the numerical analysis course you linked. Why is the "finite impulse response filters" linked to interpolation? I wonder why this is not part of the curriculum as an EE student. Oh well. Thanks! Aug 24 '11 at 18:29
• @Damian: Sampling theory, interpolation/decimation, Z transforms, the bilinear transform, and FIR/IIR filters are taught in undergraduate EE classes/labs such as signals and systems, communication systems, linear control systems, and intro to DSP. I took numerical methods as part of dual degree program in computer engineering; I don't think it should be required of EEs in general. Aug 25 '11 at 3:35

Adding to nibot's excellent explanation, just a couple of more points.

• Numerical computing environments such as MATLAB, Octave or SciPy/NumPy will save you a lot of effort compared to doing it all by yourself in a generic programming language such as C++. Juggling with double arrays and loops just doesn't compare with having data types such as complex numbers and operations such as integrals at your fingertips. (It's doable for sure, and good C++ code can be an order of magnitude faster, with good library abstractions and templates it can even be reasonably clean and clear, but it's definitely easier to begin with e.g. MATLAB.)

• MATLAB also has "toolkits" for e.g. Image Processing and Digital Signal Processing, which may help a lot, depending on what you do.

• Mitra's Digital Signal Processing is a good book to learn (in MATLAB!) the basics of discrete time, filters, transforms, etc. which is pretty much mandatory knowledge to do any decent technical algorithms.
• Yes, I've read the documentation of the Image Processing Toolboox. I seems extremely useful, but my question was geared towards implementing something like that. Basically, I wanted to know how to take a mathematical algorithm/formula and implement it (like Mathworks did with the IPT). I wanted to know about the thought pattern or some guidelines. I will look at Mitra's book. Thanks! Aug 24 '11 at 14:19
• To add to the above answer, C++ toolkits such Armadillo can greatly simplify the conversion of Matlab code into fast C++ code. Armadillo's syntax is similar to Matlab. You can also mix'n'match Matlab and C++ code via Armadillo's mex interface. Jan 15 '14 at 3:33

Numerical methods. It's usually a upper-division university course and textbook.

DSP is usually near the intersection of numerical methods and efficient implementation. If you ignore efficiency, then what you might be looking for is any numerical approximation method that might produce an "accurate enough" result for the technical paper's equation(s) of interest. Sometimes one might be dealing with sampled data, where the sampling theorems will put some bounds on both the data acquisition method (pre-filtering) and the range or quality of results you can get given that data.

Sometimes Matlab, numerical recipes, or various image/signal processing libraries will have efficient algorithms or code for the desired numerical solution. But sometimes you may have to roll your own, so it helps to know the math behind various the numerical solution methods. And that's a big subject in its own right.