Slices, Stacks and ITK

I recently had to program some registration and decided to use ITK. However, there are a number of things I could not find well documented, so I thought I'd put stuff I had to figure out in this post, in the hope it will be useful to others.

Firstly to load a TIFF series or image stack, which is well documented in the ITK Software Guide (Section 7.11), the code looks like this:
 
typedef unsigned short PixelType;
typedef itk::Image< PixelType, 3 > ImageStackType;
typedef itk::ImageSeriesReader< ImageStackType > ReaderType;
typedef itk::NumericSeriesFileNames NameGeneratorType;

///Generate Numerical File Names
NameGeneratorType::Pointer nameGenerator = NameGeneratorType::New();
    nameGenerator->SetSeriesFormat( argv[1] );
    nameGenerator->SetStartIndex( first );
    nameGenerator->SetEndIndex( last );
    nameGenerator->SetIncrementIndex( 1 );

///Read Stack as TIFFs
ReaderType::Pointer stackReader = ReaderType::New();
    stackReader->SetImageIO( itk::TIFFImageIO::New() );
    stackReader->SetFileNames( nameGenerator->GetFileNames() );
    stackReader->Update();

Here argv[1] will be an sprintf format string like 'file%03d.tif' for reading files named file001.tif, file002.tif etc. The typedefs may appear to be overly used, but results in much leaner and reusable code.

To view the stack, use the ImageToVTKFilter (found following these intructions) like:

typedef itk::ImageToVTKImageFilter< ImageStackType >  StackConnectorType;

///Export to VTK
StackConnectorType::Pointer stackConnector = StackConnectorType::New();
    stackConnector->SetInput( imageStackFiltered );
    stackConnector->Update();
    stackConnector->GetOutput()->GetExtent(bounds);
    cout << "Size: " << bounds[1] << "x" << bounds[3] << endl;

///Display
DGVImageVTK *imageStackView = new DGVImageVTK;
    imageStackView->alignImages(false); //Has no effect?
    imageStackView->setName("Image Stack");
    imageStackView->SetInput(stackConnector->GetOutput());
    imageStackView->generateImage();
    imageStackView->show();

The result of the filter is passed to my DGV Image Class which wraps VTK. You can find DGV here. You should be able to pass the connector to vtkImageViewer2 class also.

There are two main posts/email-list entries that are useful for slice traversal and processing. The first is the method using itk::PasteImageFilter. I found this to work but contained large amounts of code and was very slow.

The second method was the JoinSeriesImageFilter, which a lot better but missed some newbie ITK stuff, which I missed and resulted in all slices of the result containing the same slice. The end working result is code that looks like:

///Iterators
typedef itk::ImageSliceConstIteratorWithIndex< ImageStackType > SliceConstIteratorType;

///Extractors

typedef itk::ExtractImageFilter< ImageStackType, ImageType > ExtractFilterType;
typedef itk::JoinSeriesImageFilter< ImageType, ImageStackType > JoinSeriesFilterType;

///Traverse through slices and produce new output stack
///Setup Slice Iterators which will iterate through slices in the stack
SliceConstIteratorType inIterator( imageStackFiltered, imageStackFiltered->GetLargestPossibleRegion() );
    inIterator.SetFirstDirection( 0 ); ///x axis
    inIterator.SetSecondDirection( 1 ); ///y axis

///Setup Image Stack that will be Joined together
JoinSeriesFilterType::Pointer joinSeries = JoinSeriesFilterType::New();
    joinSeries->SetOrigin( imageStackFiltered->GetOrigin()[2] );
    joinSeries->SetSpacing( imageStackFiltered->GetSpacing()[2] );

for(inIterator.GoToBegin(); !inIterator.IsAtEnd(); inIterator.NextSlice())
{
    //cout << inIterator.GetIndex() << endl;

    ///Setup region of the slice to extract
    ImageStackType::IndexType sliceIndex = inIterator.GetIndex();
    ExtractFilterType::InputImageRegionType::SizeType sliceSize = inIterator.GetRegion().GetSize();
        sliceSize[2] = 0;
    ExtractFilterType::InputImageRegionType sliceRegion = inIterator.GetRegion();
        sliceRegion.SetSize( sliceSize );
        sliceRegion.SetIndex( sliceIndex );

    ///Pull out slice
    ExtractFilterType::Pointer inExtractor = ExtractFilterType::New(); ///Must be within loop so that smart pointer is unique
        inExtractor->SetInput( imageStackFiltered );
        inExtractor->SetExtractionRegion( sliceRegion );
        inExtractor->Update();

    ///Operate on Slice
    InvertorType::Pointer invertor2 = InvertorType::New(); ///Must be within loop so that smart pointer is unique
        invertor2->SetInput( inExtractor->GetOutput() );
        invertor2->Update();

    ///Save Slice
    joinSeries->PushBackInput( invertor2->GetOutput() );
}

