Finally some segmentation. Flattening 3D volumes to 1D array

All this time I was thinking my region-growing procedure was running out of stack space. However, after some debugging I discovered that I was doing something wrong in the way I was transforming my 3D volume to a 1D array. Just sometime ago, when I was doing this course on Advanced graphics visualizations, I realized that you could store a 3D volume in a 1D array, and essentially the same goes with a 2D array. So if your 3D volume has dimensions x, y and z, a 1D array can be used to store the scalar value at any position (x,y,z). As you would have already imagined, the 1D array runs from index 0 to (x*y*z - 1). All you need is a mapping function that maps co-ordinates (x,y,z) to an index i in your 1D array. In your application, you might also require the inverse of this, i.e. from i to (x,y,z). The equations are quite straightforward:

Mapping function from (x,y,z) to array index i:

index i = (z * X * Y) + x + (X * y);

(where X, Y and Z are: 0<= x <=X and 0 <= y <= Y and 0 <= z <= Z)

and mapping from index i to (x,y,z):

z = i/(X * Y)

y = (i - z * X * Y)/(X)

x = i - (z * X * Y) - (y * X);

I was going wrong with these equations. Finally when I had them corrected, to ensure that these were correct, I wrote a verification program to test that the equations were invertible ( i.e. fg(x) = x if f = g^-1 ).

Both the recursive and non-recursive region-growing procedures now work fine. I had them tested on an artificial image of a brain MRI. I will be posting more results on the segmentation. Here are some preliminary results.

The procedure accepts both local and global thresholds. Local thresholds are thresholds imposed on the direct neighbors of a voxel, and these threholds vary from one voxel neighborhood to another. For e.g. I am using a threshold whereby every voxel in the region is allowed to include into the region a voxel in its neighborhood which is within a certain intensity +/- 5 for instance.

It so appears now, that if we try to pinpoint a small vessel and try to segment it using a small global threshold value, the entire volume gets segmented. However, the entire volume segmented comes out much "darker" than the original volume in the Maximum Intensity Projection (MIP) rendering.

The image on the left is a MIP rendering of the segmented volume by using a seed point close to the atrium of the heart, and by using threshold values of 60 100. The image on the right is a MIP rendering of the segmented volume by using the same seed point, but using thresholds of 0 and 1000. The image on the right is essentially the entire volume and technically with nothing segmented.

Installing ITK at home

