Ganesh Swami

As I’m only taking only one course next term, you might be wondering what I’m upto. I’ll be working at the Medical Imaging Analysis Lab (MIAL) for the next eight months. To be specific, I’ll be working in the field of Computational Anatomy, which is “the use of mathematical analysis to learn how tissues grow, assume new shapes and morph into mature structures.”

My supervisor, Dr. Faisal Beg’s PhD thesis was on an innovative method to compute large deformations metric mappings via the geodesic flows of diffeomorphisms. The algorithm is called Large Deformation Diffeomorphic Metric Mapping (LDDMM) and is in use in a number of places around the world.

Let me give an example: As we are growing up, doctors draw growth curves. These are curves of height and weight as a function of age. We can then compare this with a population average to see if everything is alright. Large deviations from the population statistics is an indication of a problem, like putting on too much weight (obesity), or growing too tall without putting on enough weight.

Dr. Beg’s long term vision is to do the same thing for other organs in the body. Height and weight are metrics in a vector space, and therefore are easy to compare. What about metrics about the size of the hippocampus, or the average orientation of ventricular myocardium fibers?

Why does this excite me so much?

1. At the center of LDDMM is a semi-Lagrangian numerical scheme borrowed from the weather simulation community. Semi-Lagrangian methods enable the use of large timesteps without affecting stability. Given my interest in turbulence models and fluid dynamics, this is going to be a treat.

2. We have our own Itanium 2 supercomputer. Because of the huge MRI datasets that we work on, a non-significant amount of time will be spent in communication. We can either come up with an elegant solution to reduce communication and make use of a beowulf-style cluster, or just throw the “right” hardware at the problem in the form of a shared memory multi-processor system. I think our computer has 100 processors and can be scaled upto 512 processors without much trouble.

3. We can leverage different high performance technologies at each level - SSE at the processor level, OpenMP between processors, and MPI above that. Doesn’t make much sense to use GCC on an Itanium, but we can definitely play with the new OpenMP support in GCC 4.2. I’ll once again get to visit the problem with computing inner products efficiently with SSE. Stay tuned.

4. Tricks. I’ve never had a chance to play with a VLIW architecture like the Itanium before. This opens up a wealth of opportunities to use tricks. What’s the fastest way to add two diffusion tensors?

5. Math. The metric “distance” is computed via the Euler-Lagrange PDE that arises from the variational problem setting. I have only done this in Euclidean and Polar manifolds. It should be interesting to see the E-L equation being solved for arbitrary diffeomorphism groups. Yes, the method has strong roots in group theory.

6. More math. The E-L equation arises from the minimization of the variational problem. This is an optimization problem. The naive cost function is not convex, so we are trying to do some work here. Maybe a manifold transformation? I don’t know.

7. Tons of applications to medical devices, surgery and analysis. One of the reasons I decided to take on this offer is because of the many connections and similarities we came up with when discussing some of my other path-planning projects with Dr. Beg.

Bring it on, I say!

I recently came across two very interesting projects relating to N-body astrophysics problems. Replace celestial bodies with atoms, and you have molecular dynamics.

The first is the GRAPE (GRavity PipE) project, which aims to build specialized hardware optimized for N-body code. From their preprint on arxiv, I can see three key innovations they’ve brought to the table:

• removal of local memory from the processing elements (PE)
• removal of inter PE communication
• reduce operations are done in hardware

What impresses me is that by introducing artificial constraints, they are forced to think a certain way. The removal of local memory lets them combine many PEs on the same chip, while the removal of inter PE communication forces them to be more economical with communication. Communication will eventually kill whatever peak performance you have with a single PE.

The second is the Maya Project for dense stellar systems. This is a book about writing a N-body simulator in Ruby. If you are familiar with N-body codes, then you’d know that a majority of the computational time is spent in the inter-particle interaction computations. This is typically a small fraction of the code base. They write these parts in C and do not loose too much performance in the process. There’s a wealth of information available here.

Basics

I’ve been trying out the imaging libraries available for Java for the past couple of days. I’m happy to say that I’ve liked what I’ve seen thus far. I’m sticking to this platform to develop the server side technology for my project.

The first task was to install the latest and greatest from Sun. This included Netbeans 5.5 and JDK 6 available as a single package. I also had to get the media libraries, which isn’t a part of the standard java kit for some reason. I had absolutely no problems installing all of this on a stock Debian Linux system. Now that Java is Free, in the future all of this should just be an apt-get away.

Netbeans is a kick-ass environment to develop with. I really like how unit testing is an integral part of the development process. Unlike Eclipse, Netbeans is not sluggish at all on a machine with 512 MB of RAM.

Experiment

I didn’t want to solve “toy problems” to test this environment. I decided to write a general image magnification function. Once again, instead of basing this on a simple super-sampling interpolation algorithm, I decided to try something more sophisticated.

