From 46e15948191b583a7333683663c3975a279274a6 Mon Sep 17 00:00:00 2001 From: delmon Date: Tue, 11 Jan 2011 10:25:42 +0000 Subject: [PATCH] Add a new Transformation that wraps multiple BSplineDeformableTransform to be able to optimize them at the same time. --- .../clitkMultipleBSplineDeformableTransform.h | 276 ++++++++++ ...litkMultipleBSplineDeformableTransform.txx | 487 ++++++++++++++++++ ...pleBSplineDeformableTransformInitializer.h | 173 +++++++ ...eBSplineDeformableTransformInitializer.txx | 260 ++++++++++ 4 files changed, 1196 insertions(+) create mode 100644 registration/clitkMultipleBSplineDeformableTransform.h create mode 100644 registration/clitkMultipleBSplineDeformableTransform.txx create mode 100644 registration/clitkMultipleBSplineDeformableTransformInitializer.h create mode 100644 registration/clitkMultipleBSplineDeformableTransformInitializer.txx diff --git a/registration/clitkMultipleBSplineDeformableTransform.h b/registration/clitkMultipleBSplineDeformableTransform.h new file mode 100644 index 0000000..93277c3 --- /dev/null +++ b/registration/clitkMultipleBSplineDeformableTransform.h @@ -0,0 +1,276 @@ +/*========================================================================= + 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://oncora1.lyon.fnclcc.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 __clitkMultipleBSplineDeformableTransform_h +#define __clitkMultipleBSplineDeformableTransform_h + +#include + +#include + +#include + +namespace clitk +{ + template < + class TCoordRep = double, + unsigned int NInputDimensions = 3, + unsigned int NOutputDimensions = 3 > + class ITK_EXPORT MultipleBSplineDeformableTransform : public itk::Transform< TCoordRep, NInputDimensions, NOutputDimensions > + { + public: + + //==================================================================== + // Typedefs + //==================================================================== + typedef MultipleBSplineDeformableTransform Self; + typedef itk::Transform Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** New macro for creation of through the object factory.*/ + itkNewMacro( Self ); + + /** Run-time type information (and related methods). */ + itkTypeMacro( MultipleBSplineDeformableTransform, Transform ); + + /** Dimension of the domain space. */ + itkStaticConstMacro(OutputDimension, unsigned int, NOutputDimensions); + + /** Dimension of the input model. */ + itkStaticConstMacro(InputDimension, unsigned int, NInputDimensions); + + /** Standard scalar type for this class. */ + typedef typename Superclass::ScalarType ScalarType; + + /** Standard parameters container. */ + typedef typename Superclass::ParametersType ParametersType; + + /** Standard Jacobian container. */ + typedef typename Superclass::JacobianType JacobianType; + + /** Standard vector type for this class. */ + typedef itk::Vector InputVectorType; + typedef itk::Vector OutputVectorType; + + /** Standard covariant vector type for this class. */ + typedef itk::CovariantVector InputCovariantVectorType; + typedef itk::CovariantVector OutputCovariantVectorType; + + /** Standard vnl_vector type for this class. */ + typedef vnl_vector_fixed InputVnlVectorType; + typedef vnl_vector_fixed OutputVnlVectorType; + + /** Standard coordinate point type for this class. */ + typedef itk::Point InputPointType; + typedef itk::Point OutputPointType; + + //JV Parameters as images with OutputDimension number of components per Pixel + typedef typename ParametersType::ValueType ParametersValueType; + typedef typename itk::Vector PixelType; + typedef itk::Image CoefficientImageType; + typedef typename CoefficientImageType::Pointer CoefficientImagePointer; + + /** Typedefs for specifying the extend to the grid. */ + typedef itk::ImageRegion RegionType; + typedef typename RegionType::IndexType IndexType; + typedef typename RegionType::SizeType SizeType; + typedef typename CoefficientImageType::SpacingType SpacingType; + typedef typename CoefficientImageType::DirectionType DirectionType; + typedef typename CoefficientImageType::PointType OriginType; + typedef itk::ContinuousIndex ContinuousIndexType; + + /** Typedef of the bulk transform. */ + typedef itk::Transform BulkTransformType; + typedef BulkTransformType* BulkTransformPointer; + + typedef itk::Image< unsigned char, InputDimension> ImageLabelType; + typedef ImageLabelType* ImageLabelPointer; + //typedef typename ImageLabelType::Pointer ImageLabelPointer; + + typedef itk::NearestNeighborInterpolateImageFunction ImageLabelInterpolator; + typedef typename ImageLabelInterpolator::Pointer ImageLabelInterpolatorPointer; + + void SetLabels(ImageLabelPointer labels); + ImageLabelPointer GetLabels() {return m_labels;} + + itkGetMacro(nLabels, unsigned); + + void SetSplineOrder(const unsigned int & splineOrder); + void SetSplineOrders(const SizeType & splineOrders); + /* + itkGetMacro( SplineOrders, SizeType ); + itkGetConstMacro( SplineOrders, SizeType ); + */ + + void SetLUTSamplingFactor(const int & samplingFactor); + void SetLUTSamplingFactors(const SizeType & samplingFactors); + /* + itkGetMacro( LUTSamplingFactors, SizeType ); + itkGetConstMacro( LUTSamplingFactors,SizeType ); + */ + + void SetParameters(const ParametersType & parameters); + + void SetFixedParameters(const ParametersType & parameters); + + void SetParametersByValue(const ParametersType & parameters); + + void SetIdentity(); + + /** Get the Transformation Parameters. */ + virtual const ParametersType& GetParameters(void) const; + + /** Get the Transformation Fixed Parameters. */ + virtual const ParametersType& GetFixedParameters(void) const; + + // The coefficientImage + virtual const std::vector& GetCoefficientImages() const; + virtual void SetCoefficientImages(std::vector& images); + + /** This method specifies the region over which the grid resides. */ + virtual void SetGridRegion( const RegionType& region ); + /* + itkGetMacro( GridRegion, RegionType ); + itkGetConstMacro( GridRegion, RegionType ); + */ + + /** This method specifies the grid spacing or resolution. */ + virtual void SetGridSpacing( const SpacingType& spacing ); + /* + itkGetMacro( GridSpacing, SpacingType ); + itkGetConstMacro( GridSpacing, SpacingType ); + */ + + /** This method specifies the grid directions . */ + virtual void SetGridDirection( const DirectionType & spacing ); + /* + itkGetMacro( GridDirection, DirectionType ); + itkGetConstMacro( GridDirection, DirectionType ); + */ + + /** This method specifies the grid origin. */ + virtual void SetGridOrigin( const OriginType& origin ); + /* + itkGetMacro( GridOrigin, OriginType ); + itkGetConstMacro( GridOrigin, OriginType ); + */ + + // Set the bulk transform, real pointer + void SetBulkTransform(BulkTransformPointer b); + BulkTransformPointer GetBulkTransform(void) {return m_BulkTransform;} + + /** Transform points by a BSpline deformable transformation. */ + OutputPointType TransformPoint(const InputPointType &point ) const; + + // JV added for just the deformable part, without bulk + OutputPointType DeformablyTransformPoint(const InputPointType &point ) const; + + /** Method to transform a vector - + * not applicable for this type of transform. */ + virtual OutputVectorType TransformVector(const InputVectorType &) const + { + itkExceptionMacro(<< "Method not applicable for deformable transform." ); + return OutputVectorType(); + } + + /** Method to transform a vnl_vector - + * not applicable for this type of transform */ + virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const + { + itkExceptionMacro(<< "Method not applicable for deformable transform. "); + return OutputVnlVectorType(); + } + + /** Method to transform a CovariantVector - + * not applicable for this type of transform */ + virtual OutputCovariantVectorType TransformCovariantVector(const InputCovariantVectorType &) const + { + itkExceptionMacro(<< "Method not applicable for deformable transfrom. "); + return OutputCovariantVectorType(); + } + + /** Compute the Jacobian Matrix of the transformation at one point */ + virtual const JacobianType& GetJacobian(const InputPointType &point ) const; + + /** Return the number of parameters that completely define the Transfom */ + virtual unsigned int GetNumberOfParameters(void) const; + + /** Return the number of parameters per dimension */ + unsigned int GetNumberOfParametersPerDimension(void) const; + + virtual bool IsLinear() const { return false; } + + typedef clitk::BSplineDeformableTransform TransformType; + + protected: + + void PrintSelf(std::ostream &os, itk::Indent indent) const; + + MultipleBSplineDeformableTransform(); + virtual ~MultipleBSplineDeformableTransform(); + + void TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const; + + private: + MultipleBSplineDeformableTransform(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + /** The bulk transform. */ + BulkTransformPointer m_BulkTransform; + + /** Jacobian as OutputDimension number of images. */ + typedef typename JacobianType::ValueType JacobianValueType; + typedef typename itk::Vector JacobianPixelType; + typedef itk::Image JacobianImageType; + typename JacobianImageType::Pointer m_JacobianImage[OutputDimension]; + typedef itk::ImageRegionIterator IteratorType; + + /** Keep a pointer to the input parameters. */ + const ParametersType * m_InputParametersPointer; + + /** Internal parameters buffer. */ + ParametersType m_InternalParametersBuffer; + + unsigned int m_nLabels; + ImageLabelPointer m_labels; + ImageLabelInterpolatorPointer m_labelInterpolator; + std::vector m_trans; + std::vector m_parameters; + mutable std::vector m_CoefficientImages; + mutable int m_LastJacobian; + + void InitJacobian(); + // FIXME it seems not used + bool InsideValidRegion( const ContinuousIndexType& index ) const; + }; + +} // namespace clitk + +#if ITK_TEMPLATE_TXX +# include "clitkMultipleBSplineDeformableTransform.txx" +#endif + +#endif // __clitkMultipleBSplineDeformableTransform_h diff --git a/registration/clitkMultipleBSplineDeformableTransform.txx b/registration/clitkMultipleBSplineDeformableTransform.txx new file mode 100644 index 0000000..3a609a2 --- /dev/null +++ b/registration/clitkMultipleBSplineDeformableTransform.txx @@ -0,0 +1,487 @@ +/*========================================================================= + 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://oncora1.lyon.fnclcc.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 __clitkMultipleBSplineDeformableTransform_txx +#define __clitkMultipleBSplineDeformableTransform_txx + +// clitk +#include "clitkMultipleBSplineDeformableTransformInitializer.h" + +// itk +#include + +namespace clitk +{ + // Constructor with default arguments + template + MultipleBSplineDeformableTransform + ::MultipleBSplineDeformableTransform() : Superclass(OutputDimension, 0) + { + m_nLabels = 1; + m_labels = 0; + m_labelInterpolator = 0; + m_trans.resize(1); + m_trans[0] = TransformType::New(); + m_parameters.resize(1); + + this->m_FixedParameters.SetSize(NInputDimensions * (NInputDimensions + 5) + 4); + this->m_FixedParameters.Fill(0.0); + + m_InternalParametersBuffer = ParametersType(0); + // Make sure the parameters pointer is not NULL after construction. + m_InputParametersPointer = &m_InternalParametersBuffer; + m_BulkTransform = 0; + m_LastJacobian = -1; + } + + template + MultipleBSplineDeformableTransform + ::~MultipleBSplineDeformableTransform() + { + } + + template + void + MultipleBSplineDeformableTransform + ::SetLabels(ImageLabelPointer labels) + { + GetFixedParameters(); // Update m_FixedParameters + if (labels) + { + typedef itk::RelabelComponentImageFilter RelabelImageFilterType; + typename RelabelImageFilterType::Pointer relabelImageFilter = RelabelImageFilterType::New(); + relabelImageFilter->SetInput(labels); + relabelImageFilter->Update(); + m_labels = relabelImageFilter->GetOutput(); + m_nLabels = relabelImageFilter->GetNumberOfObjects(); + m_trans.resize(m_nLabels); + m_parameters.resize(m_nLabels); + for (unsigned i = 0; i < m_nLabels; ++i) + m_trans[i] = TransformType::New(); + m_labelInterpolator = ImageLabelInterpolator::New(); + m_labelInterpolator->SetInputImage(m_labels); + } + else + { + m_nLabels = 1; + m_trans.resize(1); + m_parameters.resize(1); + m_trans[0] = TransformType::New(); + } + this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions + 2] = (double)((size_t) m_labels); + this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions + 3] = m_nLabels; + SetFixedParameters(this->m_FixedParameters); + } + +#define LOOP_ON_LABELS(FUNC, ARGS) \ + for (unsigned i = 0; i < m_nLabels; ++i) \ + m_trans[i]->FUNC(ARGS); + + template + inline void + MultipleBSplineDeformableTransform + ::SetSplineOrder(const unsigned int & splineOrder) + { + LOOP_ON_LABELS(SetSplineOrder, splineOrder); + this->Modified(); + } + + template + inline void + MultipleBSplineDeformableTransform + ::SetSplineOrders(const SizeType & splineOrders) + { + LOOP_ON_LABELS(SetSplineOrders, splineOrders); + this->Modified(); + } + + template + inline void + MultipleBSplineDeformableTransform + ::SetLUTSamplingFactor( const int & samplingFactor) + { + LOOP_ON_LABELS(SetLUTSamplingFactor, samplingFactor); + this->Modified(); + } + + template + inline void + MultipleBSplineDeformableTransform + ::SetLUTSamplingFactors( const SizeType & samplingFactors) + { + LOOP_ON_LABELS(SetLUTSamplingFactors, samplingFactors); + this->Modified(); + } + + template + inline void + MultipleBSplineDeformableTransform + ::SetGridRegion( const RegionType & region ) + { + LOOP_ON_LABELS(SetGridRegion, region); + this->Modified(); + } + + template + inline void + MultipleBSplineDeformableTransform + ::SetGridSpacing( const SpacingType & spacing ) + { + LOOP_ON_LABELS(SetGridSpacing, spacing); + this->Modified(); + } + + template + inline void + MultipleBSplineDeformableTransform + ::SetGridDirection( const DirectionType & direction ) + { + LOOP_ON_LABELS(SetGridDirection, direction); + this->Modified(); + } + + template + inline void + MultipleBSplineDeformableTransform + ::SetGridOrigin( const OriginType& origin ) + { + LOOP_ON_LABELS(SetGridOrigin, origin); + this->Modified(); + } + + template + inline void + MultipleBSplineDeformableTransform + ::SetIdentity() + { + LOOP_ON_LABELS(SetIdentity, ); + if ( m_InputParametersPointer ) + { + ParametersType * parameters = + const_cast( m_InputParametersPointer ); + parameters->Fill( 0.0 ); + } + this->Modified(); + } + + template + inline const std::vector::CoefficientImagePointer>& + MultipleBSplineDeformableTransform + ::GetCoefficientImages() const + { + m_CoefficientImages.resize(m_nLabels); + for (unsigned i = 0; i < m_nLabels; ++i) + m_CoefficientImages[i] = m_trans[i]->GetCoefficientImage(); + return m_CoefficientImages; + } + + template + inline void + MultipleBSplineDeformableTransform + ::SetCoefficientImages(std::vector& images) + { + if (images.size() != m_nLabels) + { + itkExceptionMacro( << "Mismatched between the number of labels " + << m_nLabels + << " and the number of coefficient images " + << images.size()); + } + for (unsigned i = 0; i < m_nLabels; ++i) + m_trans[i]->SetCoefficientImage(images[i]); + } + + template + void + MultipleBSplineDeformableTransform + ::SetParameters( const ParametersType & parameters ) + { + if ( parameters.Size() != this->GetNumberOfParameters() ) + { + itkExceptionMacro(<<"Mismatched between parameters size " + << parameters.size() + << " and region size " + << this->GetNumberOfParameters() ); + } + + // Clean up buffered parameters + m_InternalParametersBuffer = ParametersType( 0 ); + + // Keep a reference to the input parameters + m_InputParametersPointer = ¶meters; + + double * dataPointer = const_cast(m_InputParametersPointer->data_block()); + unsigned tot = 0; + for (unsigned i = 0; i < m_nLabels; ++i) + { + unsigned int numberOfParameters = m_trans[i]->GetNumberOfParameters(); + m_parameters[i].SetData(dataPointer + tot, numberOfParameters); + m_trans[i]->SetParameters(m_parameters[i]); + tot += numberOfParameters; + } + InitJacobian(); + + this->Modified(); + } + + template + void + MultipleBSplineDeformableTransform + ::SetFixedParameters( const ParametersType & parameters ) + { + // BSplineSeformableTransform parameters + m_labels pointer + if ( parameters.Size() != NInputDimensions * (5 + NInputDimensions) + 4) + { + itkExceptionMacro(<< "Mismatched between parameters size " + << parameters.size() + << " and number of fixed parameters " + << NInputDimensions * (5 + NInputDimensions) + 4); + } + ImageLabelPointer tmp2 = ((ImageLabelPointer)( (size_t)parameters[(5+NInputDimensions)*NInputDimensions+2])); + if (tmp2 != m_labels) + { + if (tmp2) + { + m_labels = tmp2; + m_nLabels = parameters[(5+NInputDimensions) * NInputDimensions + 3 ]; + m_trans.resize(m_nLabels); + m_parameters.resize(m_nLabels); + for (unsigned i = 0; i < m_nLabels; ++i) + m_trans[i] = TransformType::New(); + m_labelInterpolator = ImageLabelInterpolator::New(); + m_labelInterpolator->SetInputImage(m_labels); + } + else + { + m_labels = tmp2; + m_nLabels = 1; + m_trans.resize(1); + m_parameters.resize(1); + m_trans[0] = TransformType::New(); + } + } + unsigned int numberOfParameters = NInputDimensions * (5 + NInputDimensions) + 2; + ParametersType tmp(numberOfParameters); + for (unsigned j = 0; j < numberOfParameters; ++j) + tmp.SetElement(j, parameters.GetElement(j)); + LOOP_ON_LABELS(SetFixedParameters, tmp); + this->m_FixedParameters = parameters; + this->Modified(); + } + + template + void + MultipleBSplineDeformableTransform + ::SetParametersByValue( const ParametersType & parameters ) + { + if ( parameters.Size() != this->GetNumberOfParameters() ) + { + itkExceptionMacro(<<"Mismatched between parameters size " + << parameters.size() + << " and region size " + << this->GetNumberOfParameters() ); + } + + // Copy it + m_InternalParametersBuffer = parameters; + m_InputParametersPointer = &m_InternalParametersBuffer; + + for (unsigned i = 0, tot = 0; i < m_nLabels; ++i) + { + unsigned int numberOfParameters = m_trans[i]->GetNumberOfParameters(); + ParametersType tmp(numberOfParameters); + for (unsigned j = 0; j < numberOfParameters; ++j, ++tot) + tmp.SetElement(j, parameters.GetElement(tot)); + m_trans[i]->SetParametersByValue(tmp); + } + InitJacobian(); + this->Modified(); + } + + + template + inline void + MultipleBSplineDeformableTransform + ::SetBulkTransform(BulkTransformPointer b) + { + m_BulkTransform = b; + LOOP_ON_LABELS(SetBulkTransform, b); + } + +#undef LOOP_ON_LABELS + + template + inline unsigned int + MultipleBSplineDeformableTransform + ::GetNumberOfParameters(void) const + { + unsigned int sum = 0; + for (unsigned i = 0; i < m_nLabels; ++i) + sum += m_trans[i]->GetNumberOfParameters(); + return sum; + } + + template + inline unsigned int + MultipleBSplineDeformableTransform + ::GetNumberOfParametersPerDimension(void) const + { + unsigned int sum = 0; + for (unsigned i = 0; i < m_nLabels; ++i) + sum += m_trans[i]->GetNumberOfParametersPerDimension(); + return sum; + } + + template + inline const typename MultipleBSplineDeformableTransform::ParametersType & + MultipleBSplineDeformableTransform + ::GetParameters() const + { + if (NULL == m_InputParametersPointer) + { + itkExceptionMacro( <<"Cannot GetParameters() because m_InputParametersPointer is NULL. Perhaps SetCoefficientImage() has been called causing the NULL pointer." ); + } + + return (*m_InputParametersPointer); + } + + template + inline const typename MultipleBSplineDeformableTransform::ParametersType & + MultipleBSplineDeformableTransform + ::GetFixedParameters() const + { + const ParametersType& trans0Para = m_trans[0]->GetFixedParameters(); + for (unsigned i = 0; i < trans0Para.Size() ; ++i) + this->m_FixedParameters.SetElement(i, trans0Para.GetElement(i)); + this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions + 2] = (double)((size_t) m_labels); + this->m_FixedParameters[(5+NInputDimensions)*NInputDimensions + 3] = m_nLabels; + return this->m_FixedParameters; + } + + template + inline void + MultipleBSplineDeformableTransform + ::PrintSelf(std::ostream &os, itk::Indent indent) const + { + this->Superclass::PrintSelf(os, indent); + os << indent << "nLabels" << m_nLabels << std::endl; + itk::Indent ind = indent.GetNextIndent(); + + for (unsigned i = 0; i < m_nLabels; ++i) + { + os << indent << "Label " << i << std::endl; + m_trans[i]->Print(os, ind); + } + } + + template + inline bool + MultipleBSplineDeformableTransform + ::InsideValidRegion( const ContinuousIndexType& index ) const + { + bool res = false; + for (unsigned i = 0; i < m_nLabels; ++i) + res |= m_trans[i]->InsideValidRegion(index); + return res; + } + + template + inline typename MultipleBSplineDeformableTransform::OutputPointType + MultipleBSplineDeformableTransform + ::TransformPoint(const InputPointType &inputPoint) const + { + int lidx = 0; + if (m_labels) + lidx = m_labelInterpolator->Evaluate(inputPoint) - 1; + if (lidx == -1) + return inputPoint; + return m_trans[lidx]->TransformPoint(inputPoint); + } + + template + inline typename MultipleBSplineDeformableTransform::OutputPointType + MultipleBSplineDeformableTransform + ::DeformablyTransformPoint(const InputPointType &inputPoint) const + { + int lidx = 0; + if (m_labels) + lidx = m_labelInterpolator->Evaluate(inputPoint) - 1; + if (lidx == -1) + return inputPoint; + return m_trans[lidx]->DeformablyTransformPoint(inputPoint); + } + + template + inline const typename MultipleBSplineDeformableTransform::JacobianType & + MultipleBSplineDeformableTransform + ::GetJacobian( const InputPointType & point ) const + { + /* + for (unsigned i = 0; i < m_nLabels; ++i) + m_trans[i]->ResetJacobian(); + */ + if (m_LastJacobian != -1) + m_trans[m_LastJacobian]->ResetJacobian(); + + int lidx = 0; + if (m_labels) + lidx = m_labelInterpolator->Evaluate(point) - 1; + if (lidx == -1) + { + m_LastJacobian = lidx; + return this->m_Jacobian; + } + + m_trans[lidx]->GetJacobian(point); + m_LastJacobian = lidx; + + return this->m_Jacobian; + } + + template + inline void + MultipleBSplineDeformableTransform + ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const + { + m_trans[0]->TransformPointToContinuousIndex(point, index); + } + + template + inline void + MultipleBSplineDeformableTransform::InitJacobian() + { + unsigned numberOfParameters = this->GetNumberOfParameters(); + this->m_Jacobian.set_size(OutputDimension, numberOfParameters); + JacobianPixelType * jacobianDataPointer = reinterpret_cast(this->m_Jacobian.data_block()); + memset(jacobianDataPointer, 0, numberOfParameters * sizeof (JacobianPixelType)); + + unsigned tot = 0; + for (unsigned int j = 0; j < OutputDimension; j++) + { + for (unsigned i = 0; i < m_nLabels; ++i) + { + unsigned numberOfPixels = m_trans[i]->SetJacobianImageData(jacobianDataPointer, j); + jacobianDataPointer += numberOfPixels; + tot += numberOfPixels; + } + } + assert(tot == numberOfParameters); + } + +} // namespace clitk + +#endif // __clitkMultipleBSplineDeformableTransform_txx diff --git a/registration/clitkMultipleBSplineDeformableTransformInitializer.h b/registration/clitkMultipleBSplineDeformableTransformInitializer.h new file mode 100644 index 0000000..63a30c7 --- /dev/null +++ b/registration/clitkMultipleBSplineDeformableTransformInitializer.h @@ -0,0 +1,173 @@ +/*========================================================================= + 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://oncora1.lyon.fnclcc.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 __clitkMultipleBSplineDeformableTransformInitializer_h +#define __clitkMultipleBSplineDeformableTransformInitializer_h +#include "clitkResampleBSplineDeformableTransformImageFilter.h" + +#include "itkObject.h" +#include "itkObjectFactory.h" +#include + +namespace clitk +{ + + template < class TTransform, class TImage > + class ITK_EXPORT MultipleBSplineDeformableTransformInitializer : public itk::Object + { + public: + /** Standard class typedefs. */ + typedef MultipleBSplineDeformableTransformInitializer Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** New macro for creation of through a Smart Pointer. */ + itkNewMacro( Self ); + + /** Run-time type information (and related methods). */ + itkTypeMacro( MultipleBSplineDeformableTransformInitializer, Object ); + + // Typedefs + typedef TTransform TransformType; + typedef typename TransformType::Pointer TransformPointer; + typedef typename TransformType::RegionType RegionType; + typedef typename RegionType::SizeType SizeType; + typedef typename TransformType::SpacingType SpacingType; + itkStaticConstMacro(InputDimension, unsigned int, TransformType::InputDimension); + typedef TImage ImageType; + typedef typename ImageType::ConstPointer ImagePointer; + typedef typename TransformType::CoefficientImageType CoefficientImageType; + typedef typename TransformType::ParametersType ParametersType; + + // Set and Get + itkBooleanMacro(Verbose); + itkSetMacro( Verbose, bool); + itkGetConstReferenceMacro( Verbose, bool); + itkSetObjectMacro( Transform, TransformType ); + itkGetConstObjectMacro( Transform, TransformType ); + itkSetObjectMacro( Image, ImageType ); + itkGetConstObjectMacro( Image, ImageType ); + void SetSplineOrder(const unsigned int & splineOrder) + { + SizeType s; + s.Fill(splineOrder); + this->SetSplineOrders(s); + } + void SetSplineOrders(const SizeType & splineOrders) + { + m_SplineOrders=splineOrders; + } + void SetNumberOfControlPointsInsideTheImage( SizeType & n ) + { + m_NumberOfControlPointsInsideTheImage=n; + m_NumberOfControlPointsIsGiven=true; + this->Modified(); + } + void SetNumberOfControlPointsInsideTheImage( int * n) + { + SizeType s; + for (unsigned int i=0;iSetNumberOfControlPointsInsideTheImage( s ); + } + void SetNumberOfControlPointsInsideTheImage( unsigned int & n ) + { + SizeType s; + s.Fill( n );; + this->SetNumberOfControlPointsInsideTheImage( s ); + } + void SetControlPointSpacing( SpacingType n ) + { + m_ControlPointSpacing= n; + m_ControlPointSpacingIsGiven=true; + this->Modified(); + } + void SetControlPointSpacing( double*& n ) + { + SpacingType s( n ); + this->SetControlPointSpacing(s); + } + void SetControlPointSpacing( double n ) + { + SpacingType s; + s.Fill( n ); + this->SetControlPointSpacing(s); + } + void SetSamplingFactors( SizeType n ) + { + m_SamplingFactors=n; + m_SamplingFactorIsGiven=true; + this->Modified(); + } + void SetSamplingFactors( int *& n) + { + SizeType s; + for (unsigned int i=0;i SetSamplingFactors( s ); + } + void SetSamplingFactors( unsigned int n ) + { + SizeType s; + s.Fill( n ); + this-> SetSamplingFactors( s ); + } + virtual void InitializeTransform(); + void SetInitialParameters(std::vector coefficientImage, ParametersType& params); + void SetInitialParameters(const std::string & s, ParametersType& params); + + // For easy acces, declared public + std::vector m_NumberOfControlPointsInsideTheImageArray; + std::vector m_SamplingFactorsArray; + std::vector m_ControlPointSpacingArray; + + SpacingType m_ControlPointSpacing; + SizeType m_SamplingFactors; + SizeType m_SplineOrders; + SpacingType m_ChosenSpacing; + SizeType m_NumberOfControlPointsInsideTheImage; + bool m_NumberOfControlPointsIsGiven; + bool m_ControlPointSpacingIsGiven; + bool m_SamplingFactorIsGiven; + bool m_TransformRegionIsGiven; + + typename TransformType::ParametersType* m_Parameters; + + protected: + MultipleBSplineDeformableTransformInitializer(); + ~MultipleBSplineDeformableTransformInitializer(){}; + + private: + MultipleBSplineDeformableTransformInitializer(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + bool m_Verbose; + TransformPointer m_Transform; + ImagePointer m_Image; + + }; //class MultipleBSplineDeformableTransformInitializer + + +} // namespace clitk + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkMultipleBSplineDeformableTransformInitializer.txx" +#endif + +#endif /* __clitkMultipleBSplineDeformableTransformInitializer_h */ diff --git a/registration/clitkMultipleBSplineDeformableTransformInitializer.txx b/registration/clitkMultipleBSplineDeformableTransformInitializer.txx new file mode 100644 index 0000000..1ee4be9 --- /dev/null +++ b/registration/clitkMultipleBSplineDeformableTransformInitializer.txx @@ -0,0 +1,260 @@ +/*========================================================================= + 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://oncora1.lyon.fnclcc.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 __clitkMultipleBSplineDeformableTransformInitializer_txx +#define __clitkMultipleBSplineDeformableTransformInitializer_txx +#include "clitkMultipleBSplineDeformableTransformInitializer.h" + +namespace clitk +{ + + template < class TTransform, class TImage > + MultipleBSplineDeformableTransformInitializer + ::MultipleBSplineDeformableTransformInitializer() + { + this->m_NumberOfControlPointsInsideTheImage.Fill( 5 ); + this->m_ControlPointSpacing.Fill(64.); + this->m_ChosenSpacing.Fill(64.); + m_NumberOfControlPointsIsGiven=false; + m_ControlPointSpacingIsGiven=false; + m_TransformRegionIsGiven=false; + m_SamplingFactorIsGiven=false; + m_SplineOrders.Fill(3); + m_Parameters=NULL; + } + + template < class TTransform, class TImage > + void + MultipleBSplineDeformableTransformInitializer + ::InitializeTransform() + { + // Sanity check: + // The image is required for the region and spacing + if( ! this->m_Image ) + { + itkExceptionMacro( "Reference Image has not been set" ); + return; + } + + // A pointer to the transform is required + if( ! this->m_Transform ) + { + itkExceptionMacro( "Transform has not been set" ); + return; + } + + // If the image come from a filter, then update that filter. + if( this->m_Image->GetSource() ) + { + this->m_Image->GetSource()->Update(); + } + + + //-------------------------------------- + // Spacing & Size on image + //-------------------------------------- + SpacingType fixedImageSpacing = m_Image->GetSpacing(); + SizeType fixedImageSize=m_Image->GetLargestPossibleRegion().GetSize(); + typename RegionType::SizeType gridBorderSize; + typename RegionType::SizeType totalGridSize; + + // Only spacing given: adjust if necessary + if (m_ControlPointSpacingIsGiven && !m_NumberOfControlPointsIsGiven) + { + for(unsigned int r=0; rTransformIndexToPhysicalPoint(m_Image->GetLargestPossibleRegion().GetIndex(),fixedImageOrigin); + typename ImageType::DirectionType gridDirection = m_Image->GetDirection(); + SizeType orderShift; + + // Shift depends on order + for(unsigned int r=0; rm_Transform->SetSplineOrders(m_SplineOrders); + this->m_Transform->SetGridRegion( gridRegion ); + this->m_Transform->SetGridOrigin( gridOrigin ); + this->m_Transform->SetGridSpacing( m_ControlPointSpacing ); + this->m_Transform->SetGridDirection( gridDirection ); + this->m_Transform->SetLUTSamplingFactors(m_SamplingFactors); + + } + + template < class TTransform, class TImage > + void + MultipleBSplineDeformableTransformInitializer + ::SetInitialParameters( const std::string& s, ParametersType& params) + { + //--------------------------------------- + // Read Initial parameters + //--------------------------------------- + typedef itk::ImageFileReader CoefficientReaderType; + typename CoefficientReaderType::Pointer coeffReader = CoefficientReaderType::New(); + std::vector coefficientImages; + unsigned nLabels = m_Transform->GetnLabels(); + + int dotpos = s.length() - 1; + while (dotpos >= 0 && s[dotpos] != '.') + dotpos--; + + for (unsigned i = 0; i < nLabels; ++i) + { + std::ostringstream osfname; + osfname << s.substr(0, dotpos) << '_' << i << s.substr(dotpos); + coeffReader->SetFileName(osfname.str()); + if (m_Verbose) + std::cout << "Reading initial coefficients from file: " << osfname << "..." << std::endl; + coeffReader->Update(); + coefficientImages[i] = coeffReader->GetOutput(); + } + this->SetInitialParameters(coefficientImages, params); + } + + template < class TTransform, class TImage > + void + MultipleBSplineDeformableTransformInitializer + ::SetInitialParameters(std::vector coefficientImages, ParametersType& params) + { + // Keep a reference + m_Parameters=¶ms; + + // Resample + typedef ResampleBSplineDeformableTransformImageFilter ResamplerType; + typename ResamplerType::Pointer resampler; + unsigned int i = 0; + for (unsigned l = 0; l < coefficientImages.size(); ++l) + { + typename ResamplerType::Pointer resampler = ResamplerType::New(); + resampler->SetInput(coefficientImages[l]); + resampler->SetOutputParametersFromImage(m_Transform->GetCoefficientImages()[l]); + if (m_Verbose) + std::cout << "Resampling initial coefficients..." << std::endl; + resampler->Update(); + + // Copy parameters into the existing array, I know its crappy + typedef itk::ImageRegionConstIterator Iterator; + Iterator it(resampler->GetOutput(), resampler->GetOutput()->GetLargestPossibleRegion()); + it.GoToBegin(); + while (!it.IsAtEnd()) + { + for (unsigned int j = 0; j < InputDimension; ++j) + params[i + j] = it.Get()[j]; + ++it; + i += InputDimension; + } + } + + // JV pass the reference to the array !! + this->m_Transform->SetParameters(params); + } + + + +} // namespace itk + +#endif -- 2.47.1