able to optimize them at the same time.
--- /dev/null
+/*=========================================================================
+ 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 <clitkBSplineDeformableTransform.h>
+
+#include <itkNearestNeighborInterpolateImageFunction.h>
+
+#include <vector>
+
+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<TCoordRep, NInputDimensions, NOutputDimensions> Self;
+ typedef itk::Transform<TCoordRep, NInputDimensions, NOutputDimensions> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> 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<TCoordRep,
+ InputDimension> InputVectorType;
+ typedef itk::Vector<TCoordRep,
+ OutputDimension> OutputVectorType;
+
+ /** Standard covariant vector type for this class. */
+ typedef itk::CovariantVector<TCoordRep,
+ InputDimension> InputCovariantVectorType;
+ typedef itk::CovariantVector<TCoordRep,
+ OutputDimension> OutputCovariantVectorType;
+
+ /** Standard vnl_vector type for this class. */
+ typedef vnl_vector_fixed<TCoordRep,
+ InputDimension> InputVnlVectorType;
+ typedef vnl_vector_fixed<TCoordRep,
+ OutputDimension> OutputVnlVectorType;
+
+ /** Standard coordinate point type for this class. */
+ typedef itk::Point<TCoordRep,
+ InputDimension> InputPointType;
+ typedef itk::Point<TCoordRep,
+ OutputDimension> OutputPointType;
+
+ //JV Parameters as images with OutputDimension number of components per Pixel
+ typedef typename ParametersType::ValueType ParametersValueType;
+ typedef typename itk::Vector<ParametersValueType, OutputDimension> PixelType;
+ typedef itk::Image<PixelType, InputDimension> CoefficientImageType;
+ typedef typename CoefficientImageType::Pointer CoefficientImagePointer;
+
+ /** Typedefs for specifying the extend to the grid. */
+ typedef itk::ImageRegion<InputDimension> 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<TCoordRep, InputDimension> ContinuousIndexType;
+
+ /** Typedef of the bulk transform. */
+ typedef itk::Transform<ScalarType, InputDimension, OutputDimension> BulkTransformType;
+ typedef BulkTransformType* BulkTransformPointer;
+
+ typedef itk::Image< unsigned char, InputDimension> ImageLabelType;
+ typedef ImageLabelType* ImageLabelPointer;
+ //typedef typename ImageLabelType::Pointer ImageLabelPointer;
+
+ typedef itk::NearestNeighborInterpolateImageFunction<ImageLabelType, TCoordRep> 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<CoefficientImagePointer>& GetCoefficientImages() const;
+ virtual void SetCoefficientImages(std::vector<CoefficientImagePointer>& 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<TCoordRep,InputDimension, OutputDimension > 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<JacobianValueType,OutputDimension> JacobianPixelType;
+ typedef itk::Image<JacobianPixelType, OutputDimension> JacobianImageType;
+ typename JacobianImageType::Pointer m_JacobianImage[OutputDimension];
+ typedef itk::ImageRegionIterator<JacobianImageType> 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<typename TransformType::Pointer> m_trans;
+ std::vector<ParametersType> m_parameters;
+ mutable std::vector<CoefficientImagePointer> 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
--- /dev/null
+/*=========================================================================
+ 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 <itkRelabelComponentImageFilter.h>
+
+namespace clitk
+{
+ // Constructor with default arguments
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::~MultipleBSplineDeformableTransform()
+ {
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetLabels(ImageLabelPointer labels)
+ {
+ GetFixedParameters(); // Update m_FixedParameters
+ if (labels)
+ {
+ typedef itk::RelabelComponentImageFilter<ImageLabelType, ImageLabelType> 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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetSplineOrder(const unsigned int & splineOrder)
+ {
+ LOOP_ON_LABELS(SetSplineOrder, splineOrder);
+ this->Modified();
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetSplineOrders(const SizeType & splineOrders)
+ {
+ LOOP_ON_LABELS(SetSplineOrders, splineOrders);
+ this->Modified();
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetLUTSamplingFactor( const int & samplingFactor)
+ {
+ LOOP_ON_LABELS(SetLUTSamplingFactor, samplingFactor);
+ this->Modified();
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetLUTSamplingFactors( const SizeType & samplingFactors)
+ {
+ LOOP_ON_LABELS(SetLUTSamplingFactors, samplingFactors);
+ this->Modified();
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetGridRegion( const RegionType & region )
+ {
+ LOOP_ON_LABELS(SetGridRegion, region);
+ this->Modified();
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetGridSpacing( const SpacingType & spacing )
+ {
+ LOOP_ON_LABELS(SetGridSpacing, spacing);
+ this->Modified();
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetGridDirection( const DirectionType & direction )
+ {
+ LOOP_ON_LABELS(SetGridDirection, direction);
+ this->Modified();
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetGridOrigin( const OriginType& origin )
+ {
+ LOOP_ON_LABELS(SetGridOrigin, origin);
+ this->Modified();
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetIdentity()
+ {
+ LOOP_ON_LABELS(SetIdentity, );
+ if ( m_InputParametersPointer )
+ {
+ ParametersType * parameters =
+ const_cast<ParametersType *>( m_InputParametersPointer );
+ parameters->Fill( 0.0 );
+ }
+ this->Modified();
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline const std::vector<typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::CoefficientImagePointer>&
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetCoefficientImages(std::vector<CoefficientImagePointer>& 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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<double *>(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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::SetBulkTransform(BulkTransformPointer b)
+ {
+ m_BulkTransform = b;
+ LOOP_ON_LABELS(SetBulkTransform, b);
+ }
+
+#undef LOOP_ON_LABELS
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline unsigned int
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::GetNumberOfParameters(void) const
+ {
+ unsigned int sum = 0;
+ for (unsigned i = 0; i < m_nLabels; ++i)
+ sum += m_trans[i]->GetNumberOfParameters();
+ return sum;
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline unsigned int
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::GetNumberOfParametersPerDimension(void) const
+ {
+ unsigned int sum = 0;
+ for (unsigned i = 0; i < m_nLabels; ++i)
+ sum += m_trans[i]->GetNumberOfParametersPerDimension();
+ return sum;
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::ParametersType &
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::ParametersType &
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline bool
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::OutputPointType
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::OutputPointType
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline const typename MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>::JacobianType &
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>
+ ::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<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>
+ ::TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const
+ {
+ m_trans[0]->TransformPointToContinuousIndex(point, index);
+ }
+
+ template<class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
+ inline void
+ MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions,NOutputDimensions>::InitJacobian()
+ {
+ unsigned numberOfParameters = this->GetNumberOfParameters();
+ this->m_Jacobian.set_size(OutputDimension, numberOfParameters);
+ JacobianPixelType * jacobianDataPointer = reinterpret_cast<JacobianPixelType *>(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
--- /dev/null
+/*=========================================================================
+ 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 <iostream>
+
+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<Self> Pointer;
+ typedef itk::SmartPointer<const Self> 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;i<InputDimension;i++)
+ s[i]=n[i];
+ this->SetNumberOfControlPointsInsideTheImage( 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<InputDimension;i++)
+ s[i]=n[i];
+ this-> SetSamplingFactors( s );
+ }
+ void SetSamplingFactors( unsigned int n )
+ {
+ SizeType s;
+ s.Fill( n );
+ this-> SetSamplingFactors( s );
+ }
+ virtual void InitializeTransform();
+ void SetInitialParameters(std::vector<typename CoefficientImageType::Pointer> coefficientImage, ParametersType& params);
+ void SetInitialParameters(const std::string & s, ParametersType& params);
+
+ // For easy acces, declared public
+ std::vector<SizeType> m_NumberOfControlPointsInsideTheImageArray;
+ std::vector<SizeType> m_SamplingFactorsArray;
+ std::vector<SpacingType> 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 */
--- /dev/null
+/*=========================================================================
+ 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<TTransform, TImage >
+ ::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<TTransform, TImage >
+ ::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; r<InputDimension; r++)
+ {
+ // JV
+ m_ChosenSpacing[r]= m_ControlPointSpacing[r];
+ m_ControlPointSpacing[r]= ( round(m_ChosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
+ m_NumberOfControlPointsInsideTheImage[r] = ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] );
+ if ( ( ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
+ == ( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
+ {
+ m_NumberOfControlPointsInsideTheImage[r] +=1;
+ DD("Adding control point");
+ }
+ }
+ }
+
+ // Only number of CP given: adjust if necessary
+ // JV -1 ?
+ else if (!m_ControlPointSpacingIsGiven && m_NumberOfControlPointsIsGiven)
+ {
+ for(unsigned int r=0; r<InputDimension; r++)
+ {
+ m_ChosenSpacing[r]=fixedImageSpacing[r]*( (double)(fixedImageSize[r]) /
+ (double)(m_NumberOfControlPointsInsideTheImage[r]) );
+ m_ControlPointSpacing[r]= fixedImageSpacing[r]* ceil( (double)(fixedImageSize[r] - 1) /
+ (double)(m_NumberOfControlPointsInsideTheImage[r] - 1) );
+ }
+ }
+
+ // Both or none of Spacing and number of CP given: no adjustment adjust, just warnings
+ else
+ {
+ for(unsigned int r=0; r<InputDimension; r++)
+ {
+ m_ChosenSpacing[r]= m_ControlPointSpacing[r];
+ if (m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]< fixedImageSize[r]*fixedImageSpacing[r])
+ {
+ std::cout<<"WARNING: Specified control point region ("<<m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]
+ <<"mm) does not cover the fixed image region ("<< fixedImageSize[r]*fixedImageSpacing[r]
+ <<"mm) for dimension "<<r<<"!" <<std::endl
+ <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
+ }
+ if ( fmod(m_ControlPointSpacing[r], fixedImageSpacing[r]) )
+ {
+ std::cout<<"WARNING: Specified control point spacing for dimension "<<r
+ <<" does not allow exact representation of BLUT FFD!"<<std::endl
+ <<"Spacing ratio is non-integer: "<<m_ControlPointSpacing[r]/ fixedImageSpacing[r]<<std::endl
+ <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
+ }
+ }
+ }
+
+ if (m_Verbose) std::cout<<"The chosen control point spacing "<<m_ChosenSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<m_ControlPointSpacing<<"..."<<std::endl;
+ if (m_Verbose) std::cout<<"The number of (internal) control points is "<<m_NumberOfControlPointsInsideTheImage<<"..."<<std::endl;
+
+
+ //--------------------------------------
+ // Border size
+ //--------------------------------------
+ for(unsigned int r=0; r<InputDimension; r++) gridBorderSize[r]=m_SplineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
+ totalGridSize = m_NumberOfControlPointsInsideTheImage + gridBorderSize;
+ if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
+ RegionType gridRegion;
+ gridRegion.SetSize(totalGridSize);
+
+
+ //--------------------------------------
+ // Origin
+ //--------------------------------------
+ typedef typename TransformType::OriginType OriginType;
+ OriginType fixedImageOrigin, gridOrigin;
+
+ // In case of non-zero index
+ m_Image->TransformIndexToPhysicalPoint(m_Image->GetLargestPossibleRegion().GetIndex(),fixedImageOrigin);
+ typename ImageType::DirectionType gridDirection = m_Image->GetDirection();
+ SizeType orderShift;
+
+ // Shift depends on order
+ for(unsigned int r=0; r<InputDimension; r++)
+ {
+ orderShift[r] = m_SplineOrders[r] / 2;
+ gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
+ }
+ if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
+
+
+ //--------------------------------------
+ // LUT sampling factors
+ //--------------------------------------
+ for (unsigned int i=0; i< InputDimension; i++)
+ {
+ if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
+ << m_ControlPointSpacing[i]/ fixedImageSpacing[i]<<"..."<<std::endl;
+ if ( !m_SamplingFactorIsGiven) m_SamplingFactors[i]= (int) ( m_ControlPointSpacing[i]/ fixedImageSpacing[i]);
+ if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<m_SamplingFactors[i]<<"..."<<std::endl;
+ }
+
+
+ //--------------------------------------
+ // Set
+ //--------------------------------------
+ this->m_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<TTransform, TImage >
+ ::SetInitialParameters( const std::string& s, ParametersType& params)
+ {
+ //---------------------------------------
+ // Read Initial parameters
+ //---------------------------------------
+ typedef itk::ImageFileReader<CoefficientImageType> CoefficientReaderType;
+ typename CoefficientReaderType::Pointer coeffReader = CoefficientReaderType::New();
+ std::vector<typename CoefficientImageType::Pointer> 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<TTransform, TImage >
+ ::SetInitialParameters(std::vector<typename CoefficientImageType::Pointer> coefficientImages, ParametersType& params)
+ {
+ // Keep a reference
+ m_Parameters=¶ms;
+
+ // Resample
+ typedef ResampleBSplineDeformableTransformImageFilter<CoefficientImageType,CoefficientImageType> 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<CoefficientImageType> 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
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff