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

Discrete Geometry 3D Viewer

I thought I'd write list of features of my Discrete Geometry Viewer (DGV) program. I wrote this as a request from the New Projects section of the Linux Journal. The goal of DGV is to make visualisation of data as simple as possible (not something VTK is renowned for) while being intuitive and responsive.

In DGV, the image is placed into a 3D scene using OpenGL (via the VTK library), so it uses your graphics card to render the image. U can pan and zoom, adjust the gamma and rotate it. U can save what you're viewing or just export the data/image. U can also view the raw data in table form (may not be as useful to non-scientific users). From the table u can do surface and line plots. U can also add, subtract etc. images using the Operate feature. Finally, one can see the histogram of the image, which tells u statistical info on the image (output into console also).

The main feature of DGV is Filtering and Convolutions via Fast Fourier Transforms. This allows one to see the frequencies present within the image and edit them via the table. One can do crude Low-pass (smoothing), High-pass (edge detection) or band-pass filtering. DGV also features surface and line animations.

Future advancements will include:
- Pixel values within viewer
- Saving animations
- Python Shell rather than simple console output. See my project called QPythonShell which allows one to embed a Python shell into Qt applications.
- More file formats
- More transforms, like the Number Theoretic Transforms (via my new Number Theoretic Transform C library).

Links:
Homepage, Sourceforge, Google Code and Qt Apps.org.

Cheers
Shakes - L3mming