I have been trying to get ITK (Imperial college's ITK NOT Insight Toolkit) running at my home using Visual .NET 2003. Apparently ITK has problems compiling in .NET 2005. The errors which show up are correctable and it becomes evident that .NET 2005 has improved their cl compiler. There are some errors that has to do with assigning a const char* to a char* to which .NET 2003 doesnt complain which it should have.

Getting down to important things to remember when installing ITK, one has to make sure that VTK and FLTK is compiled using the same compiler or else it might throw some errors. FLTK also needs to be compiled with the FLTK_USE_ZLIB, FLTK_USE_JPEG, FLTK_USE_PNG options set to off. This is assuming that our system does not have ZLIB, PNG, JPEG libraries and we make FLTK build the libraries for us.

It is also important to compile ITK under Release mode. There are different compile modes in Visual .NET 2003 (i.e. Debug, Release, etc). The IDE usually defaults itself to the Debug mode. Setting .NET 2003 compiler to the Release mode makes the compilation process a lot quicker. Technically we should be able to compile it in Debug mode, but I havent tested it.

If you are going to be using VTK volume rendering libraries, strangely the vtkVolumeRendering library (vtkVolumeRendering.lib file in the vtkVolumeRendering folder) is not loaded as an additional library by default to the visual .net 2003 linker. One can easily manually do this, by going to Project properties -> Linker -> Input.

One last word: watch out for the errors, ignore the warnings. Some of these errors can actually be fixed if looked into more closely.

More on Region growing for vessel segmentation

The problem of using a trivial region growing algorithm which uses a one or two-fold threshold as its growing-criteria is that it is bound to "explode" in the case of heart-vessels. What usually follows this explosion is a "stack-overflow" exception which results from the large amount of stack space taken up by the recursive region growing process. I have been experimenting with a recursive region-growing algorithm, and it runs out of stack when run on the visual .net C++ compiler despite increasing the stack size to several hundred megabytes (using the /STACK option in the linker, or the /F option of the compiler). What actually happens is that the region explodes into the atrium of the heart where it ultimately runs into a stack exception. However, the algorithm segments correctly without a stack-exception if the lower and upper-thresholds are set to point to high-intensities (scalar value of 80-150).

Grey-value intensities in contrast enhanced MR angiography is directly proportional to the amount of blood flow in a region. The atrium of the heart presents high intensities, the reason of which is quite obvious. However, the venous drainages also present high voxel intensities in certain regions. I presume this is related with the high blood-pressure in these areas due to the narrowing of the vessels towards the drainages. Although we would expect the entire drainage to be of high-intensity, however, this is not the case. Close to the endings, we come across low intensities. So now this variation in intensity poses a challenging problem when using thresholds, i.e if we were to select a threshold for region-growing segmentation of the venous endings, we would have no choice but to pick a low-threshold value (in-order to capture the minute vein endings) and a high-threshold in order to capture the parts of the vein-endings where the blood flow is high. Since the Atrium of the heart has a high-intensity value, this causes the algorithm to explode into the atrium, causing a stack-exception at some point. The stack-exception can also be equally attributed to the lower-threshold value, since this threshold value causes the algorithm to further leave the atrium and flow into the numerous vessels which branch out of the atrium.

Following is an illustration of a veinous drainage as seen through an MRI contrast-enhanced angiography. This is only a slice through the entire volume captured in the MRI.

After meeting with Prof. Daniel today, he suggested to me to further investigate the possibility of increasing the stack space and to see if the recursive algorithm can be made to work. The best possible result is to get the entire MRI segmented by using low-high threshold values. I also suggested to him how people have been performing brain vessels segmentation using an atlas-based region growing approach as described in "Region growing segmenation of brain vessels: an atlas based approach" by N. Passat et. al. Brain-vessel segmentation seems to be different from heart-vessel segmentation; in essence brain vessels dont branch out of an "atrium" like structure.

Prof. Daniel suggested to me to use a non-recursive approach to region-growing whereby we could keep the growing-region to a spherical-like shape at any time instant. A recursive process tends to grow out the region in one direction at a time. Making the region grow out spherically could allow us to develop better cost functions which for example could utilise the compactness (surface area / volume) of the growing-region.

The 9-month transfer report is due soon.

Experimenting with Region Growing methods

I have slacked off a bit due to tasmia's arrival to London and my exams. I havent been working for the past 4 days (which incl. 2 days of the weekend). I was able to install Visual studio 2005 at home, so which means I will be able to work on my experiements from home. My Google interview is soon, and so I will have to prepare for that as well.

I will be reading Region growing based techniques for volume visualization by R. Huang et. al. They have used simple statistics for region-growing criteria selection. I might get better results than the ones I have below (using simple thresholding).

I have used a recursive algorithm which uses a two-part threshold. The recursion is on the 6 neighbors of a voxel (2 in each direction -> 2 in x, 2 in y and 2 in z). The image above is a MIP rendering of the volume. The segmented region is show in dark one-color green.

Final exams and preparation week

I will soon be posting on the 3D segmentation I was able to do, only to a certain extent, using a recursive region growing algorithm. The criterion for region growing is simple and I am using an upper and lower threshold. This is known to produce "holes" in the MIP rendering of the segmented region. It seemed quite clear later on as to why such local discontinuity occurs in the rendering. Since we are doing a Maximum Intensity Projection, the maximum intensity function around a neighborhood of the segmented region is not guaranteed to be continuous or approximately continuous.

Using a recursive region growing algorithm with a criterion like the one I am using at the moment (thresholds) is just the beginning. It is popular to use other criterions, such as the Fischer's criterion. Most of these really boil down to using some sort of statistical mesaures. I am looking to first try and see how well it works if I used a simple idea of a pixel fulfilling a group if it is within some n number of standard deviations.

Apart from research work, I have been busy preparing for finals for the Machine vision and Advanced graphics and Visualization courses I took this year. Imperial exams are very thought provoking, much like the ones I "reckon" I used to get at U. Toronto.