///----------
///Write out multi-page TIFF of the result
joinSeries->Update();
WriterType::Pointer writer = WriterType::New();
    writer->SetFileName( "registered_stack.tif" );
    writer->SetInput( joinSeries->GetOutput() );

try
{
    writer->Update();
}
catch( itk::ExceptionObject & err )
{
    std::cerr << "Write Output Exception caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
}

The iterator goes through all slices without needing to know anything about dimensions or sizes. You will have to tell the ExtractFilter the slice to extract, this is described using ImageRegions.

The operation applied to the slices is simply an inversion of the greyscales, especially useful if the slices are negative images. Note that the pointers are declared within the loop, this is important because the pointers remain valid only within the loop and automatically deleted since we are using SmartPointer's. The effect is that ExtractFilter then points correctly to the current slice. The result is written as a multi-page TIFF file, which you can open in ImageJ etc.

The above operation was used as a test, which will be replaced by registration code. If you are applying filters slice by slice, look into the SliceBySliceImageFilter that can be found in the Code/Review branch of ITK atm. To invert greyscales and then rescale intensities on each slice, you get:

typedef itk::InvertIntensityImageFilter< ImageType >  InvertorType;
typedef itk::RescaleIntensityImageFilter< ImageType >  RescaleIntensityType;

typedef itk::SliceBySliceImageFilter< ImageStackType, ImageStackType, InvertorType, RescaleIntensityType > SliceFilterType;

///Filter Stack
InvertorType::Pointer invertor = InvertorType::New(); ///Invert Greyscales
RescaleIntensityType::Pointer rescaler = RescaleIntensityType::New(); ///Normalise image values
    rescaler->SetOutputMinimum( 0 );
    rescaler->SetOutputMaximum( ImageMax );
    rescaler->SetInput( invertor->GetOutput() );

///Apply the filters to each slice of the stack
SliceFilterType::Pointer sliceFilter = SliceFilterType::New();
    sliceFilter->SetInput( 0, imageStack );
    sliceFilter->SetInput( 1, imageStack );
    sliceFilter->SetInputFilter( invertor );
    sliceFilter->SetOutputFilter( rescaler );
    sliceFilter->Update();

Thats it for the moment. More on registration later and code release later.

Hope that helps.
Cheers Shakes - L3mming

Packaging Qt Applications for Ubuntu/Debian

I have recently attempted to get my Discrete Geometry 3D Viewer capable of building Debian packages, so that one may install DGV without needing to compile it or worry about dependencies for other Ubuntu Distros (like the question I got from a potential user).

I believed there would be a substantial discussion on how to do this for Qt applications, but I could only find the Maemo Guide useful. It gets even more difficult if you want to do multiple binaries from a single source. Hence, I have documented my findings of this topic in this blog.

Useful Links that I used:
Complete Ubuntu Packaging Guide
Qt App Maemo Guide
GPG Guide
Pbuilder Howto
Using Local Packages

Firstly, my situation is the following. I have three side-by-side dynamic libraries, dgv-base, dgv-contrib and dgv-vtk. The libraries are divided based on dependencies of Qt, None, and dgv-base & Qt & VTK respectively by design. Then there's the actual DGV application which depend on these libraries.

First the multiple libraries. Begin by renaming the source directory, with the name of the package and the version (with a dash as the separator, this is important). For example

libdgv-0.15

Assuming that your source is in the state you wish to distribute it, create the tarballs so that you have two tarballs

libdgv-0.15.tar.gz
libdgv_0.15.orig.tar.gz

Note the underscore between the library name and the version. This is very important.

Change to the directory of the source and execute

dh_make -e your.maintainer@address -c GPL

where the "-c GPL" assumes you are using GPLv3 license and "your.maintainer@address" is your email address. This will ask you a series of questions, where you should select library if you're doing a library or a single binary if you're doing a binary. I will assume a library for the aforementioned reasons.

