Slicer3:Python:pitky
From Slicer Wiki
Home < 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 )
