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 )