Proposals:VectorImage
ITK: Vector Image
What is a vector image
Create a templated n-dimensional vector image class.
This class differs from Image in that it is intended to represent multiple images. Each pixel represents 'k' measurements, each of datatype 'TPixel'. The memory organization of the resulting image is as follows:
... Pi0 Pi1 Pi2 Pi3 P(i+1)0 P(i+1)1 P(i+1)2 P(i+1)3 P(i+2)0 ...
where Pi0 represents the 0th measurement of the pixel at index i.
Conceptually, a VectorImage< double, 3 > is the same as a Image< Array< double >, 3 >. Creating the latter incurs a cost of allocating an Array for each pixel, which we would like to avoid.
Why ?
Excerpts from Gordon's email:
1. While a given DWI protocol has a fixed number of measurements, a user can reasonably want to analyze and compare DWI datasets with any number of different measurements without having to recompile.
2. The list of measurements should not add to the nominal dimension of the image, since, for example, the orientation information will not include that axis (e.g. DWI datasets have 3-D orientation information, not 4-D)
3. For reasons of memory efficiency, as well compatibility with memory-layout of other software, it would be ideal if the memory representations of the pixels' values were contiguous (so that the last DWI value for the Nth pixel proceeds the first DWI value for the N+1st pixel). It was the point that led Luis to reconsider the initial suggestion of simply creating an image of arrays (which also suffers a heavy one-time cost in allocating each pixel's array).
4. Using ITK to process the numerous existing NAMIC diffusion datasets in an efficient way. This is the most pressing case. Some people are already estimating tensors from DWI in ITK, but they're using volume-interleaved memory layout, which is opposite of what we would want. It would be very convenient if NAMIC participants could experiment with linear, non-linear, and other tensor estimation methods within a standard framework. Lacking such a framework, every NAMIC group is going to come up with their own non-interoperable methods.
5. Doing filtering or tractography in the DWI space, instead of in the space of estimated tensors. At least one published tractography method (Conturo PNAS '99) interpolate DWIs and estimate tensors at every step, and this has the advantage of dealing with noise in its own domain, instead of complicating it with the additional step of tensor estimation. Other methods (such as estimating two tensors per DWI list) also require access to the original DWIs.
Proposed Implementation.. feel free to add
The VectorImage class:
// The goal is to have a class that internally has the unfragmented memory organization
// described above, but exhibits the API of Image< Array< double >, N >. So it
// should return pixels of type Array< double > ( the container does not manage
// the memory in cases where are reference is returned, similar to returning a
// vnl_vector_ref). Same holds for other methods.
//
template <class TPixel, unsigned int VImageDimension=3 >
class ITK_EXPORT VectorImage : public Image< Array< TPixel >, VImageDimension >
{
...
// Override a few typedefs..
typedef Array< TPixel > PixelType;
typedef TPixel InternalPixelType;
typedef ImportImageContainer<
unsigned long, InternalPixelType> PixelContainer;
typedef VectorPixelAccessor<
PixelType > AccessorType;
...
// Override a few methods of image like Allocate().
// (Should this be a virtual method ?)
// This should reserve a buffer for ImageSize * k pixels, where k is the length
// of an Array. This will ensure that we have k contiguous blocks in memory
// rather an k pointers of itk::Array, which apart from fragmenting
// allocation space adds 8 bytes for storing the pointer.
//
// Override FillBuffer()
// Override Graft()
}
The VectorPixelAccessor class:
template class< TType >
class VectorPixelAccessor
{
...
typedef Array< TType > ExternalType; // and InternalType
...
// This will have Set/Get methods.
inline TType & Get( m_Buffer, m_Offset ) const
...
// Overloaded Additional methods
inline void Set( ouptut, input, m_Buffer, m_Offset )
{ // Set the value at m_Buffer + (k* m_Offset) }
...
private:
unsigned int k;
}
The ImageIterator class:
// minor change in API // The iterator does not know about that the pixels are contiguous blocks. All it knows is that there // is a X x Y x Z image with pixeltype of type VectorImage::InternalPixelType. We choose // to keep the same iterator class, and make the PixelAccessor the bearer of this information. // The method m_PixelAccessor.Get( *(m_Buffer + m_Offset )) will become // m_PixelAccessor.Get( m_Buffer, m_Offset ). // People using their custom adaptors will have broken API
Choices:
1. Above mentioned option
2. To Avoid breaking the API of existing pixel accessors, we can simply write a whole family of iterators to handle this image.
3. Remove VS6/7.0 and specialize the existing iterators for the VectorImage.
4. Please suggest
Was it worth the trouble ?
typedef itk::Image< itk::Array< double>, 3 > ArrayImageType1;
typedef itk::Image< itk::FixedArray< double, 3 >, 3 > ArrayImageType2;
ArrayImageType2::Pointer im2 = ArrayImageType2::New();
ArrayImageType1::Pointer im1 = ArrayImageType1::New();
//
{
ArrayImageType1::IndexType start;
itk::Array< double> f(3);
f.Fill( 3.0 );
start[0] = 0; // first index on X
start[1] = 0; // first index on Y
start[2] = 0; // first index on Z
ArrayImageType1::SizeType size;
size[0] = 200; // size along X
size[1] = 200; // size along Y
size[2] = 200; // size along Z
ArrayImageType1::RegionType region;
region.SetSize( size );
region.SetIndex( start );
im1->SetRegions( region );
im1->Allocate();
im1->FillBuffer( f );
}
//
{
ArrayImageType1::IndexType start;
itk::FixedArray< double, 3> f;
f.Fill(3.0);
start[0] = 0; // first index on X
start[1] = 0; // first index on Y
start[2] = 0; // first index on Z
ArrayImageType1::SizeType size;
size[0] = 200; // size along X
size[1] = 200; // size along Y
size[2] = 200; // size along Z
ArrayImageType1::RegionType region;
region.SetSize( size );
region.SetIndex( start );
im2->SetRegions( region );
im2->Allocate();
im2->FillBuffer( f );
}
2.01 s vs 0.59 s
Current status
1. Classes VectorImage, DefaultVectorPixelAccessor, VectorImageToImageAdaptor, VectorPixelAccessor added to ITK.
2. Gordon added the nrrdAxesPermute() function to NrrdIO, so within itkNrrdImageIO.cxx, it can take the DWI samples (currently along the slowest axis) and permute them up to the fastest axis. This assumes, of course, that there's an ITK image type that has a compatible memory layout
3. Added an extra level of indirection, the Functors (DefaultPixelAccessorFunctor and DefaultVectorPixelAccessorFunctor) to get iterators to work correctly with both images in combination with the DefaultVectorPixelAccessor and DefaultPixelAccessor.