This step will create a directory called "debian" with all the Debian package configuration files. Remove the example files, they are not needed for what we are doing as far as I know.

rm *.ex *.EX

Please read these pages of the packaging guide for the remaining files. You need to edit the control and rules files as I have for DGV (also see Maemo guide for a simple single binary example). Things to watch out for is also given in this link (at the end). Fill the changelog and copyright files as described by the package guide, making sure to use the name of the overall package where applicable and to match the author email to your GPG name. If you do not have one, you need to create one using the GPG guide to sign your packages.

For multiple binaries, there is one last step. You need to create the ".install" and ".dir" file for each package. The former has the list of files to be installed, but in a wild card format (e.g. "usr/lib/lib*.so"). The latter is where the files are to be installed (e.g. "usr/lib").Note that there is no "/" in front as per normal Linux directories from root. This is because they will be relative paths and when the real package is built, it will build it from the root. See the debian files I created for DGV here.

Now to build your package, we will use Personal Builder (pbuilder), a way to build your package using only the minimal (initial) Ubuntu base/setup plus your dependencies. This is the best part that makes constructing packages most useful. If your dependencies are correct, then the user just installs the package and Synaptic or apt-get just installs everything you need to get it going. Once a pbuilder environment is setup (which can be made to suit for building on different architectures and Ubuntu distros), the package is built with your configuration and placed into

/var/cache/pbuilder/result

Only downside is that all the base packages will be downloaded from the Ubuntu repository, which could take a while. If you can't use the Internet for whatever reason, you can still build it the packages using

debuild

when within the source directory. You might also want to do this initially to ensure that its all working correctly.

First create the pbuilder environment by

sudo pbuilder create --distribution $(lsb_release -cs) \
        --othermirror "deb http://archive.ubuntu.com/ubuntu $(lsb_release -cs) main restricted universe multiverse"

Then build the descriptor file

debuild -S

If you did a "debuild" already, then this step is unnecessary. Once the descriptor file is present, then build using pbuilder

sudo pbuilder build *.dsc

To build as a different distro or architecture, see this. Hopefully, all works fine and you get a few packages. There might be a few warning from Lintian. Ensure that the descriptions lines are not too long or google the lintian warning. If your warning says "empty-binary-package" then your files are installed in the wrong places. For Qt projections, I install the files of the library into "debian/tmp", then move it out into the relevant package directories by using "dh_movefiles -p$@ -Xcontrib -Xvtk usr/lib" command in the rules (see the wiki). The "-Xitem" tell the movefiles app to ignore the files with the string "item" from the file names.

Finally, on a Live CD or clean install, check you packages by installing them. Pbuilder didn't pickup the fact that I had incorrectly named libvtk5.2 as libvtk5 for one of my packages, so I recommend it.

I then build the dgv application and its all done. Hope this helps. I will post up more on how to use local packages later when I have it working.

Cheers Shakes - L3mming

Doctrate Passed and New Paper

I have passed my doctorate! And without any changes! A lot faster than I thought I would too. Yay! :D

I have also put up my latest paper entitled "Fast Digital Convolutions using Bit-shifts". It extends the Number Theoretic Transforms (NTTs) known as Rader Transforms. You can find it on ariXiv.org here. Hopefully it is well accepted into the scientific community.

These new results will eventually end up into my FTL/NTTW library.

Cheers Shakes - L3mming

Ph.D Submitted

I have just submitted my thesis and it feels great. Its taken 3 and half years but its done! :D Still got the examination... but that's out of my hand now.

Thanks to my supervisors and my colleagues.

Expect to see releases for DGV and FTL libraries soon.

Cheers
Shakes - L3mming

Latex/Xelatex Watermarking

I was interested in getting a watermark to display on my thesis title page, but the documentation for the packages and examples were lacking.

Here's how to do it. Assuming your image is called 'ribs.png':

\usepackage{watermark}

in the preamble, then the following where you want the watermark.

\thiswatermark{%Put watermark on this page only
%\centering
\put(0,-660){\includegraphics[width=\textwidth]{ribs}}
}

You use the coordinate arguments of '\put' to place it in the right position. For transparency, I just editted the image in The Gimp image editor. Simply go to Tools->Colour Tools->Levels and adjust the main bar.

My title page looks like:

HTH
Cheers
Shakes - L3mming