diff --git a/Utilities/InsightJournal/ijBSplineKernelFunction.h b/Utilities/InsightJournal/ijBSplineKernelFunction.h
deleted file mode 100644
index e2a2e233b5b71fbde6b18dea22a0081870eac2c5..0000000000000000000000000000000000000000
--- a/Utilities/InsightJournal/ijBSplineKernelFunction.h
+++ /dev/null
@@ -1,175 +0,0 @@
-/*=========================================================================
-
- Program: Insight Segmentation & Registration Toolkit
- Module: $RCSfile: ijBSplineKernelFunction.h,v $
- Language: C++
- Date: $Date: $
- Version: $Revision: $
-
- Copyright (c) Insight Software Consortium. All rights reserved.
- See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-#ifndef _ijBSplineKernelFunction_h
-#define _ijBSplineKernelFunction_h
-
-#include "itkKernelFunction.h"
-#include "vnl/vnl_math.h"
-#include "vnl/vnl_real_polynomial.h"
-#include "vnl/vnl_matrix.h"
-
-namespace ij
-{
-
-/** \class BSplineKernelFunction
- * \brief BSpline kernel used for density estimation and nonparameteric
- * regression.
- *
- * This class enscapsulates BSpline kernel for
- * density estimation or nonparameteric regression.
- * See documentation for KernelFunction for more details.
- *
- * This class is templated over the spline order to cohere with
- * the previous incarnation of this class. One can change the
- * order during an instantiation's existence. Note that
- * other authors have defined the B-spline order as being the
- * degree of spline + 1. In the ITK context (e.g. in this
- * class), the spline order is equivalent to the degree of
- * the spline.
- *
- * \sa KernelFunction
- *
- * \ingroup Functions
- */
-template <unsigned int VSplineOrder = 3>
-class ITK_EXPORT BSplineKernelFunction
-: public itk::KernelFunction
-{
-public:
- /** Standard class typedefs. */
- typedef BSplineKernelFunction Self;
- typedef itk::KernelFunction Superclass;
- typedef itk::SmartPointer<Self> Pointer;
-
- /** Method for creation through the object factory. */
- itkNewMacro(Self);
-
- /** Run-time type information (and related methods). */
- itkTypeMacro( BSplineKernelFunction, KernelFunction );
-
- typedef double RealType;
- typedef vnl_vector<RealType> VectorType;
- typedef vnl_real_polynomial PolynomialType;
- typedef vnl_matrix<RealType> MatrixType;
-
- /** Get/Sets the Spline Order */
- void SetSplineOrder( unsigned int );
- itkGetMacro( SplineOrder, unsigned int );
-
- /** Evaluate the function. */
- inline RealType Evaluate( const RealType & u ) const
- {
- RealType absValue = vnl_math_abs( u );
- int which;
- if ( this->m_SplineOrder % 2 == 0 )
- {
- which = static_cast< int>( absValue+0.5 );
- }
- else
- {
- which = static_cast< int>( absValue );
- }
- if ( which < static_cast<int>(this->m_BSplineShapeFunctions.rows()) )
- {
- return PolynomialType(
- this->m_BSplineShapeFunctions.get_row( which ) ).evaluate( absValue );
- }
- else
- {
- return itk::NumericTraits<RealType>::Zero;
- }
- }
-
- /** Evaluate the derivative. */
- inline RealType EvaluateDerivative( const double & u ) const
- {
- RealType absValue = vnl_math_abs( u );
- int which;
- if ( this->m_SplineOrder % 2 == 0 )
- {
- which = static_cast<unsigned int>( absValue+0.5 );
- }
- else
- {
- which = static_cast<unsigned int>( absValue );
- }
- if ( which < this->m_BSplineShapeFunctions.rows() )
- {
- RealType der = PolynomialType(
- this->m_BSplineShapeFunctions.get_row( which ) ).devaluate( absValue );
- if ( u < itk::NumericTraits<RealType>::Zero )
- {
- return -der;
- }
- else
- {
- return der;
- }
- }
- else
- {
- return itk::NumericTraits<RealType>::Zero;
- }
- }
-
- /**
- * For a specific order, return the (this->m_SplineOrder+1) pieces of
- * the single basis function centered at zero.
- */
- MatrixType GetShapeFunctions();
-
- /**
- * For a specific order, generate and return the (this->m_SplineOrder+1)
- * pieces of the different basis functions in the [0, 1] interval.
- */
- MatrixType GetShapeFunctionsInZeroToOneInterval();
-
-protected:
- BSplineKernelFunction();
- ~BSplineKernelFunction();
- void PrintSelf( std::ostream& os, itk::Indent indent ) const;
-
-private:
- BSplineKernelFunction( const Self& ); //purposely not implemented
- void operator=( const Self& ); //purposely not implemented
-
- /**
- * For a specific order, generate the (this->m_SplineOrder+1) pieces of
- * the single basis function centered at zero.
- */
- void GenerateBSplineShapeFunctions( unsigned int );
-
- /**
- * Use the CoxDeBoor recursion relation to generate the piecewise
- * polynomials which compose the basis function.
- * See, for example, L. Piegl, L. Tiller, "The NURBS Book,"
- * Springer 1997, p. 50.
- */
- PolynomialType CoxDeBoor( unsigned short, VectorType, unsigned int, unsigned int );
-
- MatrixType m_BSplineShapeFunctions;
- unsigned int m_SplineOrder;
-
-};
-
-} // end namespace itk
-
-#ifndef OTB_MANUAL_INSTANTIATION
-#include "ijBSplineKernelFunction.txx"
-#endif
-
-#endif
diff --git a/Utilities/InsightJournal/ijBSplineKernelFunction.txx b/Utilities/InsightJournal/ijBSplineKernelFunction.txx
deleted file mode 100644
index b9c9591b941c7a8d82a9d6ffa3afef30ab9d75f5..0000000000000000000000000000000000000000
--- a/Utilities/InsightJournal/ijBSplineKernelFunction.txx
+++ /dev/null
@@ -1,185 +0,0 @@
-/*=========================================================================
-
- Program: Insight Segmentation & Registration Toolkit
- Module: $RCSfile: itkBSplineKernelFunction.txx,v $
- Language: C++
- Date: $Date: $
- Version: $Revision: $
-
- Copyright (c) Insight Software Consortium. All rights reserved.
- See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-#ifndef _ijBSplineKernelFunction_txx
-#define _ijBSplineKernelFunction_txx
-
-#include "ijBSplineKernelFunction.h"
-
-namespace ij
-{
-
-template <unsigned int VSplineOrder>
-BSplineKernelFunction<VSplineOrder>
-::BSplineKernelFunction()
-{
- this->m_SplineOrder = VSplineOrder;
- this->GenerateBSplineShapeFunctions( this->m_SplineOrder+1 );
-}
-
-template <unsigned int VSplineOrder>
-BSplineKernelFunction<VSplineOrder>
-::~BSplineKernelFunction()
-{
-}
-
-template <unsigned int VSplineOrder>
-void
-BSplineKernelFunction<VSplineOrder>
-::SetSplineOrder( unsigned int order )
-{
- if ( order != this->m_SplineOrder )
- {
- this->m_SplineOrder = order;
- this->GenerateBSplineShapeFunctions( this->m_SplineOrder+1 );
- this->Modified();
- }
-}
-
-template <unsigned int VSplineOrder>
-void
-BSplineKernelFunction<VSplineOrder>
-::GenerateBSplineShapeFunctions( unsigned int order )
-{
- unsigned int NumberOfPieces = static_cast<unsigned int>( 0.5*( order+1 ) );
- this->m_BSplineShapeFunctions.set_size( NumberOfPieces, order );
-
- VectorType knots( order+1 );
- for ( unsigned int i = 0; i < knots.size(); i++)
- {
- knots[i] = -0.5*static_cast<RealType>( order ) + static_cast<RealType>( i );
- }
-
- for ( unsigned int i = 0; i < NumberOfPieces; i++ )
- {
- PolynomialType poly = this->CoxDeBoor(order, knots,
- 0, static_cast<unsigned int>( 0.5*( order ) ) + i );
- this->m_BSplineShapeFunctions.set_row( i, poly.coefficients() );
- }
-}
-
-template <unsigned int VSplineOrder>
-typename BSplineKernelFunction<VSplineOrder>::PolynomialType
-BSplineKernelFunction<VSplineOrder>
-::CoxDeBoor( unsigned short order, VectorType knots,
- unsigned int whichBasisFunction, unsigned int whichPiece )
-{
- VectorType tmp(2);
- PolynomialType poly1(0.0), poly2(0.0);
- RealType den;
- unsigned short p = order-1;
- unsigned short i = whichBasisFunction;
-
- if ( p == 0 && whichBasisFunction == whichPiece )
- {
- PolynomialType poly(1.0);
- return poly;
- }
-
- // Term 1
- den = knots(i+p)-knots(i);
- if ( den == itk::NumericTraits<RealType>::Zero )
- {
- PolynomialType poly(0.0);
- poly1 = poly;
- }
- else
- {
- tmp(0) = 1.0;
- tmp(1) = -knots(i);
- tmp /= den;
- poly1 = PolynomialType(tmp) * this->CoxDeBoor( order-1, knots, i, whichPiece );
- }
-
- // Term 2
- den = knots(i+p+1)-knots(i+1);
- if ( den == itk::NumericTraits<RealType>::Zero )
- {
- PolynomialType poly(0.0);
- poly2 = poly;
- }
- else
- {
- tmp(0) = -1.0;
- tmp(1) = knots(i+p+1);
- tmp /= den;
- poly2 = PolynomialType(tmp) * this->CoxDeBoor( order-1, knots, i+1, whichPiece );
- }
- return ( poly1 + poly2 );
-}
-
-template <unsigned int VSplineOrder>
-typename BSplineKernelFunction<VSplineOrder>::MatrixType
-BSplineKernelFunction<VSplineOrder>
-::GetShapeFunctionsInZeroToOneInterval()
-{
- int order = this->m_SplineOrder+1;
- unsigned int NumberOfPieces = static_cast<unsigned int>( order );
- MatrixType ShapeFunctions( NumberOfPieces, order );
-
- VectorType knots( 2*order );
- for ( unsigned int i = 0; i < knots.size(); i++ )
- {
- knots[i] = -static_cast<RealType>( this->m_SplineOrder )
- + static_cast<RealType>( i );
- }
-
- for (unsigned int i = 0; i < NumberOfPieces; i++ )
- {
- PolynomialType poly = this->CoxDeBoor( order, knots, i, order-1 );
- ShapeFunctions.set_row( i, poly.coefficients() );
- }
- return ShapeFunctions;
-}
-
-template <unsigned int VSplineOrder>
-void
-BSplineKernelFunction<VSplineOrder>
-::PrintSelf( std::ostream& os, itk::Indent indent ) const
-{
- Superclass::PrintSelf( os, indent );
- os << indent << "Spline Order: " << this->m_SplineOrder << std::endl;
- os << indent << "Piecewise Polynomial Pieces: " << std::endl;
- RealType a(0.), b(0.);
- for ( unsigned int i = 0; i < this->m_BSplineShapeFunctions.rows(); i++ )
- {
- os << indent << indent;
- PolynomialType( this->m_BSplineShapeFunctions.get_row( i ) ).print( os );
- if ( i == 0 )
- {
- a = 0.0;
- if ( this->m_SplineOrder % 2 == 0 )
- {
- b = 0.5;
- }
- else
- {
- b = 1.0;
- }
- }
- else
- {
- a = b;
- b += 1.0;
- }
- os << ", X \\in [" << a << ", " << b << "]" << std::endl;
- }
-}
-
-
-} // namespace ij
-
-#endif
diff --git a/Utilities/InsightJournal/ijBSplineScatteredDataPointSetToImageFilter.h b/Utilities/InsightJournal/ijBSplineScatteredDataPointSetToImageFilter.h
deleted file mode 100644
index 68432151ac7d9e9db06def65c94a420adcd1c59b..0000000000000000000000000000000000000000
--- a/Utilities/InsightJournal/ijBSplineScatteredDataPointSetToImageFilter.h
+++ /dev/null
@@ -1,236 +0,0 @@
-/*=========================================================================
-
- Program: Insight Segmentation & Registration Toolkit
- Module: $RCSfile: itkBSplineScatteredDataPointSetToImageFilter.h,v $
- Language: C++
- Date: $Date: $
- Version: $Revision: $
-
- Copyright (c) Insight Software Consortium. All rights reserved.
- See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-#ifndef __ijBSplineScatteredDataPointSetToImageFilter_h
-#define __ijBSplineScatteredDataPointSetToImageFilter_h
-
-#include "ijPointSetToImageFilter.h"
-#include "ijBSplineKernelFunction.h"
-#include "itkFixedArray.h"
-#include "itkVariableSizeMatrix.h"
-#include "itkVector.h"
-#include "itkVectorContainer.h"
-
-#include "vnl/vnl_matrix.h"
-
-namespace ij
-{
-
-/** \class BSplineScatteredDataPointSetToImageFilter.h
- * \brief Image filter which provides a B-spline output approximation.
- *
- * Given an n-D image with scattered data, this filter finds
- * a fast approximation to that irregulary spaced data using uniform
- * B-splines. The traditional method of inverting the observation
- * matrix to find a least-squares fit is made obsolete. Therefore,
- * memory issues are not a concern and inverting large matrices are
- * unnecessary. The reference below gives the algorithm for 2-D images
- * and cubic splines. This class generalizes that work to encompass n-D
- * images and any *feasible* B-spline order.
- *
- * In addition to specifying the input point set, one must specify the number
- * of control points. If one wishes to use the multilevel component of
- * this algorithm, one must also specify the number of levels in the
- * hieararchy. If this is desired, the number of control points becomes
- * the number of control points for the coarsest level. The algorithm
- * then increases the number of control points at each level so that
- * the B-spline n-D grid is refined to twice the previous level. The
- * scattered data is specified by the pixel values. Pixels which
- * are not to be included in the calculation of the B-spline grid must
- * have a value equal to that specified by the variable m_BackgroundValue.
- *
- * Note that the specified number of control points must be > m_SplineOrder.
- *
- * \par REFERENCE
- * S. Lee, G. Wolberg, and S. Y. Shin, "Scattered Data Interpolation
- * with Multilevel B-Splines", IEEE Transactions on Visualization and
- * Computer Graphics, 3(3):228-244, 1997.
- *
- * N.J. Tustison and J.C. Gee, "Generalized n-D C^k Scattered Data Approximation
- * with COnfidence Values", Proceedings of the MIAR conference, August 2006.
- */
-
-template <class TInputPointSet, class TOutputImage>
-class BSplineScatteredDataPointSetToImageFilter
-: public ij::PointSetToImageFilter<TInputPointSet, TOutputImage>
-{
-public:
- typedef BSplineScatteredDataPointSetToImageFilter Self;
- typedef PointSetToImageFilter<TInputPointSet, TOutputImage> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
-
- /** Method for creation through the object factory. */
- itkNewMacro(Self);
-
- /** Extract dimension from input and output image. */
- itkStaticConstMacro( ImageDimension, unsigned int,
- TOutputImage::ImageDimension );
-
- typedef TOutputImage ImageType;
- typedef TInputPointSet PointSetType;
-
- /** Image typedef support. */
- typedef typename ImageType::PixelType PixelType;
- typedef typename ImageType::RegionType RegionType;
- typedef typename ImageType::SizeType SizeType;
- typedef typename ImageType::IndexType IndexType;
- typedef typename ImageType::PointType ContinuousIndexType;
-
- /** PointSet typedef support. */
- typedef typename PointSetType::PointType PointType;
- typedef typename PointSetType::PixelType PointDataType;
- typedef typename PointSetType::PointDataContainer PointDataContainerType;
-
- /** Other typedef */
- typedef double RealType;
- typedef itk::VectorContainer<unsigned, RealType> WeightsContainerType;
- typedef itk::FixedArray<unsigned,
- itkGetStaticConstMacro( ImageDimension )> ArrayType;
- typedef itk::VariableSizeMatrix<RealType> GradientType;
- typedef itk::Image<RealType,
- itkGetStaticConstMacro( ImageDimension )> RealImageType;
- typedef itk::Image<PointDataType,
- itkGetStaticConstMacro( ImageDimension )> PointDataImageType;
-
- /** Interpolation kernel type (default spline order = 3) */
- typedef ij::BSplineKernelFunction<3> KernelType;
-
- /** Helper functions */
-
- void SetNumberOfLevels( unsigned int );
- void SetNumberOfLevels( ArrayType );
- itkGetConstReferenceMacro( NumberOfLevels, ArrayType );
-
- void SetSplineOrder( unsigned int );
- void SetSplineOrder( ArrayType );
- itkGetConstReferenceMacro( SplineOrder, ArrayType );
-
- itkSetMacro( NumberOfControlPoints, ArrayType );
- itkGetConstReferenceMacro( NumberOfControlPoints, ArrayType );
- itkGetConstReferenceMacro( CurrentNumberOfControlPoints, ArrayType );
-
- itkSetMacro( CloseDimension, ArrayType );
- itkGetConstReferenceMacro( CloseDimension, ArrayType );
-
- itkSetMacro( GenerateOutputImage, bool );
- itkGetConstReferenceMacro( GenerateOutputImage, bool );
- itkBooleanMacro( GenerateOutputImage );
-
- void SetPointWeights( typename WeightsContainerType::Pointer weights );
-
- /**
- * Get the control point lattice.
- */
- itkGetConstMacro( PhiLattice, typename PointDataImageType::Pointer );
-
- /**
- * Evaluate the resulting B-spline object at a specified
- * point or index within the image domain.
- */
- void EvaluateAtPoint( PointType, PointDataType & );
- void EvaluateAtIndex( IndexType, PointDataType & );
- void EvaluateAtContinuousIndex( ContinuousIndexType, PointDataType & );
-
- /**
- * Evaluate the resulting B-spline object at a specified
- * parameteric point. Note that the parameterization over
- * each dimension of the B-spline object is [0, 1).
- */
- void Evaluate( PointType, PointDataType & );
-
- /**
- * Evaluate the gradient of the resulting B-spline object at a specified
- * point or index within the image domain.
- */
- void EvaluateGradientAtPoint( PointType, GradientType & );
- void EvaluateGradientAtIndex( IndexType, GradientType & );
- void EvaluateGradientAtContinuousIndex( ContinuousIndexType, GradientType & );
-
- /**
- * Evaluate the gradient of the resulting B-spline object
- * at a specified parameteric point. Note that the
- * parameterization over each dimension of the B-spline
- * object is [0, 1).
- */
- void EvaluateGradient( PointType, GradientType & );
-
-protected:
- BSplineScatteredDataPointSetToImageFilter();
- virtual ~BSplineScatteredDataPointSetToImageFilter();
- void PrintSelf( std::ostream& os, itk::Indent indent ) const;
-
- void GenerateData();
-
-private:
- BSplineScatteredDataPointSetToImageFilter(const Self&); //purposely not implemented
- void operator=(const Self&); //purposely not implemented
-
- void GenerateControlLattice();
- void RefineControlLattice();
- void UpdatePointSet();
- void GenerateOutputImage();
- void GenerateOutputImageFast();
- void CollapsePhiLattice( PointDataImageType *,
- PointDataImageType *,
- RealType, unsigned int );
-
-
- bool m_DoMultilevel;
- bool m_GenerateOutputImage;
- bool m_UsePointWeights;
- unsigned int m_MaximumNumberOfLevels;
- unsigned int m_CurrentLevel;
- ArrayType m_NumberOfControlPoints;
- ArrayType m_CurrentNumberOfControlPoints;
- ArrayType m_CloseDimension;
- ArrayType m_SplineOrder;
- ArrayType m_NumberOfLevels;
- typename WeightsContainerType::Pointer m_PointWeights;
- typename KernelType::Pointer m_Kernel[ImageDimension];
- typename KernelType::MatrixType m_RefinedLatticeCoefficients[ImageDimension];
- typename PointDataImageType::Pointer m_PhiLattice;
- typename PointDataImageType::Pointer m_PsiLattice;
- typename PointDataContainerType::Pointer m_InputPointData;
- typename PointDataContainerType::Pointer m_OutputPointData;
-
- inline typename RealImageType::IndexType
- IndexToSubscript(unsigned int index, typename RealImageType::SizeType size)
- {
- typename RealImageType::IndexType k;
- k[0] = 1;
-
- for ( unsigned int i = 1; i < ImageDimension; i++ )
- {
- k[i] = size[ImageDimension-i-1]*k[i-1];
- }
- typename RealImageType::IndexType sub;
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- sub[ImageDimension-i-1] = static_cast<unsigned int>( index/k[ImageDimension-i-1] );
- index %= k[ImageDimension-i-1];
- }
- return sub;
- }
-};
-} // end namespace ij
-
-#ifndef OTB_MANUAL_INSTANTIATION
-#include "ijBSplineScatteredDataPointSetToImageFilter.txx"
-#endif
-
-#endif
-
diff --git a/Utilities/InsightJournal/ijBSplineScatteredDataPointSetToImageFilter.txx b/Utilities/InsightJournal/ijBSplineScatteredDataPointSetToImageFilter.txx
deleted file mode 100644
index 7cd4a7279915ef9c700382977973212b80330fc1..0000000000000000000000000000000000000000
--- a/Utilities/InsightJournal/ijBSplineScatteredDataPointSetToImageFilter.txx
+++ /dev/null
@@ -1,970 +0,0 @@
-#ifndef __ijBSplineScatteredDataPointSetToImageFilter_txx
-#define __ijBSplineScatteredDataPointSetToImageFilter_txx
-
-#include "ijBSplineScatteredDataPointSetToImageFilter.h"
-#include "itkImageRegionConstIteratorWithIndex.h"
-#include "itkImageRegionIterator.h"
-#include "itkImageRegionIteratorWithIndex.h"
-#include "itkImageDuplicator.h"
-#include "itkCastImageFilter.h"
-#include "itkNumericTraits.h"
-#include "itkTimeProbe.h"
-
-#include "vnl/vnl_math.h"
-#include "vnl/algo/vnl_matrix_inverse.h"
-#include "vnl/vnl_vector.h"
-namespace ij
-{
-
-template <class TInputPointSet, class TOutputImage>
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::BSplineScatteredDataPointSetToImageFilter()
-{
- this->m_SplineOrder.Fill( 3 );
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- this->m_NumberOfControlPoints[i] = ( this->m_SplineOrder[i]+1 );
- this->m_Kernel[i] = KernelType::New();
- this->m_Kernel[i]->SetSplineOrder( this->m_SplineOrder[i] );
- }
-
- this->m_CloseDimension.Fill( 0 );
- this->m_DoMultilevel = false;
- this->m_GenerateOutputImage = true;
- this->m_NumberOfLevels.Fill( 1 );
- this->m_MaximumNumberOfLevels = 1;
-
- this->m_PhiLattice = PointDataImageType::New();
- this->m_PsiLattice = PointDataImageType::New();
- this->m_InputPointData = PointDataContainerType::New();
- this->m_OutputPointData = PointDataContainerType::New();
-
- this->m_PointWeights = WeightsContainerType::New();
- this->m_UsePointWeights = false;
-}
-
-template <class TInputPointSet, class TOutputImage>
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::~BSplineScatteredDataPointSetToImageFilter()
-{
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::SetSplineOrder( unsigned int order )
-{
- this->m_SplineOrder.Fill( order );
- this->SetSplineOrder( this->m_SplineOrder );
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::SetSplineOrder( ArrayType order )
-{
- itkDebugMacro(<< "Setting m_SplineOrder to " << order);
-
- this->m_SplineOrder = order;
- for ( int i = 0; i < ImageDimension; i++ )
- {
- if ( this->m_SplineOrder[i] == 0 )
- {
- itkExceptionMacro( "The spline order in each dimension must be greater than 0" );
- }
-
- this->m_Kernel[i] = KernelType::New();
- this->m_Kernel[i]->SetSplineOrder( this->m_SplineOrder[i] );
-
- if ( this->m_DoMultilevel )
- {
- typename KernelType::MatrixType R;
- typename KernelType::MatrixType C;
- C = this->m_Kernel[i]->GetShapeFunctionsInZeroToOneInterval();
- R = C;
- for ( unsigned int j = 0; j < C.cols(); j++ )
- {
- RealType c = pow( 2.0, static_cast<RealType>( C.cols()-j-1 ) );
- for ( unsigned int k = 0; k < C.rows(); k++)
- {
- R(k, j) *= c;
- }
- }
- R = R.transpose();
- R.flipud();
- C = C.transpose();
- C.flipud();
- this->m_RefinedLatticeCoefficients[i]
- = ( vnl_svd<RealType>( R ).solve( C ) ).extract( 2, C.cols() );
- }
- }
- this->Modified();
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::SetNumberOfLevels( unsigned int levels )
-{
- this->m_NumberOfLevels.Fill( levels );
- this->SetNumberOfLevels( this->m_NumberOfLevels );
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::SetNumberOfLevels( ArrayType levels )
-{
- this->m_NumberOfLevels = levels;
- this->m_MaximumNumberOfLevels = 1;
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- if ( this->m_NumberOfLevels[i] == 0 )
- {
- itkExceptionMacro( "The number of levels in each dimension must be greater than 0" );
- }
- if ( this->m_NumberOfLevels[i] > this->m_MaximumNumberOfLevels )
- {
- this->m_MaximumNumberOfLevels = this->m_NumberOfLevels[i];
- }
- }
-
- itkDebugMacro( "Setting m_NumberOfLevels to " <<
- this->m_NumberOfLevels );
- itkDebugMacro( "Setting m_MaximumNumberOfLevels to " <<
- this->m_MaximumNumberOfLevels );
-
- if ( this->m_MaximumNumberOfLevels > 1 )
- {
- this->m_DoMultilevel = true;
- }
- else
- {
- this->m_DoMultilevel = false;
- }
- this->SetSplineOrder( this->m_SplineOrder );
- this->Modified();
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::SetPointWeights( typename WeightsContainerType::Pointer weights )
-{
- this->m_UsePointWeights = true;
- this->m_PointWeights = weights;
- this->Modified();
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::GenerateData()
-{
- if ( this->GetInput()->GetNumberOfPoints() == 0 )
- {
- itkExceptionMacro( "The number of points must be greater than 0." );
- }
- if ( this->m_UsePointWeights &&
- ( this->m_PointWeights->Size() != this->GetInput()->GetNumberOfPoints() ) )
- {
- itkExceptionMacro( "The number of weight points and input points must be equal." );
- }
-
-
- /**
- * Create the output image
- */
-
- for ( unsigned int i = 0; i < ImageDimension; i++)
- {
- if ( this->m_Size[i] == 0 )
- {
- itkExceptionMacro( "Size must be specified." );
- }
- }
- this->GetOutput()->SetOrigin( this->m_Origin );
- this->GetOutput()->SetSpacing( this->m_Spacing );
- this->GetOutput()->SetRegions( this->m_Size );
- this->GetOutput()->Allocate();
-
- this->m_InputPointData->Initialize();
- this->m_OutputPointData->Initialize();
-
- typename PointSetType::PointDataContainer::ConstIterator It;
- It = this->GetInput()->GetPointData()->Begin();
- while ( It != this->GetInput()->GetPointData()->End() )
- {
- if ( !this->m_UsePointWeights )
- {
- this->m_PointWeights->InsertElement( It.Index(), 1.0 );
- }
- this->m_InputPointData->InsertElement( It.Index(), It.Value() );
- this->m_OutputPointData->InsertElement( It.Index(), It.Value() );
- ++It;
- }
-
- this->m_CurrentLevel = 0;
- this->m_CurrentNumberOfControlPoints = this->m_NumberOfControlPoints;
- this->GenerateControlLattice();
- this->UpdatePointSet();
-
- itkDebugMacro( << "Current Level = " << this->m_CurrentLevel );
- itkDebugMacro( << " Number of control points = "
- << this->m_CurrentNumberOfControlPoints );
-
- if ( this->m_DoMultilevel )
- {
- this->m_PsiLattice->SetRegions( this->m_PhiLattice->GetLargestPossibleRegion() );
- this->m_PsiLattice->Allocate();
- PointDataType P( 0.0 );
- this->m_PsiLattice->FillBuffer( P );
- }
-
- for ( this->m_CurrentLevel = 1;
- this->m_CurrentLevel < this->m_MaximumNumberOfLevels; this->m_CurrentLevel++ )
- {
-
- itk::ImageRegionIterator<PointDataImageType> ItPsi( this->m_PsiLattice,
- this->m_PsiLattice->GetLargestPossibleRegion() );
- itk::ImageRegionIterator<PointDataImageType> ItPhi( this->m_PhiLattice,
- this->m_PhiLattice->GetLargestPossibleRegion() );
- for ( ItPsi.GoToBegin(), ItPhi.GoToBegin();
- !ItPsi.IsAtEnd(); ++ItPsi, ++ItPhi )
- {
- ItPsi.Set(ItPhi.Get()+ItPsi.Get());
- }
- this->RefineControlLattice();
-
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- if ( this->m_CurrentLevel < this->m_NumberOfLevels[i] )
- {
- this->m_CurrentNumberOfControlPoints[i] =
- 2*this->m_CurrentNumberOfControlPoints[i]-this->m_SplineOrder[i];
- }
- }
-
- itkDebugMacro( << "Current Level = " << this->m_CurrentLevel );
- itkDebugMacro( << " No. of control points = "
- << this->m_CurrentNumberOfControlPoints );
-
- RealType avg_p = 0.0;
- RealType totalWeight = 0.0;
-
- typename PointDataContainerType::Iterator ItIn, ItOut;
- ItIn = this->m_InputPointData->Begin();
- ItOut = this->m_OutputPointData->Begin();
- while ( ItIn != this->m_InputPointData->End() )
- {
- this->m_InputPointData->InsertElement( ItIn.Index(), ItIn.Value()-ItOut.Value() );
- if ( this->GetDebug() )
- {
- RealType weight = this->m_PointWeights->GetElement( ItIn.Index() );
- avg_p += ( ItIn.Value()-ItOut.Value() ).GetNorm() * weight;
- totalWeight += weight;
- }
- ++ItIn;
- ++ItOut;
- }
- itkDebugMacro( << "The average weighted difference norm of the point set is "
- << avg_p / totalWeight );
-
- this->GenerateControlLattice();
- this->UpdatePointSet();
- }
-
- if (this->m_DoMultilevel)
- {
- itk::ImageRegionIterator<PointDataImageType>
- ItPsi( this->m_PsiLattice, this->m_PsiLattice->GetLargestPossibleRegion() );
- itk::ImageRegionIterator<PointDataImageType>
- ItPhi( this->m_PhiLattice, this->m_PhiLattice->GetLargestPossibleRegion() );
- for ( ItPsi.GoToBegin(), ItPhi.GoToBegin(); !ItPsi.IsAtEnd(); ++ItPsi, ++ItPhi )
- {
- ItPsi.Set( ItPhi.Get()+ItPsi.Get() );
- }
-
- typedef itk::ImageDuplicator<PointDataImageType> ImageDuplicatorType;
- typename ImageDuplicatorType::Pointer Duplicator = ImageDuplicatorType::New();
- Duplicator->SetInputImage( this->m_PsiLattice );
- Duplicator->Update();
- this->m_PhiLattice = Duplicator->GetOutput();
- this->UpdatePointSet();
- }
-
- if ( this->m_GenerateOutputImage )
- {
- //this->GenerateOutputImage();
- this->GenerateOutputImageFast();
- }
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::RefineControlLattice()
-{
- ArrayType NumberOfNewControlPoints = this->m_CurrentNumberOfControlPoints;
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- if ( this->m_CurrentLevel < this->m_NumberOfLevels[i] )
- {
- NumberOfNewControlPoints[i] = 2*NumberOfNewControlPoints[i]-this->m_SplineOrder[i];
- }
- }
- typename RealImageType::RegionType::SizeType size;
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- if ( this->m_CloseDimension[i] )
- {
- size[i] = NumberOfNewControlPoints[i] - this->m_SplineOrder[i];
- }
- else
- {
- size[i] = NumberOfNewControlPoints[i];
- }
- }
-
- typename PointDataImageType::Pointer RefinedLattice = PointDataImageType::New();
- RefinedLattice->SetRegions( size );
- RefinedLattice->Allocate();
- PointDataType data;
- data.Fill( 0.0 );
- RefinedLattice->FillBuffer( data );
-
- typename PointDataImageType::IndexType idx;
- typename PointDataImageType::IndexType idx_Psi;
- typename PointDataImageType::IndexType tmp;
- typename PointDataImageType::IndexType tmp_Psi;
- typename PointDataImageType::IndexType off;
- typename PointDataImageType::IndexType off_Psi;
- typename PointDataImageType::RegionType::SizeType size_Psi;
-
- size.Fill( 2 );
- int N = 1;
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- N *= ( this->m_SplineOrder[i]+1 );
- size_Psi[i] = this->m_SplineOrder[i] + 1;
- }
-
- itk::ImageRegionIteratorWithIndex<PointDataImageType>
- It( RefinedLattice, RefinedLattice->GetLargestPossibleRegion() );
-
- It.GoToBegin();
- while ( !It.IsAtEnd() )
- {
- idx = It.GetIndex();
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- if ( this->m_CurrentLevel < this->m_NumberOfLevels[i] )
- {
- idx_Psi[i] = static_cast<unsigned int>( 0.5*idx[i] );
- }
- else
- {
- idx_Psi[i] = static_cast<unsigned int>( idx[i] );
- }
- }
- for ( unsigned int i = 0;
- i < pow( 2.0, static_cast<RealType>( ImageDimension ) ); i++ )
- {
- PointDataType sum( 0.0 );
- PointDataType val;
- off = this->IndexToSubscript( i, size );
-
- bool OutOfBoundary = false;
- for ( unsigned int j = 0; j < ImageDimension; j++ )
- {
- tmp[j] = idx[j] + off[j];
- if ( tmp[j] >= static_cast<long>(NumberOfNewControlPoints[j]) && (!this->m_CloseDimension[j]) )
- {
- OutOfBoundary = true;
- break;
- }
- if ( this->m_CloseDimension[j] )
- {
- tmp[j] %= RefinedLattice->GetLargestPossibleRegion().GetSize()[j];
- }
- }
- if ( OutOfBoundary )
- {
- continue;
- }
-
- for ( unsigned int j = 0; j < static_cast<unsigned int>(N); j++ )
- {
- off_Psi = this->IndexToSubscript( j, size_Psi );
-
- bool OutOfBoundary = false;
- for ( unsigned int k = 0; k < ImageDimension; k++ )
- {
- tmp_Psi[k] = idx_Psi[k] + off_Psi[k];
- if ( tmp_Psi[k] >= static_cast<long>(this->m_CurrentNumberOfControlPoints[k] )
- && (!this->m_CloseDimension[k]) )
- {
- OutOfBoundary = true;
- break;
- }
- if ( this->m_CloseDimension[k] )
- {
- tmp_Psi[k] %= this->m_PsiLattice->GetLargestPossibleRegion().GetSize()[k];
- }
- }
- if ( OutOfBoundary )
- {
- continue;
- }
- RealType coeff = 1.0;
- for ( unsigned int k = 0; k < ImageDimension; k++ )
- {
- coeff *= this->m_RefinedLatticeCoefficients[k]( off[k], off_Psi[k] );
- }
- val = this->m_PsiLattice->GetPixel( tmp_Psi );
- val *= coeff;
- sum += val;
- }
- RefinedLattice->SetPixel( tmp, sum );
- }
-
- bool IsEvenIndex = false;
- while ( !IsEvenIndex && !It.IsAtEnd() )
- {
- ++It;
- idx = It.GetIndex();
- IsEvenIndex = true;
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- if ( idx[i] % 2 )
- {
- IsEvenIndex = false;
- }
- }
- }
- }
-
- typedef itk::ImageDuplicator<PointDataImageType> ImageDuplicatorType;
- typename ImageDuplicatorType::Pointer Duplicator = ImageDuplicatorType::New();
- Duplicator->SetInputImage( RefinedLattice );
- Duplicator->Update();
- this->m_PsiLattice = Duplicator->GetOutput();
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::GenerateControlLattice()
-{
- typename RealImageType::RegionType::SizeType size;
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- if ( this->m_CloseDimension[i] )
- {
- size[i] = this->m_CurrentNumberOfControlPoints[i]-this->m_SplineOrder[i];
- }
- else
- {
- size[i] = this->m_CurrentNumberOfControlPoints[i];
- }
- }
-
- this->m_PhiLattice->SetRegions( size );
- this->m_PhiLattice->Allocate();
- this->m_PhiLattice->FillBuffer( 0.0 );
-
- typename RealImageType::Pointer omega;
- omega = RealImageType::New();
- omega->SetRegions( size );
- omega->Allocate();
- omega->FillBuffer( 0.0 );
-
- typename PointDataImageType::Pointer delta;
- delta = PointDataImageType::New();
- delta->SetRegions( size );
- delta->Allocate();
- delta->FillBuffer( 0.0 );
-
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- size[i] = this->m_SplineOrder[i]+1;
- }
-
- typename RealImageType::Pointer w;
- w = RealImageType::New();
- w->SetRegions( size );
- w->Allocate();
-
- typename PointDataImageType::Pointer phi;
- phi = PointDataImageType::New();
- phi->SetRegions( size );
- phi->Allocate();
-
- itk::ImageRegionIteratorWithIndex<RealImageType>
- Itw( w, w->GetBufferedRegion() );
- itk::ImageRegionIteratorWithIndex<PointDataImageType>
- Itp( phi, phi->GetBufferedRegion() );
-
- vnl_vector<RealType> p( ImageDimension );
- vnl_vector<RealType> r( ImageDimension );
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- r[i] = static_cast<RealType>( this->m_CurrentNumberOfControlPoints[i]-this->m_SplineOrder[i] )
- /static_cast<RealType>( this->m_Size[i]*this->m_Spacing[i]+1e-10 );
- }
-
- typename PointDataContainerType::ConstIterator It;
- for ( It = this->m_InputPointData->Begin(); It != this->m_InputPointData->End(); ++It )
- {
- PointType point;
- this->GetInput()->GetPoint( It.Index(), &point );
-
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- p[i] = ( point[i] - this->m_Origin[i] ) * r[i];
- }
-
- RealType w2_sum = 0.0;
- for ( Itw.GoToBegin(); !Itw.IsAtEnd(); ++Itw )
- {
- RealType B = 1.0;
- typename RealImageType::IndexType idx = Itw.GetIndex();
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- RealType u = static_cast<RealType>( p[i] -
- static_cast<unsigned>( p[i] ) - idx[i] )
- + 0.5*static_cast<RealType>( this->m_SplineOrder[i]-1 );
- B *= this->m_Kernel[i]->Evaluate( u );
- }
- Itw.Set( B );
- w2_sum += B*B;
- }
-
- for ( Itp.GoToBegin(), Itw.GoToBegin(); !Itp.IsAtEnd(); ++Itp, ++Itw )
- {
- typename RealImageType::IndexType idx = Itw.GetIndex();
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- idx[i] += static_cast<unsigned>( p[i] );
- if ( this->m_CloseDimension[i] )
- {
- idx[i] %= delta->GetLargestPossibleRegion().GetSize()[i];
- }
- }
- RealType wc = this->m_PointWeights->GetElement( It.Index() );
- RealType t = Itw.Get();
- PointDataType data = It.Value();
- data *= ( t/w2_sum );
- Itp.Set( data );
- data *= ( t*t*wc );
- delta->SetPixel( idx, delta->GetPixel( idx ) + data );
- omega->SetPixel( idx, omega->GetPixel( idx ) + wc*t*t );
- }
- }
-
- itk::ImageRegionIterator<PointDataImageType>
- Itl( this->m_PhiLattice, this->m_PhiLattice->GetBufferedRegion() );
- itk::ImageRegionIterator<PointDataImageType>
- Itd( delta, delta->GetBufferedRegion() );
- itk::ImageRegionIterator<RealImageType>
- Ito( omega, omega->GetBufferedRegion() );
-
- for ( Itl.GoToBegin(), Ito.GoToBegin(), Itd.GoToBegin();
- !Itl.IsAtEnd(); ++Itl, ++Ito, ++Itd )
- {
- PointDataType P;
- P.Fill( 0 );
- if ( Ito.Get() != 0 )
- {
- P = Itd.Get()/Ito.Get();
- for ( unsigned int i = 0; i < P.Size(); i++ )
- {
- if ( vnl_math_isnan( P[i] ) || vnl_math_isinf( P[i] ) )
- {
- P[i] = 0;
- }
- }
- Itl.Set( P );
- }
- }
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::UpdatePointSet()
-{
- typename PointDataContainerType::ConstIterator ItIn;
-
- ItIn = this->m_InputPointData->Begin();
- while ( ItIn != this->m_InputPointData->End() )
- {
- PointType point;
- PointDataType dataOut;
- point.Fill(0);
-
- this->GetInput()->GetPoint( ItIn.Index(), &point );
- this->EvaluateAtPoint( point, dataOut );
- this->m_OutputPointData->InsertElement( ItIn.Index(), dataOut );
- ++ItIn;
- }
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::GenerateOutputImage()
-{
- itk::ImageRegionIteratorWithIndex<ImageType>
- It( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
-
- for ( It.GoToBegin(); !It.IsAtEnd(); ++It )
- {
- PointDataType data;
- this->EvaluateAtIndex( It.GetIndex(), data );
- It.Set( data );
- }
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::GenerateOutputImageFast()
-{
- typename PointDataImageType::Pointer collapsedPhiLattices[ImageDimension+1];
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- collapsedPhiLattices[i] = PointDataImageType::New();
- collapsedPhiLattices[i]->SetOrigin( this->m_PhiLattice->GetOrigin() );
- collapsedPhiLattices[i]->SetSpacing( this->m_PhiLattice->GetSpacing() );
- typename PointDataImageType::SizeType size;
- size.Fill( 1 );
- for ( unsigned int j = 0; j < i; j++ )
- {
- size[j] = this->m_PhiLattice->GetLargestPossibleRegion().GetSize()[j];
- }
- typename PointDataImageType::IndexType index;
- index.Fill(0);
-
- typename PointDataImageType::RegionType region;
- region.SetIndex(index);
- region.SetSize(size);
-
- collapsedPhiLattices[i]->SetRegions( region );
- collapsedPhiLattices[i]->Allocate();
- collapsedPhiLattices[i]->FillBuffer(typename PointDataImageType::PixelType(0.));
- }
- collapsedPhiLattices[ImageDimension] = this->m_PhiLattice;
- ArrayType totalNumberOfSpans;
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- totalNumberOfSpans[i] = this->m_CurrentNumberOfControlPoints[i] - this->m_SplineOrder[i];
- }
- itk::FixedArray<RealType, ImageDimension> U;
- itk::FixedArray<RealType, ImageDimension> currentU;
- currentU.Fill( -1 );
- typename ImageType::IndexType startIndex
- = this->GetOutput()->GetRequestedRegion().GetIndex();
- typename PointDataImageType::IndexType startPhiIndex
- = this->GetOutput()->GetLargestPossibleRegion().GetIndex();
-
- itk::ImageRegionIteratorWithIndex<ImageType>
- It( this->GetOutput(), this->GetOutput()->GetRequestedRegion() );
- for ( It.GoToBegin(); !It.IsAtEnd(); ++It )
- {
- typename ImageType::IndexType idx = It.GetIndex();
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- U[i] = static_cast<RealType>( totalNumberOfSpans[i] ) *
- static_cast<RealType>( idx[i] - startIndex[i] ) /
- static_cast<RealType>( this->m_Size[i] - 1.0 );
- if ( U[i] == 1.0 )
- {
- U[i] -= vnl_math::float_eps;
- }
- }
- for ( int i = ImageDimension-1; i >= 0; i-- )
- {
- if ( U[i] != currentU[i] )
- {
- for ( int j = i; j >= 0; j-- )
- {
- this->CollapsePhiLattice( collapsedPhiLattices[j+1], collapsedPhiLattices[j], U[j], j );
- currentU[j] = U[j];
- }
- break;
- }
- }
- It.Set( collapsedPhiLattices[0]->GetPixel( startPhiIndex ) );
- }
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::CollapsePhiLattice( PointDataImageType *lattice,
- PointDataImageType *collapsedLattice,
- RealType u, unsigned int dimension )
-{
- itk::ImageRegionIteratorWithIndex<PointDataImageType> It
- ( collapsedLattice, collapsedLattice->GetLargestPossibleRegion() );
-
- for ( It.GoToBegin(); !It.IsAtEnd(); ++It )
- {
- PointDataType data;
- data.Fill( 0.0 );
- typename PointDataImageType::IndexType idx = It.GetIndex();
- for ( unsigned int i = 0; i < this->m_SplineOrder[dimension] + 1; i++ )
- {
- idx[dimension] = static_cast<unsigned int>( u ) + i;
- if ( this->m_CloseDimension[dimension] )
- {
- idx[dimension] %= lattice->GetLargestPossibleRegion().GetSize()[dimension];
- }
- RealType B = this->m_Kernel[dimension]->Evaluate
- ( u - idx[dimension] + 0.5*static_cast<RealType>( this->m_SplineOrder[dimension] - 1 ) );
-
- if(lattice->GetLargestPossibleRegion().IsInside(idx))
- {
- data += ( lattice->GetPixel( idx ) * B );
- }
- }
- It.Set( data );
- }
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::EvaluateAtPoint( PointType point, PointDataType &data )
-{
- for ( unsigned int i = 0; i < ImageDimension; i++)
- {
- point[i] -= this->m_Origin[i];
- point[i] /= ( static_cast<RealType>( ( this->m_Size[i] - 1 )*this->m_Spacing[i] ) );
- }
- this->Evaluate( point, data );
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::EvaluateAtIndex( IndexType idx, PointDataType &data )
-{
- PointType point;
- this->GetOutput()->TransformIndexToPhysicalPoint( idx, point );
- this->EvaluateAtPoint( point, data );
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::EvaluateAtContinuousIndex( ContinuousIndexType idx, PointDataType &data )
-{
- PointType point;
- this->GetOutput()->TransformContinuousIndexToPhysicalPoint( idx, point );
- this->EvaluateAtPoint( point, data );
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::Evaluate( PointType params, PointDataType &data )
-{
- vnl_vector<RealType> p( ImageDimension );
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- if ( params[i] < 0.0 || params[i] > 1.0 )
- {
- itkExceptionMacro( "The specifed point is outside the image domain." );
- }
- if ( params[i] == 1.0 )
- {
- params[i] = 1.0 - vnl_math::float_eps;
- }
- p[i] = static_cast<RealType>( params[i] )
- * static_cast<RealType>( this->m_CurrentNumberOfControlPoints[i]
- - this->m_SplineOrder[i] );
- }
-
- typename RealImageType::RegionType::SizeType size;
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- size[i] = this->m_SplineOrder[i] + 1;
- }
- typename RealImageType::Pointer w;
- w = RealImageType::New();
- w->SetRegions( size );
- w->Allocate();
-
- PointDataType val;
- data.Fill( 0.0 );
-
- itk::ImageRegionIteratorWithIndex<RealImageType> Itw( w, w->GetLargestPossibleRegion() );
-
- for ( Itw.GoToBegin(); !Itw.IsAtEnd(); ++Itw )
- {
- RealType B = 1.0;
- typename RealImageType::IndexType idx = Itw.GetIndex();
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- RealType u = p[i] - static_cast<RealType>( static_cast<unsigned>(p[i]) + idx[i] )
- + 0.5*static_cast<RealType>( this->m_SplineOrder[i] - 1 );
- B *= this->m_Kernel[i]->Evaluate( u );
- }
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- idx[i] += static_cast<unsigned int>( p[i] );
- if ( this->m_CloseDimension[i] )
- {
- idx[i] %= this->m_PhiLattice->GetLargestPossibleRegion().GetSize()[i];
- }
- }
- if ( this->m_PhiLattice->GetLargestPossibleRegion().IsInside( idx ) )
- {
- val = this->m_PhiLattice->GetPixel( idx );
- val *= B;
- data += val;
- }
- }
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::EvaluateGradientAtPoint( PointType point, GradientType &gradient )
-{
- for ( unsigned int i = 0; i < ImageDimension; i++)
- {
- point[i] -= this->m_Origin[i];
- point[i] /= ( static_cast<RealType>( ( this->m_Size[i] - 1 )*this->m_Spacing[i] ) );
- }
- this->EvaluateGradient( point, gradient );
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::EvaluateGradientAtIndex( IndexType idx, GradientType &gradient )
-{
- PointType point;
- this->GetOutput()->TransformIndexToPhysicalPoint( idx, point );
- this->EvaluateGradientAtPoint( point, gradient );
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::EvaluateGradientAtContinuousIndex( ContinuousIndexType idx, GradientType &gradient )
-{
- PointType point;
- this->GetOutput()->TransformContinuousIndexToPhysicalPoint( idx, gradient );
- this->EvaluateGradientAtPoint( point, gradient );
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::EvaluateGradient( PointType params, GradientType &gradient )
-{
- vnl_vector<RealType> p( ImageDimension );
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- if ( params[i] < 0.0 || params[i] > 1.0 )
- {
- itkExceptionMacro( "The specifed point is outside the image domain." );
- }
- if ( params[i] == 1.0 )
- {
- params[i] = 1.0 - vnl_math::float_eps;
- }
- p[i] = static_cast<RealType>( params[i] )
- * static_cast<RealType>( this->m_CurrentNumberOfControlPoints[i]
- - this->m_SplineOrder[i] );
- }
-
- typename RealImageType::RegionType::SizeType size;
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- size[i] = this->m_SplineOrder[i]+1;
- }
- typename RealImageType::Pointer w;
- w = RealImageType::New();
- w->SetRegions( size );
- w->Allocate();
-
- PointDataType val;
- gradient.SetSize( val.Size(), ImageDimension );
- gradient.Fill( 0.0 );
-
- itk::ImageRegionIteratorWithIndex<RealImageType>
- Itw( w, w->GetLargestPossibleRegion() );
-
- for ( unsigned int j = 0; j < gradient.Rows(); j++ )
- {
- for ( Itw.GoToBegin(); !Itw.IsAtEnd(); ++Itw )
- {
- RealType B = 1.0;
- typename RealImageType::IndexType idx = Itw.GetIndex();
- for ( unsigned int k = 0; k < ImageDimension; k++ )
- {
- RealType u = p[k] - static_cast<RealType>( static_cast<unsigned>( p[k] ) + idx[k] )
- + 0.5*static_cast<RealType>( this->m_SplineOrder[k] - 1 );
- if ( j == k )
- {
- B *= this->m_Kernel[k]->EvaluateDerivative( u );
- }
- else
- {
- B *= this->m_Kernel[k]->Evaluate( u );
- }
- }
- for ( unsigned int k = 0; k < ImageDimension; k++ )
- {
- idx[k] += static_cast<unsigned int>( p[k] );
- if ( this->m_CloseDimension[k] )
- {
- idx[k] %= this->m_PhiLattice->GetLargestPossibleRegion().GetSize()[k];
- }
- }
- if ( this->m_PhiLattice->GetLargestPossibleRegion().IsInside( idx ) )
- {
- val = this->m_PhiLattice->GetPixel( idx );
- val *= B;
- for ( unsigned int k = 0; k < val.Size(); k++ )
- {
- gradient( j, k ) += val[k];
- }
- }
- }
- }
-
-}
-
-/**
- * Standard "PrintSelf" method
- */
-template <class TInputPointSet, class TOutputImage>
-void
-BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
-::PrintSelf(
- std::ostream& os,
- itk::Indent indent) const
-{
- Superclass::PrintSelf( os, indent );
- for ( unsigned int i = 0; i < ImageDimension; i++ )
- {
- this->m_Kernel[i]->Print( os, indent );
- }
- os << indent << "B-spline order: "
- << this->m_SplineOrder << std::endl;
- os << indent << "Number Of control points: "
- << this->m_NumberOfControlPoints << std::endl;
- os << indent << "Close dimension: "
- << this->m_CloseDimension << std::endl;
- os << indent << "Number of levels "
- << this->m_NumberOfLevels << std::endl;
-}
-} // end namespace ij
-
-#endif
diff --git a/Utilities/InsightJournal/ijPointSetToImageFilter.h b/Utilities/InsightJournal/ijPointSetToImageFilter.h
deleted file mode 100644
index ce11e8d72e37fd132abaf17404b760a75530c4fa..0000000000000000000000000000000000000000
--- a/Utilities/InsightJournal/ijPointSetToImageFilter.h
+++ /dev/null
@@ -1,154 +0,0 @@
-/*=========================================================================
-
- Program: Insight Segmentation & Registration Toolkit
- Module: $RCSfile: itkPointSetToImageFilter.h,v $
- Language: C++
- Date: $Date: 2005/07/27 15:21:11 $
- Version: $Revision: 1.3 $
-
- Copyright (c) Insight Software Consortium. All rights reserved.
- See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-#ifndef __ijPointSetToImageFilter_h
-#define __ijPointSetToImageFilter_h
-
-#include "itkImageSource.h"
-#include "itkConceptChecking.h"
-
-namespace ij
-{
-
-/** \class PointSetToImageFilter
- * \brief Base class for filters that take a PointSet
- * as input and produce an image as output.
- * By default, if the user does not specify the size of the output image,
- * the maximum size of the point-set's bounding box is used.
- */
-template <class TInputPointSet, class TOutputImage>
-class ITK_EXPORT PointSetToImageFilter : public itk::ImageSource<TOutputImage>
-{
-public:
- /** Standard class typedefs. */
- typedef PointSetToImageFilter Self;
- typedef itk::ImageSource<TOutputImage> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
- typedef typename TOutputImage::SizeType SizeType;
- typedef TOutputImage OutputImageType;
- typedef typename OutputImageType::Pointer OutputImagePointer;
- typedef typename OutputImageType::ValueType ValueType;
-
- /** Method for creation through the object factory. */
- itkNewMacro(Self);
-
- /** Run-time type information (and related methods). */
- itkTypeMacro(PointSetToImageFilter,ImageSource);
-
- /** Superclass typedefs. */
- typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
-
- /** Some convenient typedefs. */
- typedef TInputPointSet InputPointSetType;
- typedef typename InputPointSetType::Pointer InputPointSetPointer;
- typedef typename InputPointSetType::ConstPointer InputPointSetConstPointer;
-
- /** Dimension constants */
- itkStaticConstMacro(InputPointSetDimension, unsigned int,
- InputPointSetType::PointDimension);
- itkStaticConstMacro(OutputImageDimension, unsigned int,
- TOutputImage::ImageDimension);
-
- /** Image spacing and origin typedefs */
- typedef typename TOutputImage::SpacingType SpacingType;
- typedef typename TOutputImage::PointType PointType;
-
- /** Set/Get the input point-set of this process object. */
- virtual void SetInput( const InputPointSetType *pointset);
- virtual void SetInput( unsigned int, const InputPointSetType * pointset);
- const InputPointSetType * GetInput(void);
- const InputPointSetType * GetInput(unsigned int idx);
-
- /** Set the spacing (size of a pixel) of the image. The
- * spacing is the geometric distance between image samples.
- * It is stored internally as double, but may be set from
- * float. \sa GetSpacing() */
- itkSetMacro(Spacing,SpacingType);
- virtual void SetSpacing( const double* spacing);
- virtual void SetSpacing( const float* spacing);
-
- /** Get the spacing (size of a pixel) of the image. The
- * spacing is the geometric distance between image samples.
- * The value returned is a pointer to a double array.
- * For ImageBase and Image, the default data spacing is unity. */
- itkGetConstReferenceMacro(Spacing,SpacingType);
-
- /** Set the origin of the image. The origin is the geometric
- * coordinates of the image origin. It is stored internally
- * as double but may be set from float.
- * \sa GetOrigin() */
- itkSetMacro(Origin,PointType);
- virtual void SetOrigin( const double* origin);
- virtual void SetOrigin( const float* origin);
-
- /** Get the origin of the image. The origin is the geometric
- * coordinates of the index (0,0). The value returned is a pointer
- * to a double array. For ImageBase and Image, the default origin is
- * 0. */
- itkGetConstReferenceMacro(Origin,PointType);
-
- /** Set/Get the value for pixels in the point-set.
- * By default, this filter will return an image
- * that contains values from the point-set specified as input.
- * If this "inside" value is changed to a non-null value,
- * the output produced by this filter will be a mask with inside/outside values
- * specified by the user. */
- itkSetMacro(InsideValue, ValueType);
- itkGetMacro(InsideValue, ValueType);
-
- /** Set/Get the value for pixels outside the point-set.
- * By default, this filter will return an image
- * that contains values from the point specified as input.
- * If this "outside" value is changed to a non-null value,
- * the output produced by this filter will be a mask with inside/outside values
- * specified by the user. */
- itkSetMacro(OutsideValue, ValueType);
- itkGetMacro(OutsideValue, ValueType);
-
- /** Set/Get Size */
- itkSetMacro(Size,SizeType);
- itkGetMacro(Size,SizeType);
-
-protected:
- PointSetToImageFilter();
- ~PointSetToImageFilter();
-
- virtual void GenerateOutputInformation(){}; // do nothing
- virtual void GenerateData();
-
- SizeType m_Size;
- SpacingType m_Spacing;
- PointType m_Origin;
- ValueType m_InsideValue;
- ValueType m_OutsideValue;
-
- virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
-
-private:
- PointSetToImageFilter(const Self&); //purposely not implemented
- void operator=(const Self&); //purposely not implemented
-
-
-};
-
-} // end namespace itk
-
-#ifndef OTB_MANUAL_INSTANTIATION
-#include "ijPointSetToImageFilter.txx"
-#endif
-
-#endif
diff --git a/Utilities/InsightJournal/ijPointSetToImageFilter.txx b/Utilities/InsightJournal/ijPointSetToImageFilter.txx
deleted file mode 100644
index 5eeb9774015af5a417dfe99768e8f02e831bb46d..0000000000000000000000000000000000000000
--- a/Utilities/InsightJournal/ijPointSetToImageFilter.txx
+++ /dev/null
@@ -1,282 +0,0 @@
-/*=========================================================================
-
-Program: Insight Segmentation & Registration Toolkit
-Module: $RCSfile: itkPointSetToImageFilter.txx,v $
-Language: C++
-Date: $Date: 2005/07/27 15:21:11 $
-Version: $Revision: 1.11 $
-
-Copyright (c) Insight Software Consortium. All rights reserved.
-See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
-
-This software is distributed WITHOUT ANY WARRANTY; without even
-the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
-PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-#ifndef _ijPointSetToImageFilter_txx
-#define _ijPointSetToImageFilter_txx
-
-#include "ijPointSetToImageFilter.h"
-#include <itkImageRegionIteratorWithIndex.h>
-#include <itkNumericTraits.h>
-
-namespace ij
-{
-
-/** Constructor */
-template <class TInputPointSet, class TOutputImage>
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::PointSetToImageFilter()
-{
- this->SetNumberOfRequiredInputs(1);
- m_Size.Fill(0);
- m_Spacing.Fill(1.0);
- m_Origin.Fill(0.0);
-/*
- m_InsideValue = 1;
- m_OutsideValue = 0;
-*/
-}
-
-/** Destructor */
-template <class TInputPointSet, class TOutputImage>
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::~PointSetToImageFilter()
-{
-}
-
-
-/** Set the Input PointSet */
-template <class TInputPointSet, class TOutputImage>
-void
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::SetInput(const InputPointSetType *input)
-{
- // Process object is not const-correct so the const_cast is required here
- this->itk::ProcessObject::SetNthInput(0,
- const_cast< InputPointSetType * >( input ) );
-}
-
-
-/** Connect one of the operands */
-template <class TInputPointSet, class TOutputImage>
-void
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::SetInput( unsigned int index, const TInputPointSet * pointset )
-{
- // Process object is not const-correct so the const_cast is required here
- this->itk::ProcessObject::SetNthInput(index,
- const_cast< TInputPointSet *>( pointset ) );
-}
-
-
-
-/** Get the input point-set */
-template <class TInputPointSet, class TOutputImage>
-const typename PointSetToImageFilter<TInputPointSet,TOutputImage>::InputPointSetType *
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::GetInput(void)
-{
- if (this->GetNumberOfInputs() < 1)
- {
- return 0;
- }
-
- return static_cast<const TInputPointSet * >
- (this->itk::ProcessObject::GetInput(0) );
-}
-
-/** Get the input point-set */
-template <class TInputPointSet, class TOutputImage>
-const typename PointSetToImageFilter<TInputPointSet,TOutputImage>::InputPointSetType *
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::GetInput(unsigned int idx)
-{
- return static_cast< const TInputPointSet * >
- (this->itk::ProcessObject::GetInput(idx));
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::SetSpacing(const float* v)
-{
- itk::Vector<float, OutputImageDimension> vf(v);
- SpacingType spacing;
- spacing.CastFrom(vf);
- this->SetSpacing(spacing);
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::SetSpacing(const double* v)
-{
- SpacingType spacing(v);
- this->SetSpacing(spacing);
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::SetOrigin(const float* v)
-{
- itk::Point<float,OutputImageDimension> pf(v);
- PointType origin;
- origin.CastFrom(pf);
- this->SetOrigin(origin);
-}
-
-template <class TInputPointSet, class TOutputImage>
-void
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::SetOrigin(const double* v)
-{
- PointType origin(v);
- this->SetOrigin(origin);
-}
-
-//----------------------------------------------------------------------------
-
-/** Update */
-template <class TInputPointSet, class TOutputImage>
-void
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::GenerateData(void)
-{
- unsigned int i;
- itkDebugMacro(<< "PointSetToImageFilter::Update() called");
-
- // Get the input and output pointers
- const InputPointSetType * InputPointSet = this->GetInput();
- OutputImagePointer OutputImage = this->GetOutput();
-
- // Generate the image
- double origin[InputPointSetDimension];
- SizeType size;
-
- typedef typename InputPointSetType::BoundingBoxType BoundingBoxType;
-
- const BoundingBoxType * bb = InputPointSet->GetBoundingBox();
-
- for(i=0;i<InputPointSetDimension;i++)
- {
- size[i] = (unsigned long)(bb->GetBounds()[2*i+1]-bb->GetBounds()[2*i]);
- origin[i]= 0; //bb->GetBounds()[2*i];
- }
-
- typename OutputImageType::RegionType region;
-
- // If the size of the output has been explicitly specified, the filter
- // will set the output size to the explicit size, otherwise the size from the spatial
- // PointSet's bounding box will be used as default.
-
- bool specified = false;
- for (i = 0; i < OutputImageDimension; i++)
- {
- if (m_Size[i] != 0)
- {
- specified = true;
- break;
- }
- }
-
- if (specified)
- {
- region.SetSize( m_Size );
- }
- else
- {
- region.SetSize( size );
- }
-
- OutputImage->SetRegions( region);
-
- // If the spacing has been explicitly specified, the filter
- // will set the output spacing to that explicit spacing, otherwise the spacing from
- // the point-set is used as default.
-
- specified = false;
- for (i = 0; i < OutputImageDimension; i++)
- {
- if (m_Spacing[i] != 0)
- {
- specified = true;
- break;
- }
- }
-
- if (specified)
- {
- OutputImage->SetSpacing(this->m_Spacing); // set spacing
- }
-
-
- specified = false;
- for (i = 0; i < OutputImageDimension; i++)
- {
- if (m_Origin[i] != 0)
- {
- specified = true;
- break;
- }
- }
-
- if (specified)
- {
- for (i = 0; i < OutputImageDimension; i++)
- {
- origin[i] = m_Origin[i]; // set origin
- }
- }
-
- OutputImage->SetOrigin(origin); // and origin
- OutputImage->Allocate(); // allocate the image
-
- OutputImage->FillBuffer(m_OutsideValue);
-
- typedef typename InputPointSetType::PointsContainer::ConstIterator PointIterator;
- PointIterator pointItr = InputPointSet->GetPoints()->Begin();
- PointIterator pointEnd = InputPointSet->GetPoints()->End();
-
- typename OutputImageType::IndexType index;
-
- while( pointItr != pointEnd )
- {
- if(OutputImage->TransformPhysicalPointToIndex(pointItr.Value(),index))
- {
- OutputImage->SetPixel(index,m_InsideValue);
- }
- pointItr++;
- }
-
- itkDebugMacro(<< "PointSetToImageFilter::Update() finished");
-
-
-} // end update function
-
-
-template<class TInputPointSet, class TOutputImage>
-void
-PointSetToImageFilter<TInputPointSet,TOutputImage>
-::PrintSelf(std::ostream& os, itk::Indent indent) const
-{
- Superclass::PrintSelf(os, indent);
- os << indent << "Size : " << m_Size << std::endl;
- os << indent << "Origin: " << m_Origin << std::endl;
- os << indent << "Spacing: " << m_Spacing << std::endl;
- os << indent << "Inside Value : "
- << static_cast<typename itk::NumericTraits<ValueType>::PrintType>(m_InsideValue)
- << std::endl;
- os << indent << "Outside Value : "
- << static_cast<typename itk::NumericTraits<ValueType>::PrintType>(m_OutsideValue)
- << std::endl;
-
-}
-
-
-
-} // end namespace itk
-
-#endif