From cd1565deedb43f91203452c1f87fbb8839c112e4 Mon Sep 17 00:00:00 2001 From: Simon Rit Date: Wed, 31 Aug 2011 14:11:36 +0200 Subject: [PATCH] Removed clitkBSplineDeformableRegistration which was not supported anymore. Prefer using clitkBLUTDIR --itkbspline. --- .../clitkBSplineDeformableRegistration.cxx | 48 - .../clitkBSplineDeformableRegistration.ggo | 87 -- ...ineDeformableRegistrationGenericFilter.cxx | 51 - ...plineDeformableRegistrationGenericFilter.h | 120 --- ...ineDeformableRegistrationGenericFilter.txx | 978 ------------------ 5 files changed, 1284 deletions(-) delete mode 100644 registration/clitkBSplineDeformableRegistration.cxx delete mode 100644 registration/clitkBSplineDeformableRegistration.ggo delete mode 100644 registration/clitkBSplineDeformableRegistrationGenericFilter.cxx delete mode 100644 registration/clitkBSplineDeformableRegistrationGenericFilter.h delete mode 100644 registration/clitkBSplineDeformableRegistrationGenericFilter.txx diff --git a/registration/clitkBSplineDeformableRegistration.cxx b/registration/clitkBSplineDeformableRegistration.cxx deleted file mode 100644 index eb03bde..0000000 --- a/registration/clitkBSplineDeformableRegistration.cxx +++ /dev/null @@ -1,48 +0,0 @@ -/*========================================================================= - Program: vv http://www.creatis.insa-lyon.fr/rio/vv - - Authors belong to: - - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr - - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr - - This software is distributed WITHOUT ANY WARRANTY; without even - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the copyright notices for more information. - - It is distributed under dual licence - - - BSD See included LICENSE.txt file - - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================**/ -/** - * @file clitkBSplineDeformableRegistration.cxx - * @author Jef Vandemeulebroucke - * @date February 16 10:14:53 2009 - * - * @brief Optimized BSpline-FFD(itk=accellaration+additional memory requirements) or BLUT-FFD (clitk) - * - */ - -// clitk include -#include "clitkBSplineDeformableRegistration_ggo.h" -#include "clitkBSplineDeformableRegistrationGenericFilter.h" -#include "clitkIO.h" -#include "clitkImageCommon.h" - - -int main( int argc, char *argv[] ) -{ - //======================================================= - // Init command line - //======================================================= - GGO(clitkBSplineDeformableRegistration, args_info); - CLITK_INIT; - - //filter - clitk::BSplineDeformableRegistrationGenericFilter::Pointer genericFilter=clitk::BSplineDeformableRegistrationGenericFilter::New(); - genericFilter->SetArgsInfo(args_info); - genericFilter->Update(); - - return EXIT_SUCCESS; -} diff --git a/registration/clitkBSplineDeformableRegistration.ggo b/registration/clitkBSplineDeformableRegistration.ggo deleted file mode 100644 index 198f995..0000000 --- a/registration/clitkBSplineDeformableRegistration.ggo +++ /dev/null @@ -1,87 +0,0 @@ -#File clitkBSplineDeformableRegistration.ggo -#Author: Jef Vandemeulebroucke -#Date : Tue 15 July 2008 16.35 - -package "clitk" -version "Register 2 images using optimized ITK FFD (high memory requirements) or BLUT FFD. Both implementations are provided, mainly for debugging purposes.." - -option "config" - "Config file" string no - - -section "Run Time" - -option "verbose" v "Verbose" flag off -option "threads" - "Number of threads to use for intensive algorithms (default=min(cores,8))" int no - - -section "Input" - -option "reference" r "Input reference 3D image (float)" string yes -option "target" t "Input target 2D image (float)" string yes -option "mask" m "Mask to placed over the reference image" string no - - -section "Output" - -option "vf" - "Result DVF" string yes -option "coeff" - "Result coefficient images" string no multiple -option "output" o "Deformed target image" string yes -option "before" - "Difference image before (but after rigid transform)" string no -option "after" - "Difference image after " string no - - -section "Transform (Note that for BLUT FFD only one of --control, --spacing is required. The other will be adjusted to fit the region and allow exact representation. SamplingFactor will be set accordingly" - -option "init" - "Initial coefficient images" string no multiple -option "rigid" - "Prior rigid transform matrix from reference to target space" string no -option "order" - "Spline Order FFD" int no multiple default="3" -option "control" - "Internal control points for each dimension" int no multiple -option "spacing" - "Control point spacing for each dimension (mm)" double no multiple -option "wlut" - "With LUT for BSpline calculation" flag off -option "samplingFactor" - "LUT sampling factor" int no multiple - - -section "Interpolator" - -option "interp" - "Interpolation: 0=NN, 1=Linear, 2=BSpline, 3=BLUT" int no default="1" -option "interpOrder" - "Order if BLUT or BSpline (0-5)" int no default="3" -option "interpSF" - "Sampling factor if BLUT" int no default="20" - - -section "Metric (optimized, threaded (opt) versions are marked *)" - -option "metric" - "Type: 0=Mean-Squares*, 1=Normalized CC*, 2=Histogram CC, 3=Gradient-Difference, 4=Viola-Wells MI, 5=Histogram MI, 6=Mattes' MI*, 7=Normalized MI, 8=CR, 9=SSD for BLUT FFD**" int no default="9" -option "samples" - "Specify fraction [0, 1] of samples of the reference image used for the metric (* only). Use high fraction for detailed images (eg. 0.2, 0.5), for smooth images 0.01 might be enough." float no default="1" -option "intThreshold" - "Fixed image samples intensity threshold (* only)" float no -option "subtractMean" - "1: Subtract mean for NCC calculation (narrows optimal)" flag on -option "bins" - "2,5-8: Number of histogram bins" int no default="50" -option "random" - "4,6: Samples should be taken randomly, otherwise uniformly" flag off -option "stdDev" - "4: specify the standard deviation in mm of the gaussian kernels for both PDF estimations" float no default="0.4" -option "explicitPDFDerivatives" - "6: Calculate PDF derivatives explicitly (rigid=true; FFD=false)" flag off - - -section "Optimizer" - -option "optimizer" - "0=Simplex, 1=Powell, 2=FRPR, 3=Regular Step GD, 4=VersorRigid3D, 5=Conjugated Gradient, 6=L-BFGS, 7=L-BFGS-B" int no default="7" -option "delta" - "0: Initial delta, otherwise automatic" double no -option "step" - "1,2,3,4: Initial stepsize (to be multiplied with the gradient)" double no default="2.0" -option "relax" - "3,4: Relaxation of the stepsize (multiplied each time the gradient changes sign)" double no default="0.7" -option "valueTol" - "0,1,2: Tolerance on the function" double no default="0.01" -option "stepTol" - "0,1,3,4: Tolerance on the step size" double no default="0.1" -option "gradTol" - "3,4,6,7: Tolerance on the (projected) gradient magnitude (7: 1=low->1e-10=high precision)" double no default="1e-5" -option "lineAcc" - "6: Line accuracy (eg: high=0.1, low=0.9)" double no default="0.9" -option "convFactor" - "7: Convergence factor: terminate if factor*machine_precision>reduction in cost (1e+12 low -> 1e+1 high precision) " double no default="1e+7" -option "maxIt" - "0-7: Maximum number of iterations" int no default="500" -option "maxLineIt" - "1,2: Maximum number of line iterations" int no default="50" -option "maxEval" - "6,7: Maximum number of evaluations" int no default="500" -option "maxCorr" - "7: Maximum number of corrections" int no default="5" -option "selectBound" - "7: Select the type of bound: 0=none, 1=u, 2=u&l, 3=l" int no default="0" -option "lowerBound" - "7: The lower bound" double no default="0.0" -option "upperBound" - "7: The upper bound" double no default="0.0" - - -section "Registration" - -option "levels" - "Number of resolution levels" int no default="1" - - diff --git a/registration/clitkBSplineDeformableRegistrationGenericFilter.cxx b/registration/clitkBSplineDeformableRegistrationGenericFilter.cxx deleted file mode 100644 index c6ac94c..0000000 --- a/registration/clitkBSplineDeformableRegistrationGenericFilter.cxx +++ /dev/null @@ -1,51 +0,0 @@ -/*========================================================================= - Program: vv http://www.creatis.insa-lyon.fr/rio/vv - - Authors belong to: - - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr - - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr - - This software is distributed WITHOUT ANY WARRANTY; without even - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the copyright notices for more information. - - It is distributed under dual licence - - - BSD See included LICENSE.txt file - - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================**/ -#ifndef _clitkBSplineDeformableRegistrationGenericFilter_cxx -#define _clitkBSplineDeformableRegistrationGenericFilter_cxx -#include "clitkBSplineDeformableRegistrationGenericFilter.h" - - -namespace clitk { - - clitk::BSplineDeformableRegistrationGenericFilter::BSplineDeformableRegistrationGenericFilter() - { - m_Verbose=false; - } - - - void clitk::BSplineDeformableRegistrationGenericFilter::Update() - { - //Get the image Dimension and PixelType - int Dimension; - std::string PixelType; - - clitk::ReadImageDimensionAndPixelType(m_ReferenceFileName, Dimension, PixelType); - - if(Dimension==2) UpdateWithDim<2>(PixelType); - else if(Dimension==3) UpdateWithDim<3>(PixelType); - //else if (Dimension==4)UpdateWithDim<4>(PixelType); - else - { - std::cout<<"Error, Only for 2 and 3 Dimensions!!!"< - * @date 14 March 2009 - * - * @brief - * - =================================================*/ - -// clitk -#include "clitkBSplineDeformableRegistration_ggo.h" -#include "clitkIO.h" -#include "clitkImageCommon.h" -#include "clitkDifferenceImageFilter.h" -#include "clitkTransformUtilities.h" -#include "clitkLBFGSBOptimizer.h" -#include "clitkBSplineDeformableTransform.h" -#include "clitkGenericOptimizer.h" -#include "clitkGenericInterpolator.h" -#include "clitkGenericMetric.h" - -// itk include -#include "itkMultiResolutionImageRegistrationMethod.h" -#include "itkMultiResolutionPyramidImageFilter.h" -#include "itkImage.h" -#include "itkBSplineDeformableTransform.h" -#include "itkImageFileReader.h" -#include "itkImageFileWriter.h" -#include "itkResampleImageFilter.h" -#include "itkCommand.h" -#include "itkImageMaskSpatialObject.h" -#include "itkEuler3DTransform.h" -#include "itkWarpImageFilter.h" -#include "itkLightObject.h" -#include "itkImageToImageMetric.h" -#include "itkInterpolateImageFunction.h" -#include "itkLabelStatisticsImageFilter.h" - - -namespace clitk -{ - - class ITK_EXPORT BSplineDeformableRegistrationGenericFilter : public itk::LightObject - - { - public: - typedef BSplineDeformableRegistrationGenericFilter Self; - typedef itk::LightObject Superclass; - typedef itk::SmartPointer Pointer; - typedef itk::SmartPointer ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Run-time type information (and related methods) */ - itkTypeMacro( BSplineDeformableRegistrationGenericFilter, LightObject ); - - - //==================================================================== - // Set methods - void SetArgsInfo(const args_info_clitkBSplineDeformableRegistration& a) - { - m_ArgsInfo=a; - m_ReferenceFileName=m_ArgsInfo.reference_arg; - m_Verbose=m_ArgsInfo.verbose_flag; - } - - //==================================================================== - // Update - virtual void Update(); - - protected: - //const char * GetNameOfClass() const { return "BSplineDeformableRegistrationGenericFilter"; } - - //==================================================================== - // Constructor & Destructor - BSplineDeformableRegistrationGenericFilter(); - ~BSplineDeformableRegistrationGenericFilter(){;} - - //==================================================================== - //Protected member functions - template void UpdateWithDim(std::string PixelType); - template void UpdateWithDimAndPixelType(); - - args_info_clitkBSplineDeformableRegistration m_ArgsInfo; - bool m_Verbose; - std::string m_ReferenceFileName; - - }; - -} // end namespace clitk -#ifndef ITK_MANUAL_INSTANTIATION -#include "clitkBSplineDeformableRegistrationGenericFilter.txx" -#endif - -#endif //#define _clitkBSplineDeformableRegistrationGenericFilter_h - - - - diff --git a/registration/clitkBSplineDeformableRegistrationGenericFilter.txx b/registration/clitkBSplineDeformableRegistrationGenericFilter.txx deleted file mode 100644 index a0f0d7a..0000000 --- a/registration/clitkBSplineDeformableRegistrationGenericFilter.txx +++ /dev/null @@ -1,978 +0,0 @@ -/*========================================================================= - Program: vv http://www.creatis.insa-lyon.fr/rio/vv - - Authors belong to: - - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr - - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr - - This software is distributed WITHOUT ANY WARRANTY; without even - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the copyright notices for more information. - - It is distributed under dual licence - - - BSD See included LICENSE.txt file - - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================**/ -#ifndef __clitkBSplineDeformableRegistrationGenericFilter_txx -#define __clitkBSplineDeformableRegistrationGenericFilter_txx -#include "clitkBSplineDeformableRegistrationGenericFilter.h" - - -namespace clitk -{ - - //============================================================================== - //Creating an observer class that allows us to change parameters at subsequent levels - //============================================================================== - template - class RegistrationInterfaceCommand : public itk::Command - { - public: - typedef RegistrationInterfaceCommand Self; - - - typedef itk::Command Superclass; - typedef itk::SmartPointer Pointer; - itkNewMacro( Self ); - protected: - RegistrationInterfaceCommand() {}; - public: - typedef TRegistration RegistrationType; - typedef RegistrationType * RegistrationPointer; - - // Two arguments are passed to the Execute() method: the first - // is the pointer to the object which invoked the event and the - // second is the event that was invoked. - void Execute(itk::Object * object, const itk::EventObject & event) - { - if( !(itk::IterationEvent().CheckEvent( &event )) ) - { - return; - } - RegistrationPointer registration = dynamic_cast( object ); - unsigned int numberOfLevels=registration->GetNumberOfLevels(); - unsigned int currentLevel=registration->GetCurrentLevel()+1; - std::cout< Pointer; - itkNewMacro( Self ); - protected: - CommandIterationUpdate() {}; - public: - typedef clitk::GenericOptimizer OptimizerType; - typedef const OptimizerType * OptimizerPointer; - - // We set the generic optimizer - void SetOptimizer(OptimizerPointer o){m_Optimizer=o;} - - // Execute - void Execute(itk::Object *caller, const itk::EventObject & event) - { - Execute( (const itk::Object *)caller, event); - } - - void Execute(const itk::Object * object, const itk::EventObject & event) - { - if( !(itk::IterationEvent().CheckEvent( &event )) ) - { - return; - } - - m_Optimizer->OutputIterationInfo(); - } - - OptimizerPointer m_Optimizer; - }; - - - //============================================================================== - // Update with the number of dimensions - //============================================================================== - template - void BSplineDeformableRegistrationGenericFilter::UpdateWithDim(std::string PixelType) - { - - if (m_Verbose) std::cout << "Images were detected to be "<< Dimension << "D and " << PixelType << "..." << std::endl; - - if(PixelType == "short"){ - if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed short..." << std::endl; - UpdateWithDimAndPixelType(); - } - // else if(PixelType == "unsigned_short"){ - // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_short..." << std::endl; - // UpdateWithDimAndPixelType(); - // } - - // else if (PixelType == "unsigned_char"){ - // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_char..." << std::endl; - // UpdateWithDimAndPixelType(); - // } - - // else if (PixelType == "char"){ - // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed_char..." << std::endl; - // UpdateWithDimAndPixelType(); - // } - else { - if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl; - UpdateWithDimAndPixelType(); - } - } - - - - //============================================================================== - // Update with the number of dimensions and pixeltype - //============================================================================== - template - void BSplineDeformableRegistrationGenericFilter::UpdateWithDimAndPixelType() - { - - //======================================================= - // Run-time - //======================================================= - bool threadsGiven=m_ArgsInfo.threads_given; - int threads=m_ArgsInfo.threads_arg; - - typedef itk::Image< PixelType, ImageDimension > FixedImageType; - typedef itk::Image< PixelType, ImageDimension > MovingImageType; - const unsigned int SpaceDimension = ImageDimension; - typedef double TCoordRep; - - - //======================================================= - //Input - //======================================================= - typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType; - typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType; - - typename FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New(); - typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New(); - - fixedImageReader->SetFileName( m_ArgsInfo.reference_arg ); - movingImageReader->SetFileName( m_ArgsInfo.target_arg ); - if (m_Verbose) std::cout<<"Reading images..."<Update(); - movingImageReader->Update(); - - typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput(); - typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput(); - typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion(); - - // The metric region: where should the metric be CALCULATED (depends on mask) - typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion(); - typename FixedImageType::RegionType::SizeType metricRegionSize=metricRegion.GetSize(); - typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex(); - typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin(); - - // The transform region: where should the transform be DEFINED (depends on mask) - typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion(); - typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize(); - typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex(); - typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin(); - - - //======================================================= - // If given, we connect a mask to the fixed image - //====================================================== - typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType; - typename MaskType::Pointer spatialObjectMask=NULL; - - if (m_ArgsInfo.mask_given) - { - typedef itk::Image< unsigned char, ImageDimension > ImageMaskType; - typedef itk::ImageFileReader< ImageMaskType > MaskReaderType; - typename MaskReaderType::Pointer maskReader = MaskReaderType::New(); - maskReader->SetFileName(m_ArgsInfo.mask_arg); - - try - { - maskReader->Update(); - } - catch( itk::ExceptionObject & err ) - { - std::cerr << "ExceptionObject caught while reading mask !" << std::endl; - std::cerr << err << std::endl; - return; - } - if (m_Verbose)std::cout <<"Reference image mask was read..." <SetImage( maskReader->GetOutput() ); - - // Find the bounding box of the "inside" label - typedef itk::LabelStatisticsImageFilter StatisticsImageFilterType; - typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New(); - statisticsImageFilter->SetInput(maskReader->GetOutput()); - statisticsImageFilter->SetLabelInput(maskReader->GetOutput()); - statisticsImageFilter->Update(); - typename StatisticsImageFilterType::BoundingBoxType boundingBox = statisticsImageFilter->GetBoundingBox(1); - - // Limit the transform region to the mask - for (unsigned int i=0; iTransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin); - - // Limit the metric region to the mask - metricRegion=transformRegion; - fixedImage->TransformIndexToPhysicalPoint(metricRegion.GetIndex(), metricRegionOrigin); - - } - - - - //======================================================= - // Regions - //===================================================== - if (m_Verbose) - { - // Fixed image region - std::cout<<"The fixed image has its origin at "<GetOrigin()< FixedImagePyramidType; - typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType; - typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New(); - typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New(); - - - // //======================================================= - // // Rigid Transform - // //======================================================= - // typedef itk::Euler3DTransform RigidTransformType; - // RigidTransformType::Pointer rigidTransform=RigidTransformType::New(); - - // if (m_ArgsInfo.rigid_given) - // { - // itk::Matrix rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.rigid_arg); - - // //Set the rotation - // itk::Matrix finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix); - // rigidTransform->SetMatrix(finalRotation); - - // //Set the translation - // itk::Vector finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix); - // rigidTransform->SetTranslation(finalTranslation); - - // } - - - //======================================================= - // BSpline Transform - //======================================================= - typename FixedImageType::RegionType::SizeType splineOrders ; - - //Default is cubic splines - splineOrders.Fill(3); - if (m_ArgsInfo.order_given) - for(unsigned int i=0; i TransformType; - typename TransformType::Pointer transform; - typedef itk::BSplineDeformableTransform BSplineTransformType; - typedef BSplineTransformType* BSplineTransformPointer; - typedef clitk::BSplineDeformableTransform BLUTTransformType; - typedef BLUTTransformType* BLUTTransformPointer; - - // JV parameter array is passed by reference, create outside context so it exists afterwards!!!!! - typedef typename TransformType::ParametersType ParametersType; - ParametersType parameters; - - - // CLITK BLUT transform - if(m_ArgsInfo.wlut_flag) - { - typename BLUTTransformType::Pointer bsplineTransform = BLUTTransformType::New(); - if (m_Verbose) std::cout<<"Setting the spline orders to "<SetSplineOrders(splineOrders); - - //------------------------------------------------------------------------- - // Define the region: Either the spacing or the number of CP should be given - //------------------------------------------------------------------------- - - // Region - typedef typename BSplineTransformType::RegionType RegionType; - RegionType bsplineRegion; - typename RegionType::SizeType gridSizeOnImage; - typename RegionType::SizeType gridBorderSize; - typename RegionType::SizeType totalGridSize; - - // Spacing - typedef typename BSplineTransformType::SpacingType SpacingType; - SpacingType fixedImageSpacing, chosenSpacing, adaptedSpacing; - fixedImageSpacing = fixedImage->GetSpacing(); - - // Only spacing given: adjust if necessary - if (m_ArgsInfo.spacing_given && !m_ArgsInfo.control_given) - { - for(unsigned int r=0; r(chosenSpacing[r]/fixedImageSpacing[r]) ) ); - adaptedSpacing[r]= ( itk::Math::Round(chosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ; - } - if (m_Verbose) std::cout<<"The chosen control point spacing "<GetDirection(); - SpacingType gridOriginOffset = gridDirection * adaptedSpacing; - - // Origin: 1 CP border for spatial dimensions - typedef typename BSplineTransformType::OriginType OriginType; - OriginType gridOrigin = transformRegionOrigin - gridOriginOffset; - if (m_Verbose) std::cout<<"The control point grid origin was set to "<SetGridSpacing( adaptedSpacing ); - bsplineTransform->SetGridOrigin( gridOrigin ); - bsplineTransform->SetGridRegion( bsplineRegion ); - bsplineTransform->SetGridDirection( gridDirection ); - - //Bulk transform - //if (m_Verbose) std::cout<<"Setting rigid transform..."<SetBulkTransform( rigidTransform ); - - //Vector BSpline interpolator - //bsplineTransform->SetOutputSpacing(fixedImage->GetSpacing()); - typename RegionType::SizeType samplingFactors; - for (unsigned int i=0; i< ImageDimension; i++) - { - if (m_Verbose) std::cout<<"For dimension "<GetSpacing()[i]); - if (m_Verbose) std::cout<<"Setting sampling factor "<SetLUTSamplingFactors(samplingFactors); - - //initial parameters - if (m_ArgsInfo.init_given) - { - typedef itk::ImageFileReader CoefficientReaderType; - typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New(); - coeffReader->SetFileName(m_ArgsInfo.init_arg[0]); - coeffReader->Update(); - bsplineTransform->SetCoefficientImage(coeffReader->GetOutput()); - } - else - { - //typedef typename TransformType::ParametersType ParametersType; - const unsigned int numberOfParameters = bsplineTransform->GetNumberOfParameters(); - parameters=ParametersType( numberOfParameters ); - parameters.Fill( 0.0 ); - bsplineTransform->SetParameters( parameters ); - } - - // Mask - if (spatialObjectMask) bsplineTransform->SetMask( spatialObjectMask ); - - // Pass - transform=bsplineTransform; - } - - //ITK BSpline transform - else - { - typename BSplineTransformType::Pointer bsplineTransform = BSplineTransformType::New(); - - // Define the region - typedef typename BSplineTransformType::RegionType RegionType; - RegionType bsplineRegion; - typename RegionType::SizeType gridSizeOnImage; - typename RegionType::SizeType gridBorderSize; - typename RegionType::SizeType totalGridSize; - -#if 0 - //Set the number of control points - for(unsigned int r=0; rGetSpacing(); - typename FixedImageType::SizeType fixedImageSize = fixedImageRegion.GetSize(); - if (m_ArgsInfo.spacing_given) - { - - for(unsigned int r=0; r(fixedImageSize[r] - 1) / - static_cast(gridSizeOnImage[r] - 1); - } - } - if (m_Verbose) std::cout<<"The control points spacing was set to "<GetDirection(); - SpacingType gridOriginOffset = gridDirection * spacing; - - // Origin - typedef typename BSplineTransformType::OriginType OriginType; - OriginType origin = fixedImage->GetOrigin(); - OriginType gridOrigin = origin - gridOriginOffset; - - // Set - bsplineTransform->SetGridSpacing( spacing ); - bsplineTransform->SetGridOrigin( gridOrigin ); - bsplineTransform->SetGridRegion( bsplineRegion ); - bsplineTransform->SetGridDirection( gridDirection ); - - // Bulk transform - // if (m_Verbose) std::cout<<"Setting rigid transform..."<SetBulkTransform( rigidTransform ); - - // Initial parameters - if (m_ArgsInfo.init_given) - { - typedef itk::ImageFileReader CoefficientReaderType; - typename BSplineTransformType::ImageType::Pointer coeffImages[SpaceDimension]; - for(unsigned int i=0; iSetFileName(m_ArgsInfo.init_arg[i]); - coeffReader->Update(); - coeffImages[i]=coeffReader->GetOutput(); - } - bsplineTransform->SetCoefficientImage(coeffImages); - } - else - { - const unsigned int numberOfParameters = bsplineTransform->GetNumberOfParameters(); - parameters=ParametersType( numberOfParameters ); - parameters.Fill( 0.0 ); - bsplineTransform->SetParameters( parameters ); - } - - // Pass - transform=bsplineTransform; -#else - // Spacing - typedef typename BSplineTransformType::SpacingType SpacingType; - SpacingType fixedImageSpacing, chosenSpacing, adaptedSpacing; - fixedImageSpacing = fixedImage->GetSpacing(); - - // Only spacing given: adjust if necessary - if (m_ArgsInfo.spacing_given && !m_ArgsInfo.control_given) - { - for(unsigned int r=0; rGetDirection(); - SpacingType gridOriginOffset = gridDirection * adaptedSpacing; - - // Origin: 1 CP border for spatial dimensions - typedef typename BSplineTransformType::OriginType OriginType; - OriginType gridOrigin = transformRegionOrigin - gridOriginOffset; - if (m_Verbose) std::cout<<"The control point grid origin was set to "<SetGridSpacing( adaptedSpacing ); - bsplineTransform->SetGridOrigin( gridOrigin ); - bsplineTransform->SetGridRegion( bsplineRegion ); - bsplineTransform->SetGridDirection( gridDirection ); - - //Bulk transform - //if (m_Verbose) std::cout<<"Setting rigid transform..."<SetBulkTransform( rigidTransform ); - - //Vector BSpline interpolator - //bsplineTransform->SetOutputSpacing(fixedImage->GetSpacing()); - - // Initial parameters - if (m_ArgsInfo.init_given) - { - typedef itk::ImageFileReader CoefficientReaderType; - typename BSplineTransformType::ImageType::Pointer coeffImages[SpaceDimension]; - for(unsigned int i=0; iSetFileName(m_ArgsInfo.init_arg[i]); - coeffReader->Update(); - coeffImages[i]=coeffReader->GetOutput(); - } - bsplineTransform->SetCoefficientImage(coeffImages); - } - else - { - const unsigned int numberOfParameters = bsplineTransform->GetNumberOfParameters(); - parameters=ParametersType( numberOfParameters ); - parameters.Fill( 0.0 ); - bsplineTransform->SetParameters( parameters ); - } - - // Pass - transform=bsplineTransform; -#endif - - } - - - //======================================================= - // Interpolator - //======================================================= - typedef clitk::GenericInterpolator GenericInterpolatorType; - typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New(); - genericInterpolator->SetArgsInfo(m_ArgsInfo); - typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep > InterpolatorType; - typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer(); - - - //======================================================= - // Metric - //======================================================= - typedef clitk::GenericMetric GenericMetricType; - typename GenericMetricType::Pointer genericMetric=GenericMetricType::New(); - genericMetric->SetArgsInfo(m_ArgsInfo); - genericMetric->SetFixedImage(fixedImage); - if (spatialObjectMask) genericMetric->SetFixedImageMask( spatialObjectMask ); - genericMetric->SetFixedImageRegion(metricRegion); - typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType; - typename MetricType::Pointer metric=genericMetric->GetMetricPointer(); - -#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS - if (threadsGiven) metric->SetNumberOfThreads( threads ); -#else - if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."< GenericOptimizerType; - GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New(); - genericOptimizer->SetArgsInfo(m_ArgsInfo); - genericOptimizer->SetMaximize(genericMetric->GetMaximize()); - genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters()); - typedef itk::SingleValuedNonLinearOptimizer OptimizerType; - OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer(); - - - //======================================================= - // Registration - //======================================================= - typedef itk::MultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType; - typename RegistrationType::Pointer registration = RegistrationType::New(); - registration->SetMetric( metric ); - registration->SetOptimizer( optimizer ); - registration->SetInterpolator( interpolator ); - registration->SetTransform (transform); - if(threadsGiven) registration->SetNumberOfThreads(threads); - registration->SetFixedImage( fixedImage ); - registration->SetMovingImage( movingImage ); - registration->SetFixedImageRegion( metricRegion ); - registration->SetFixedImagePyramid( fixedImagePyramid ); - registration->SetMovingImagePyramid( movingImagePyramid ); - registration->SetInitialTransformParameters( transform->GetParameters() ); - registration->SetNumberOfLevels(m_ArgsInfo.levels_arg); - if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<SetOptimizer(genericOptimizer); - optimizer->AddObserver( itk::IterationEvent(), observer ); - - - // Output level info - typedef RegistrationInterfaceCommand CommandType; - typename CommandType::Pointer command = CommandType::New(); - registration->AddObserver( itk::IterationEvent(), command ); - } - - - //======================================================= - // Let's go - //======================================================= - if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl; - - try - { - registration->StartRegistration(); - } - catch( itk::ExceptionObject & err ) - { - std::cerr << "ExceptionObject caught while registering!" << std::endl; - std::cerr << err << std::endl; - return; - } - - - //======================================================= - // Get the result - //======================================================= - OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters(); - transform->SetParameters( finalParameters ); - - - //======================================================= - // Get the BSpline coefficient images and write them - //======================================================= - if (m_ArgsInfo.coeff_given) - { - if(m_ArgsInfo.wlut_flag) - { - BLUTTransformPointer bsplineTransform=dynamic_cast(registration->GetTransform()); - typedef itk::Image, ImageDimension> CoefficientImageType; - typename CoefficientImageType::Pointer coefficientImage =bsplineTransform->GetCoefficientImage(); - typedef itk::ImageFileWriter CoeffWriterType; - typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New(); - coeffWriter->SetInput(coefficientImage); - coeffWriter->SetFileName(m_ArgsInfo.coeff_arg[0]); - coeffWriter->Update(); - } - else - { - BSplineTransformPointer bsplineTransform=dynamic_cast(registration->GetTransform()); - typedef itk::Image CoefficientImageType; -#if ITK_VERSION_MAJOR > 3 - typename BSplineTransformType::CoefficientImageArray coefficientImages = bsplineTransform->GetCoefficientImage(); -#else - typename CoefficientImageType::Pointer *coefficientImages =bsplineTransform->GetCoefficientImage(); -#endif - typedef itk::ImageFileWriter CoeffWriterType; - for (unsigned int i=0;iSetInput(coefficientImages[i]); - coeffWriter->SetFileName(m_ArgsInfo.coeff_arg[i]); - coeffWriter->Update(); - } - } - } - - - //======================================================= - // Generate the DVF - //======================================================= - typedef itk::Vector< float, SpaceDimension > DisplacementType; - typedef itk::Image< DisplacementType, ImageDimension > DeformationFieldType; - - typename DeformationFieldType::Pointer field = DeformationFieldType::New(); - field->SetRegions( fixedImageRegion ); - field->SetOrigin( fixedImage->GetOrigin() ); - field->SetSpacing( fixedImage->GetSpacing() ); - field->SetDirection( fixedImage->GetDirection() ); - field->Allocate(); - - typedef itk::ImageRegionIteratorWithIndex< DeformationFieldType > FieldIterator; - FieldIterator fi( field, fixedImageRegion ); - fi.GoToBegin(); - - typename TransformType::InputPointType fixedPoint; - typename TransformType::OutputPointType movingPoint; - typename DeformationFieldType::IndexType index; - - DisplacementType displacement; - while( ! fi.IsAtEnd() ) - { - index = fi.GetIndex(); - field->TransformIndexToPhysicalPoint( index, fixedPoint ); - movingPoint = transform->TransformPoint( fixedPoint ); - displacement = movingPoint - fixedPoint; - fi.Set( displacement ); - ++fi; - } - - - //======================================================= - // Write the DVF - //======================================================= - typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType; - typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New(); - fieldWriter->SetFileName( m_ArgsInfo.vf_arg ); - fieldWriter->SetInput( field ); - try - { - fieldWriter->Update(); - } - catch( itk::ExceptionObject & excp ) - { - std::cerr << "Exception thrown writing the DVF" << std::endl; - std::cerr << excp << std::endl; - return; - } - - - //======================================================= - // Resample the moving image - //======================================================= - typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType > WarpFilterType; - typename WarpFilterType::Pointer warp = WarpFilterType::New(); - - warp->SetDeformationField( field ); - warp->SetInput( movingImageReader->GetOutput() ); - warp->SetOutputOrigin( fixedImage->GetOrigin() ); - warp->SetOutputSpacing( fixedImage->GetSpacing() ); - warp->SetOutputDirection( fixedImage->GetDirection() ); - warp->SetEdgePaddingValue( 0.0 ); - warp->Update(); - - - //======================================================= - // Write the warped image - //======================================================= - typedef itk::ImageFileWriter< FixedImageType > WriterType; - typename WriterType::Pointer writer = WriterType::New(); - writer->SetFileName( m_ArgsInfo.output_arg ); - writer->SetInput( warp->GetOutput() ); - - try - { - writer->Update(); - } - catch( itk::ExceptionObject & err ) - { - std::cerr << "ExceptionObject caught !" << std::endl; - std::cerr << err << std::endl; - return; - } - - - //======================================================= - // Calculate the difference after the deformable transform - //======================================================= - typedef clitk::DifferenceImageFilter< FixedImageType, FixedImageType> DifferenceFilterType; - if (m_ArgsInfo.after_given) - { - typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New(); - difference->SetValidInput( fixedImage ); - difference->SetTestInput( warp->GetOutput() ); - - try - { - difference->Update(); - } - catch( itk::ExceptionObject & err ) - { - std::cerr << "ExceptionObject caught calculating the difference !" << std::endl; - std::cerr << err << std::endl; - return; - } - - typename WriterType::Pointer differenceWriter=WriterType::New(); - differenceWriter->SetInput(difference->GetOutput()); - differenceWriter->SetFileName(m_ArgsInfo.after_arg); - differenceWriter->Update(); - - } - - - //======================================================= - // Calculate the difference before the deformable transform - //======================================================= - if( m_ArgsInfo.before_given ) - { - - typename FixedImageType::Pointer moving=FixedImageType::New(); - if (m_ArgsInfo.rigid_given) - { - typedef itk::ResampleImageFilter ResamplerType; - typename ResamplerType::Pointer resampler=ResamplerType::New(); - resampler->SetInput(movingImage); - resampler->SetOutputOrigin(fixedImage->GetOrigin()); - resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize()); - resampler->SetOutputSpacing(fixedImage->GetSpacing()); - resampler->SetDefaultPixelValue( 0. ); - //resampler->SetTransform(rigidTransform); - resampler->Update(); - moving=resampler->GetOutput(); - } - else - moving=movingImage; - - typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New(); - difference->SetValidInput( fixedImage ); - difference->SetTestInput( moving ); - - try - { - difference->Update(); - } - catch( itk::ExceptionObject & err ) - { - std::cerr << "ExceptionObject caught calculating the difference !" << std::endl; - std::cerr << err << std::endl; - return; - } - - typename WriterType::Pointer differenceWriter=WriterType::New(); - writer->SetFileName( m_ArgsInfo.before_arg ); - writer->SetInput( difference->GetOutput() ); - writer->Update( ); - } - - return; - } -} - -#endif // __clitkBSplineDeformableRegistrationGenericFilter_txx -- 2.45.1