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