Slicer3:Python:pitky

This is part of a NA-MIC Summer Project Week 2009 Project.

Example of an Image Writer

""" A simple example to show how to use weave with ITK.  This lets one
pass numpy arrays to ITK for processing.

based on vtk_example.py from weave distribution.

Author: Steve Pieper
Copyright (c) 2009, Steve Pieper
License: BSD Style.
"""

d = '/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/'
packages = ['matplotlib', 'neuroimaging', 'numpy', '', 'scipy']
for p in packages:
  sys.path.append(d+p)

import scipy.weave as weave
import numpy

import sys
import time
import os

# from ITKConfig.cmake - TODO: extract automatically, along with build flags etc (or invoke cmake via distutils?!)
inc_dirs = ['/Users/pieper/slicer3/latest/Slicer3-lib/Insight-build', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Code/Algorithms', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Code/BasicFilters', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Code/Common', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Code/Numerics', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Code/IO', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Code/Numerics/FEM', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Code/Numerics/Statistics', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Code/Numerics/NeuralNetworks', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Code/SpatialObject', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/MetaIO', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/NrrdIO', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight-build/Utilities/NrrdIO', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/DICOMParser', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight-build/Utilities/DICOMParser', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight-build/Utilities/expat', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/expat', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/nifti/niftilib', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/nifti/znzlib', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/itkExtHdrs', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight-build/Utilities', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/vxl/v3p/netlib', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/vxl/vcl', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/vxl/core', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight-build/Utilities/vxl/v3p/netlib', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight-build/Utilities/vxl/vcl', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight-build/Utilities/vxl/core', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight-build/Utilities/gdcm', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Utilities/gdcm/src', '/Users/pieper/slicer3/latest/Slicer3-lib/Insight/Code/Review']

# The ITK library directories.
lib_dirs = ['/Users/pieper/slicer3/latest/Slicer3-lib/Insight-build/bin']


def writer_test():
    """A simple example of how you can use ITK code in weave
       to write a volume to a file.
    """

    fileName = '/tmp/rand.nrrd'
    array = numpy.arange(0,24,dtype='float64').reshape(2,3,4)

    code=r"""
    typedef itk::Image<double,3> ImageType; // TODO: construct from dtype and shape
    ImageType::Pointer image = ImageType::New();
    image->GetPixelContainer()->SetImportPointer( array, Narray[0]*Narray[1]*Narray[2], false );

    ImageType::RegionType region;
    ImageType::IndexType index;
    ImageType::SizeType size;
    index[0] = index[1] = index[2] = 0;
    size[0] = Narray[0]; size[1] = Narray[1]; size[2] = Narray[2];
    region.SetIndex(index);
    region.SetSize(size);
    image->SetLargestPossibleRegion(region);
    image->SetBufferedRegion(region);


    itk::ImageFileWriter< ImageType >::Pointer writer = WriterType::New();
    writer->SetFileName( fileName );
    writer->SetInput( image );
    writer->Update();
    """

    weave.inline(code, ['fileName', 'array'],
                  include_dirs = inc_dirs,
                  library_dirs = lib_dirs,
                  headers=['"itkImage.h"', '"itkImageFileWriter.h"'],
                  libraries=['itkCommon', 'itkIO'])


TODO = """
- apply a filter (volume in and out, where output size may be different from input)
- compile on windows
- interface with cmake?
- how to deploy built libs in cpacked binary
-- how to trigger weave compilation during build process for packaging
"""


if __name__ == "__main__":
    writer_test()

Example of Nonlinear Registration

""" A simple example to show how to use weave with ITK.  This lets one
pass numpy arrays to ITK for processing.

based on pitky.py from steve Pieper

Author: Demian Wassermann
Copyright (c) 2010, Demian Wassermann
License: BSD Style.
"""

import scipy.weave as weave
import numpy as np

