From f9c73d25ecff9bfb90c9b040acab74dc102ceac2 Mon Sep 17 00:00:00 2001 From: Romulo Pinho Date: Mon, 6 Feb 2012 15:13:30 +0100 Subject: [PATCH] Convert coeffs to VF with itk bspline - conversion with BLUT has bugs (grid appears in VF's jacobian with anisotropic spacing) - disabled 4D convertion (spatio temporal) - conversion with BLUT still possible (type = 1, not default) - input coefficient image is always in BLUT style - input is converted to is converted to ITK style (one image per component) if type = 0 --- ...kConvertBSplineDeformableTransformToVF.ggo | 9 +- ...neDeformableTransformToVFGenericFilter.cxx | 464 ++++++++++-------- ...lineDeformableTransformToVFGenericFilter.h | 13 +- 3 files changed, 274 insertions(+), 212 deletions(-) diff --git a/registration/clitkConvertBSplineDeformableTransformToVF.ggo b/registration/clitkConvertBSplineDeformableTransformToVF.ggo index e20defb..9c715f3 100644 --- a/registration/clitkConvertBSplineDeformableTransformToVF.ggo +++ b/registration/clitkConvertBSplineDeformableTransformToVF.ggo @@ -8,7 +8,7 @@ option "verbose" v "Verbose" flag off section "IO" -option "input" i "Input coefficient image filename" string yes +option "input" i "Input BLUT-coefficient image filename" string yes option "output" o "Output image filename" string yes @@ -20,6 +20,7 @@ option "spacing" - "Spacing for the output image" double multiple no default section "Transform" -option "order" - "Spline order" int multiple no -option "mask" - "Mask image filename" string no -option "shape" - "Transform shape: 0=egg, 1=diamond" int no default="0" +option "type" t "Type (0: itk bspline; 1: BLUT)" int no default="0" +option "order" - "1: Spline order" int multiple no +option "mask" - "1: Mask image filename" string no +option "shape" - "1: Transform shape: 0=egg, 1=diamond" int no default="0" diff --git a/registration/clitkConvertBSplineDeformableTransformToVFGenericFilter.cxx b/registration/clitkConvertBSplineDeformableTransformToVFGenericFilter.cxx index 169de13..98e5f75 100644 --- a/registration/clitkConvertBSplineDeformableTransformToVFGenericFilter.cxx +++ b/registration/clitkConvertBSplineDeformableTransformToVFGenericFilter.cxx @@ -27,7 +27,24 @@ * ===================================================*/ +// clitk include +#include "clitkIO.h" +#include "clitkCommon.h" +#include "clitkImageCommon.h" +#include "clitkConvertBSplineDeformableTransformToVF_ggo.h" +#include "clitkBSplineDeformableTransform.h" +#include "clitkTransformToDeformationFieldSource.h" +#include "clitkShapedBLUTSpatioTemporalDeformableTransform.h" +#include "itkImageMaskSpatialObject.h" + #include "clitkConvertBSplineDeformableTransformToVFGenericFilter.h" +#include "clitkVectorImageToImageFilter.h" +#if ITK_VERSION_MAJOR >= 4 +#include "itkTransformToDisplacementFieldSource.h" +#else +#include "itkTransformToDeformationFieldSource.h" +#endif +#include "itkBSplineDeformableTransform.h" namespace clitk @@ -49,26 +66,26 @@ namespace clitk template<> void ConvertBSplineDeformableTransformToVFGenericFilter::UpdateWithDim<3>(std::string PixelType, int Components) - { +{ // Components if (Components !=3) - { - std::cerr<<"Number of components is "< InputPixelType; typedef itk::Vector OutputPixelType; typedef itk::Image InputImageType; typedef itk::Image OutputImageType; - + // Read the input typedef itk::ImageFileReader InputReaderType; InputReaderType::Pointer reader = InputReaderType::New(); @@ -89,105 +106,162 @@ namespace clitk //Output image info if (m_ArgsInfo.like_given) - { - typedef itk::ImageFileReader ReaderType; - ReaderType::Pointer reader2=ReaderType::New(); - reader2->SetFileName(m_ArgsInfo.like_arg); - reader2->Update(); - - OutputImageType::Pointer image=reader2->GetOutput(); - filter->SetOutputParametersFromImage(image); - } + { +/* typedef itk::ImageFileReader ReaderType; + ReaderType::Pointer reader2=ReaderType::New(); + reader2->SetFileName(m_ArgsInfo.like_arg); + reader2->Update(); + + OutputImageType::Pointer image=reader2->GetOutput(); + filter->SetOutputParametersFromImage(image);*/ + + typedef itk::ImageIOBase ImageIOType; + typename ImageIOType::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(m_ArgsInfo.like_arg, itk::ImageIOFactory::ReadMode); + imageIO->SetFileName(m_ArgsInfo.like_arg); + imageIO->ReadImageInformation(); + + typename ConvertorType::SizeType output_size; + typename ConvertorType::SpacingType output_spacing; + typename ConvertorType::OriginType output_origin; + typename ConvertorType::DirectionType output_direction; + for (unsigned int i = 0; i < Dimension; i++) { + output_size[i] = imageIO->GetDimensions(i); + output_spacing[i] = imageIO->GetSpacing(i); + output_origin[i] = imageIO->GetOrigin(i); + for (unsigned int j = 0; j < Dimension; j++) + output_direction[i][j] = imageIO->GetDirection(i)[j]; + } + + filter->SetOutputOrigin(output_origin); + filter->SetOutputSpacing(output_spacing); + filter->SetOutputSize(output_size); + filter->SetOutputDirection(output_direction); + } else - { - unsigned int i=0; - if(m_ArgsInfo.origin_given) - { - OutputImageType::PointType origin; - for(i=0;iSetOutputOrigin(origin); - } - if (m_ArgsInfo.spacing_given) - { - OutputImageType::SpacingType spacing; - for(i=0;iSetOutputSpacing(spacing); - } - if (m_ArgsInfo.spacing_given) - { - OutputImageType::SizeType size; - for(i=0;iSetOutputSize(size); - } - } + { + unsigned int i=0; + if (m_ArgsInfo.origin_given) + { + OutputImageType::PointType origin; + for (i=0;iSetOutputOrigin(origin); + } + if (m_ArgsInfo.spacing_given) + { + OutputImageType::SpacingType spacing; + for (i=0;iSetOutputSpacing(spacing); + } + if (m_ArgsInfo.spacing_given) + { + OutputImageType::SizeType size; + for (i=0;iSetOutputSize(size); + } + } if (m_Verbose) - { - std::cout<< "Setting output origin to "<GetOutputOrigin()<<"..."<GetOutputSpacing()<<"..."<GetOutputSize()<<"..."<GetOutputOrigin()<<"..."<GetOutputSpacing()<<"..."<GetOutputSize()<<"..."< TransformType; - TransformType::Pointer transform=TransformType::New(); - - // Spline orders: Default is cubic splines - InputImageType::RegionType::SizeType splineOrders ; - splineOrders.Fill(3); - if (m_ArgsInfo.order_given) - for(unsigned int i=0; i BLUTTransformType; + BLUTTransformType::Pointer blut_transform=BLUTTransformType::New(); + + typedef itk::BSplineDeformableTransform< double, Dimension, Dimension> ITKTransformType; + ITKTransformType::Pointer itk_transform=ITKTransformType::New(); + + typedef itk::Transform< double, Dimension, Dimension> GenericTransformType; + typename GenericTransformType::Pointer transform; + // Mask typedef itk::ImageMaskSpatialObject< Dimension > MaskType; MaskType::Pointer spatialObjectMask=NULL; if (m_ArgsInfo.mask_given) + { + typedef itk::Image< unsigned char, Dimension > ImageMaskType; + typedef itk::ImageFileReader< ImageMaskType > MaskReaderType; + 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 <<"Mask was read..." <SetImage( maskReader->GetOutput() ); + blut_transform->SetMask(spatialObjectMask); + } + + if (m_ArgsInfo.type_arg != 0 ) { // using BLUT + // Spline orders: Default is cubic splines + if (m_Verbose) { + std::cout << "Using clitk::BLUT." << std::endl; + std::cout << "Setting spline orders and sampling factors." << std::endl; + } + InputImageType::RegionType::SizeType splineOrders ; + splineOrders.Fill(3); + if (m_ArgsInfo.order_given) + for (unsigned int i=0; i ImageMaskType; - typedef itk::ImageFileReader< ImageMaskType > MaskReaderType; - 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 <<"Mask was read..." <SetImage( maskReader->GetOutput() ); + samplingFactors[i]= (int) ( input->GetSpacing()[i]/ filter->GetOutputSpacing()[i]); + if (m_Verbose) std::cout<<"Setting sampling factor "<GetSpacing()[i]/ filter->GetOutputSpacing()[i]); - if (m_Verbose) std::cout<<"Setting sampling factor "<SetSplineOrders(splineOrders); - transform->SetMask(spatialObjectMask); - transform->SetLUTSamplingFactors(samplingFactors); - transform->SetCoefficientImage(input); + blut_transform->SetSplineOrders(splineOrders); + blut_transform->SetLUTSamplingFactors(samplingFactors); + blut_transform->SetCoefficientImage(input); + + transform = blut_transform; + } + else { // using ITK transform + if (m_Verbose) { + std::cout << "Using itk::BSpline" << std::endl; + std::cout << "Extracting components from input coefficient image and creating one coefficient image per-component" << std::endl; + } + + typedef double PixelType; + typedef itk::Vector CoefficientType; + typedef itk::Image CoefficientImageType; + typedef clitk::VectorImageToImageFilter FilterType; + typename FilterType::Pointer component_filter[Dimension]; + + typename ITKTransformType::ImagePointer coefficient_images[Dimension]; + for (unsigned int i=0; i < Dimension; i++) { + component_filter[i] = FilterType::New(); + component_filter[i]->SetInput(input); + component_filter[i]->SetComponentIndex(i); + component_filter[i]->Update(); + coefficient_images[i] = component_filter[i]->GetOutput(); + } + itk_transform->SetCoefficientImage(coefficient_images); + + transform = itk_transform; + } + filter->SetTransform(transform); @@ -196,13 +270,13 @@ namespace clitk // ----------------------------------------------- if (m_Verbose)std::cout<< "Converting the BSpline transform..."<Update(); - } + { + filter->Update(); + } catch (itk::ExceptionObject) - { - std::cerr<<"Error: Exception thrown during execution convertion filter!"<GetOutput(); @@ -218,26 +292,26 @@ namespace clitk } - +/* //------------------------------------------------------------------- // Update with the number of dimensions //------------------------------------------------------------------- template<> - void - ConvertBSplineDeformableTransformToVFGenericFilter::UpdateWithDim<4>(std::string PixelType, int Components) - { +void +ConvertBSplineDeformableTransformToVFGenericFilter::UpdateWithDim<4>(std::string PixelType, int Components) +{ // Components if (Components !=3) - { - std::cerr<<"Number of components is "< OutputPixelType; typedef itk::Image InputImageType; typedef itk::Image OutputImageType; - + // Read the input typedef itk::ImageFileReader InputReaderType; InputReaderType::Pointer reader = InputReaderType::New(); @@ -262,53 +336,53 @@ namespace clitk //Output image info if (m_ArgsInfo.like_given) - { - typedef itk::ImageFileReader ReaderType; - ReaderType::Pointer reader2=ReaderType::New(); - reader2->SetFileName(m_ArgsInfo.like_arg); - reader2->Update(); - - OutputImageType::Pointer image=reader2->GetOutput(); - filter->SetOutputParametersFromImage(image); - } + { + typedef itk::ImageFileReader ReaderType; + ReaderType::Pointer reader2=ReaderType::New(); + reader2->SetFileName(m_ArgsInfo.like_arg); + reader2->Update(); + + OutputImageType::Pointer image=reader2->GetOutput(); + filter->SetOutputParametersFromImage(image); + } else - { - unsigned int i=0; - if(m_ArgsInfo.origin_given) - { - OutputImageType::PointType origin; - for(i=0;iSetOutputOrigin(origin); - } - if (m_ArgsInfo.spacing_given) - { - OutputImageType::SpacingType spacing; - for(i=0;iSetOutputSpacing(spacing); - } - if (m_ArgsInfo.spacing_given) - { - OutputImageType::SizeType size; - for(i=0;iSetOutputSize(size); - } - } + { + unsigned int i=0; + if (m_ArgsInfo.origin_given) + { + OutputImageType::PointType origin; + for (i=0;iSetOutputOrigin(origin); + } + if (m_ArgsInfo.spacing_given) + { + OutputImageType::SpacingType spacing; + for (i=0;iSetOutputSpacing(spacing); + } + if (m_ArgsInfo.spacing_given) + { + OutputImageType::SizeType size; + for (i=0;iSetOutputSize(size); + } + } //Output image info if (m_Verbose) - { - std::cout<< "Setting output origin to "<GetOutputOrigin()<<"..."<GetOutputSpacing()<<"..."<GetOutputSize()<<"..."<GetOutputOrigin()<<"..."<GetOutputSpacing()<<"..."<GetOutputSize()<<"..."< TransformType; + typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform< double, Dimension, Dimension > TransformType; TransformType::Pointer transform=TransformType::New(); transform->SetTransformShape(m_ArgsInfo.shape_arg); @@ -316,46 +390,46 @@ namespace clitk InputImageType::RegionType::SizeType splineOrders ; splineOrders.Fill(3); if (m_ArgsInfo.order_given) - for(unsigned int i=0; i MaskType; MaskType::Pointer spatialObjectMask=NULL; if (m_ArgsInfo.mask_given) - { - typedef itk::Image< unsigned char, Dimension > ImageMaskType; - typedef itk::ImageFileReader< ImageMaskType > MaskReaderType; - 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 <<"Mask was read..." <SetImage( maskReader->GetOutput() ); - } + { + typedef itk::Image< unsigned char, Dimension > ImageMaskType; + typedef itk::ImageFileReader< ImageMaskType > MaskReaderType; + 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 <<"Mask was read..." <SetImage( maskReader->GetOutput() ); + } // Samplingfactors - InputImageType::SizeType samplingFactors; + InputImageType::SizeType samplingFactors; for (unsigned int i=0; i< Dimension; i++) - { - samplingFactors[i]= (int) ( input->GetSpacing()[i]/ filter->GetOutputSpacing()[i]); - if (m_Verbose) std::cout<<"Setting sampling factor "<GetSpacing()[i]/ filter->GetOutputSpacing()[i]); + if (m_Verbose) std::cout<<"Setting sampling factor "<SetSplineOrders(splineOrders); @@ -370,13 +444,13 @@ namespace clitk // ----------------------------------------------- if (m_Verbose)std::cout<< "Converting the BSpline transform..."<Update(); - } + { + filter->Update(); + } catch (itk::ExceptionObject) - { - std::cerr<<"Error: Exception thrown during execution convertion filter!"<GetOutput(); @@ -390,10 +464,8 @@ namespace clitk writer->SetInput(output); writer->Update(); - } - - - +} +*/ //----------------------------------------------------------- // Update //----------------------------------------------------------- @@ -408,7 +480,7 @@ namespace clitk // Call UpdateWithDim //if(Dimension==2) UpdateWithDim<2>(PixelType, Components); if(Dimension==3) UpdateWithDim<3>(PixelType, Components); - else if (Dimension==4) UpdateWithDim<4>(PixelType, Components); + //else if (Dimension==4) UpdateWithDim<4>(PixelType, Components); else { std::cout<<"Error, Only for 3 Dimensions!!!"<= 4 - #include "itkTransformToDisplacementFieldSource.h" -#else - #include "itkTransformToDeformationFieldSource.h" -#endif namespace clitk -- 2.47.1