In signal processing, a simple way to “zoom” into signals is by padding the frequency spectra of the signal with zeros. This method of zooming reduces the computational complexity from to due to the use of the Fast Fourier Transform.

Outline

We’ll assume that our domain of interest is always a power of two. Without loss of generality, we’ll also assume that the width and height are the same. Here’s how you’d use the FFT to do the zooming. I’ve also included the Matlab/Octave command equivalents in brackets to make it easier to understand:

1. Take the 2d-fft of the image. (fft2)
2. Make sure the DC components are in the center. If not shift-cycle it. (fftshift)
3. Add a border of zeros around the image. The final dimensions of the image is still a power of two.
4. Take the 2d-inverse fft of the image. (ifft2) Again, you might have to shift-cycle the image as in step 2.
5. This is the zoomed in image.

In programs like Matlab, the finite precision of the floating points leads to non-zero complex values in the inverse transform. These can be safely ignored.

Step three is the key step. What it does is add a bunch of zeros in the middle of the frequency spectrum. You could have added the pad right in the middle of the image, but as we’ll see in the next section, it’s simpler to add a border of zeros to the boundaries of the image.

Methodology

The Java Advanced Imaging (JAI) APIs define something known as an operator. This is the core of your signal processing. You have point operators, area operators, color quantization operators, edge extraction operators, statistical operators…you get the idea. All of these operators are constructed in a general way using a ParameterBlock. Here’s an example of using the AddConst operator for adding a constant “2.0″ to every pixel of an image:

ParameterBlock pb = new ParameterBlock ();
pb.addSource (src);
double [] val = {2.0};
pb.add (val);
PlanarImage dst = JAI.create ("AddConst", pb);

PlanarImage is the most general class for representing a 2-dimensional image. The sweet thing about JAI is that it builds a directed acyclic graph of operations, much like GEGL. This means that if you have an image reader, the image isn’t actually read until the pixel value is written to disk much later in the pipeline of operations. This is very powerful.

Using the outline, we build a general pipeline using the operators fileload, dft, PeriodicShift, border, PeriodicShift, idft, MultiplyConst, filesave. Done. We require the multiplication of the image by a constant because of the scaling factors in the forward and inverse transforms. Below, we see two images of sizes 256×256 and the zoomed in 512×512 pixel images. This means I zero-padded the frequency spectra by 128 pixels on each side.

The performance was not bad at all. During installation, I had to drop in a .so file, so I’m pretty sure the math intensive parts of JAI are optimized in a native language, and not Java.

This is pretty exciting because it emphasizes the UNIX approach to solving problems. Just like how you can combine sort, wc, grep and maybe other small tools to create a powerful end-product, various operators can be combined to build really cool image processing pipelines.

Suppose you want to acquire data in , then Nyquist’s sampling theorem places a lower bound on the number of samples required. Mathematically, with a basis set we can represent any function as an infinite sum:

where for an orthonormal basis, the coefficients are given as:

In reality, we represent the samples in an appropriate basis and ignore values below some threshold. So, if we are going to approximate signals anyway, why not make a smaller number of measurements rather than each point as required by the sampling theorem? This is exactly the point of Compressed Sensing.

One of the pioneers of Compressed Sampling, David Donobo gives an example:

For certain classes of images with m pixels, choosing an appropriate set of linear functions, with a reasonable margin of reconstruction error lets us get away with only n non-adaptive non-pixel samples for faithful recovery. The size of n is dramatically smaller than m.

The bound on n is .

A news article in SIAM News further explains how this reconstruction is accomplished by a norm minimization. It also talks about how cutting down the number of samples required has profound advantages — savings in money, time or even tissue damage (as in the case of X-rays.) It also has an example of a numerical experiment where they perfectly reconstruct a 512×512 image from just 512 Fourier coefficients, that is with over 95% of the data missing.

The compressed sensing website at Rice is a great resource.

Most search engines rank results based on their “authority.” Google for example, does this by the number of incoming high quality links. Considering the fact that Wikipedia these days is a source of high quality material, you find Wikipedia as the top result of many search terms.

Google’s results these days are less than optimal. Spammers have become quite sophisticated in their methods. Results are cluttered with iframe results, shell websites, and other useless listings.

Wikiasari is a new venture by the same guys behind Wikipedia. This is exactly what I’ve been looking forward to – an engine that prefers Wikipedia over other sources.

The Times quotes Jimmy Wales, the founder of Wikipedia:

Google is very good at many types of search, but in many instances it produces nothing but spam and useless crap. Try searching for the term ‘Tampa hotels’, for example, and you will not get any useful results.

The engine is built over Lucene, which is also behind Beagle, the Free indexer for Linux desktops. They also say the index will be released under the GFDL.