# from ITKConfig.cmake - TODO: extract automatically, along with build flags etc (or invoke cmake via distutils?!)
inc_dirs = ['/Users/demian/software/Slicer3-lib/Insight-build', '/Users/demian/software/Slicer3-lib/Insight/Code/Algorithms', '/Users/demian/software/Slicer3-lib/Insight/Code/BasicFilters', '/Users/demian/software/Slicer3-lib/Insight/Code/Common', '/Users/demian/software/Slicer3-lib/Insight/Code/Numerics', '/Users/demian/software/Slicer3-lib/Insight/Code/IO', '/Users/demian/software/Slicer3-lib/Insight/Code/Numerics/FEM', '/Users/demian/software/Slicer3-lib/Insight/Code/Numerics/Statistics', '/Users/demian/software/Slicer3-lib/Insight/Code/Numerics/NeuralNetworks', '/Users/demian/software/Slicer3-lib/Insight/Code/SpatialObject', '/Users/demian/software/Slicer3-lib/Insight/Utilities/MetaIO', '/Users/demian/software/Slicer3-lib/Insight/Utilities/NrrdIO', '/Users/demian/software/Slicer3-lib/Insight-build/Utilities/NrrdIO', '/Users/demian/software/Slicer3-lib/Insight/Utilities/DICOMParser', '/Users/demian/software/Slicer3-lib/Insight-build/Utilities/DICOMParser', '/Users/demian/software/Slicer3-lib/Insight-build/Utilities/expat', '/Users/demian/software/Slicer3-lib/Insight/Utilities/expat', '/Users/demian/software/Slicer3-lib/Insight/Utilities/nifti/niftilib', '/Users/demian/software/Slicer3-lib/Insight/Utilities/nifti/znzlib', '/Users/demian/software/Slicer3-lib/Insight/Utilities/itkExtHdrs', '/Users/demian/software/Slicer3-lib/Insight-build/Utilities', '/Users/demian/software/Slicer3-lib/Insight/Utilities', '/Users/demian/software/Slicer3-lib/Insight/Utilities/vxl/v3p/netlib', '/Users/demian/software/Slicer3-lib/Insight/Utilities/vxl/vcl', '/Users/demian/software/Slicer3-lib/Insight/Utilities/vxl/core', '/Users/demian/software/Slicer3-lib/Insight-build/Utilities/vxl/v3p/netlib', '/Users/demian/software/Slicer3-lib/Insight-build/Utilities/vxl/vcl', '/Users/demian/software/Slicer3-lib/Insight-build/Utilities/vxl/core', '/Users/demian/software/Slicer3-lib/Insight-build/Utilities/gdcm', '/Users/demian/software/Slicer3-lib/Insight/Utilities/gdcm/src', '/Users/demian/software/Slicer3-lib/Insight/Code/Review']

# The ITK library directories.
lib_dirs = ['/Users/demian/software/Slicer3-lib/Insight-build/bin']


def demons_registration( fixed, moving, output, deformation_field, sigmaDef = 1.5, sigmaUp = 0.0, numIterations = 10, maxStepLength = 2.):

    if any( im.dtype != np.dtype('float64') for im in (fixed,moving,output,deformation_field) ):
        raise ValueError('Images must be of dtype float64')

    if any( len(im.shape) != 3 for im in (fixed,moving,output,deformation_field[...,0])): 
        raise ValueError('Images must be 3-dimensional')

    if any( im.shape != fixed.shape for im in (fixed,moving,output,deformation_field[...,0])): 
        raise ValueError('Images must be all of the same size')

    if deformation_field.shape[-1] != 3:
        raise ValueError('deformation field must be a 3-component image')

    sigmaDef = float( sigmaDef)
    sigmaUp = float( sigmaUp )
    numIterations = int( numIterations )
    maxStepLength = float( maxStepLength )

    """A simple example of how you can use ITK code in weave
       to write a volume to a file.
    """

    code=r"""
    #line 45 "demons.py"

    #define Dimension                              (3)
    typedef double                                 PixelType;
    typedef itk::Image< PixelType, Dimension >     ImageType;

    typedef itk::Vector< PixelType, Dimension >    VectorPixelType;
    typedef itk::Image
     < VectorPixelType, Dimension >                DeformationFieldType;

    typedef itk::DiffeomorphicDemonsRegistrationFilter
       < ImageType, ImageType, DeformationFieldType >
                                                   ActualRegistrationFilterType;

    typedef ActualRegistrationFilterType::GradientType GradientType;


    ImageType::Pointer fixedImage = ImageType::New();
    fixedImage->GetPixelContainer()->SetImportPointer( fixed, Nfixed[0]*Nfixed[1]*Nfixed[2], false );


    ImageType::RegionType regionFixed;
    ImageType::IndexType indexFixed;
    ImageType::SizeType sizeFixed;
    indexFixed[0] = indexFixed[1] = indexFixed[2] = 0;

    sizeFixed[0] = Nfixed[0]; sizeFixed[1] = Nfixed[1]; sizeFixed[2] = Nfixed[2];
    regionFixed.SetIndex(indexFixed);
    regionFixed.SetSize(sizeFixed);
    fixedImage->SetLargestPossibleRegion(regionFixed);
    fixedImage->SetBufferedRegion(regionFixed);

    ImageType::Pointer movingImage = ImageType::New();
    movingImage->GetPixelContainer()->SetImportPointer( moving, Nmoving[0]*Nmoving[1]*Nmoving[2], false );

    ImageType::RegionType regionMoving;
    ImageType::IndexType indexMoving;
    ImageType::SizeType sizeMoving;
    indexMoving[0] = indexMoving[1] = indexMoving[2] = 0;

    sizeMoving[0] = Nmoving[0]; sizeMoving[1] = Nmoving[1]; sizeMoving[2] = Nmoving[2];
    regionMoving.SetIndex(indexMoving);
    regionMoving.SetSize(sizeMoving);
    movingImage->SetLargestPossibleRegion(regionMoving);
    movingImage->SetBufferedRegion(regionMoving);


    ImageType::Pointer outputImage = ImageType::New();
    outputImage->GetPixelContainer()->SetImportPointer( output, Noutput[0]*Noutput[1]*Noutput[2], false );

    ImageType::RegionType regionOutput;
    ImageType::IndexType indexOutput;
    ImageType::SizeType sizeOutput;
    indexOutput[0] = indexOutput[1] = indexOutput[2] = 0;

    sizeOutput[0] = Noutput[0]; sizeOutput[1] = Noutput[1]; sizeOutput[2] = Noutput[2];
    regionOutput.SetIndex(indexOutput);
    regionOutput.SetSize(sizeOutput);
    outputImage->SetLargestPossibleRegion(regionOutput);
    outputImage->SetBufferedRegion(regionOutput);

    
// As the GraftOutput didn't work for the deformation field, this paragraph can be ignored.
    DeformationFieldType::Pointer outputDeformationFieldImage = DeformationFieldType::New();
    outputDeformationFieldImage->GetPixelContainer()->SetImportPointer( 
        reinterpret_cast<VectorPixelType*>(deformation_field), 
        Ndeformation_field[0]*Ndeformation_field[1]*Ndeformation_field[2], 
    false );

    ImageType::RegionType regionOutputDeformation;
    ImageType::IndexType indexOutputDeformation;
    ImageType::SizeType sizeOutputDeformation;
    indexOutputDeformation[0] = indexOutputDeformation[1] = indexOutputDeformation[2] = 0;

    sizeOutputDeformation[0] = Ndeformation_field[0]; sizeOutputDeformation[1] = Ndeformation_field[1]; sizeOutputDeformation[2] = Ndeformation_field[2];
    regionOutputDeformation.SetIndex(indexOutputDeformation);
    regionOutputDeformation.SetSize(sizeOutputDeformation);
    outputDeformationFieldImage->SetLargestPossibleRegion(regionOutputDeformation);
    outputDeformationFieldImage->SetBufferedRegion(regionOutputDeformation);

    ActualRegistrationFilterType::Pointer filter
       = ActualRegistrationFilterType::New();

    filter->SetMaximumUpdateStepLength( maxStepLength );
    filter->SetUseGradientType( 
       static_cast<GradientType>(0) );

    filter->GraftOutput( outputDeformationFieldImage );
    filter->UseFirstOrderExpOn();

    if ( sigmaDef > 0.1 )
    {
        filter->SmoothDeformationFieldOn();
        filter->SetStandardDeviations( sigmaDef );
    }
    else
    {
        filter->SmoothDeformationFieldOff();
    }

    if ( sigmaUp > 0.1 )
    {
        filter->SmoothUpdateFieldOn();
        filter->SetUpdateFieldStandardDeviations( sigmaUp );
    }
    else
    {
        filter->SmoothUpdateFieldOff();
    }


    filter->SetFixedImage( fixedImage );
    filter->SetMovingImage( movingImage );
    filter->SetNumberOfIterations( numIterations );

    //This is made for memory efficiency and doesn't work :( 
    filter->GraftOutput( outputDeformationFieldImage );

    filter->UpdateLargestPossibleRegion();

    // warp the result
    typedef itk::WarpImageFilter
     < ImageType, ImageType, DeformationFieldType >  WarperType;
    WarperType::Pointer warper = WarperType::New();
    warper->SetInput( movingImage );
    warper->SetOutputSpacing( fixedImage->GetSpacing() );
    warper->SetOutputOrigin( fixedImage->GetOrigin() );
    warper->SetOutputDirection( fixedImage->GetDirection() );
    warper->SetDeformationField( filter->GetOutput() );
    warper->GraftOutput( outputImage );
    warper->UpdateLargestPossibleRegion();

   // copy resulting deformation field to output

   const size_t numPix = Ndeformation_field[0]*Ndeformation_field[1]*Ndeformation_field[2];
   const size_t bufferSize = sizeof( PixelType )  * 3 * numPix;

   memcpy( reinterpret_cast<void*>(deformation_field), reinterpret_cast<void*>(filter->GetOutput()->GetBufferPointer()), bufferSize);

    """

    weave.inline(code, ['fixed', 'moving','sigmaDef','sigmaUp','numIterations', 'maxStepLength', 'output', 'deformation_field' ],
                  include_dirs = inc_dirs,
                  library_dirs = lib_dirs,
                  headers=['"itkImage.h"', '"itkImageFileWriter.h"','"itkDiffeomorphicDemonsRegistrationFilter.h"'],
                  libraries=['itkCommon', 'itkIO'])



if __name__ == '__main__':
    import nifti

    dtype_ = np.dtype('float64') # 64-bit float or C++ double

    fixed = nifti.NiftiImage('fixed.nii.gz').getDataArray().astype(dtype_)
    moving = nifti.NiftiImage('moving.nii.gz').getDataArray().astype(dtype_)

    output = np.empty( fixed.shape, dtype = dtype_ )
    deformation_field = np.empty( fixed.shape + (3,), dtype = dtype_ )

    demons_registration( fixed, moving, output, deformation_field )