From fafcf4aacc8997f06a0bc9bf80db3bfd148f5e22 Mon Sep 17 00:00:00 2001
From: Julien Michel <julien.michel@c-s.fr>
Date: Wed, 16 May 2007 08:48:58 +0000
Subject: [PATCH] =?UTF-8?q?Int=C3=A9gration=20du=20code=20d'interpolation?=
=?UTF-8?q?=20par=20splines=20de=20points=20=C3=A9parts,=20contribution=20?=
=?UTF-8?q?de=20Insight=20Journal.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
Testing/Utilities/CMakeLists.txt | 39 +-
...ScatteredDataPointSetToImageFilterTest.cxx | 122 +++
Testing/Utilities/ossimIntegrationTest.cxx | 4 +-
Testing/Utilities/otbUtilitiesTests.cxx | 31 +
Utilities/CMakeLists.txt | 2 +-
Utilities/InsightJournal/CMakeLists.txt | 18 +
Utilities/InsightJournal/foo.cxx | 0
.../InsightJournal/ijBSplineKernelFunction.h | 175 ++++
.../ijBSplineKernelFunction.txx | 185 ++++
...SplineScatteredDataPointSetToImageFilter.h | 228 +++++
...lineScatteredDataPointSetToImageFilter.txx | 854 ++++++++++++++++++
.../InsightJournal/ijPointSetToImageFilter.h | 154 ++++
.../ijPointSetToImageFilter.txx | 282 ++++++
otbIncludeDirectories.cmake | 2 +
14 files changed, 2089 insertions(+), 7 deletions(-)
create mode 100755 Testing/Utilities/ijBSplineScatteredDataPointSetToImageFilterTest.cxx
create mode 100644 Testing/Utilities/otbUtilitiesTests.cxx
create mode 100755 Utilities/InsightJournal/CMakeLists.txt
create mode 100644 Utilities/InsightJournal/foo.cxx
create mode 100755 Utilities/InsightJournal/ijBSplineKernelFunction.h
create mode 100755 Utilities/InsightJournal/ijBSplineKernelFunction.txx
create mode 100755 Utilities/InsightJournal/ijBSplineScatteredDataPointSetToImageFilter.h
create mode 100755 Utilities/InsightJournal/ijBSplineScatteredDataPointSetToImageFilter.txx
create mode 100755 Utilities/InsightJournal/ijPointSetToImageFilter.h
create mode 100755 Utilities/InsightJournal/ijPointSetToImageFilter.txx
diff --git a/Testing/Utilities/CMakeLists.txt b/Testing/Utilities/CMakeLists.txt
index 7b64214725..84204d7835 100644
--- a/Testing/Utilities/CMakeLists.txt
+++ b/Testing/Utilities/CMakeLists.txt
@@ -1,14 +1,45 @@
IF( NOT OTB_DISABLE_CXX_TESTING )
-# ------- ossim integration test ----------
+SET(BASELINE ${OTB_DATA_ROOT}/Baseline/OTB/Images)
+SET(BASELINE_FILES ${OTB_DATA_ROOT}/Baseline/OTB/Files)
+SET(INPUTDATA ${OTB_DATA_ROOT}/Input)
+#Images de teledetection (grosses images )
+SET(IMAGEDATA ${OTB_DATA_ROOT}/LargeInput )
+SET(TEMP ${OTBTesting_BINARY_DIR}/Temporary)
+
+#Tolerance sur diff pixel image
+SET(TOL 0.0)
+
+
+SET(UTILITIES_TESTS ${CXX_TEST_PATH}/otbUtilitiesTests)
+
+
+
+ADD_TEST(utOssimIntegrationTest ${UTILITIES_TESTS}
+ossimIntegrationTest
+)
+
+ADD_TEST(utIjBSplineScatteredDataPointSetToImageFilterTest ${UTILITIES_TESTS}
+--compare-image ${TOL}
+ ${BASELINE}/ijBSplineScatteredDataPointSetToImageFilterTestOutput.hdr
+ ${TEMP}/ijBSplineScatteredDataPointSetToImageFilterTestOutput.hdr
+ ijBSplineScatteredDataPointSetToImageFilterTest
+ ${INPUTDATA}/brain.png
+ ${TEMP}/ijBSplineScatteredDataPointSetToImageFilterTestOutput.hdr
+)
-ADD_TEST(utTuOSSIMIntegrationTest ${CXX_TEST_PATH}/utOSSIMIntegrationTest)
# ------- Fichiers sources CXX -----------------------------------
+SET(UtilitiesTests_SRCS
+ossimIntegrationTest.cxx
+ijBSplineScatteredDataPointSetToImageFilterTest.cxx
+)
+
INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}")
-ADD_EXECUTABLE(utOSSIMIntegrationTest ossimIntegrationTest.cxx)
-TARGET_LINK_LIBRARIES(utOSSIMIntegrationTest gdal ossim)
+ADD_EXECUTABLE(otbUtilitiesTests otbUtilitiesTests.cxx ${UtilitiesTests_SRCS})
+TARGET_LINK_LIBRARIES(otbUtilitiesTests OTBIO OTBCommon gdal ITKIO ITKAlgorithms ITKStatistics ITKCommon ossim)
+
ENDIF( NOT OTB_DISABLE_CXX_TESTING )
diff --git a/Testing/Utilities/ijBSplineScatteredDataPointSetToImageFilterTest.cxx b/Testing/Utilities/ijBSplineScatteredDataPointSetToImageFilterTest.cxx
new file mode 100755
index 0000000000..c528c5e2c7
--- /dev/null
+++ b/Testing/Utilities/ijBSplineScatteredDataPointSetToImageFilterTest.cxx
@@ -0,0 +1,122 @@
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkPointSet.h"
+#include "ijBSplineScatteredDataPointSetToImageFilter.h"
+
+/**
+ * In this test, we approximate a 2-D scalar field.
+ * The scattered data is derived from a segmented
+ * image. We write the output to an image for
+ * comparison.
+ */
+int ijBSplineScatteredDataPointSetToImageFilterTest( int argc, char *argv[] )
+{
+ const unsigned int ParametricDimension = 2;
+ const unsigned int DataDimension = 1;
+
+ typedef int PixelType;
+ typedef otb::Image<PixelType, ParametricDimension> InputImageType;
+ typedef float RealType;
+ typedef itk::Vector<RealType, DataDimension> VectorType;
+ typedef otb::Image<VectorType, ParametricDimension> VectorImageType;
+ typedef itk::PointSet
+ <VectorImageType::PixelType, ParametricDimension> PointSetType;
+ PointSetType::Pointer pointSet = PointSetType::New();
+
+ typedef otb::ImageFileReader<InputImageType> ReaderType;
+ ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName( argv[1] );
+ reader->Update();
+
+ itk::ImageRegionIteratorWithIndex<InputImageType>
+ It( reader->GetOutput(), reader->GetOutput()->GetLargestPossibleRegion() );
+
+ // Iterate through the input image which consists of multivalued
+ // foreground pixels (=nonzero) and background values (=zero).
+ // The foreground pixels comprise the input point set.
+
+ for ( It.GoToBegin(); !It.IsAtEnd(); ++It )
+ {
+ if ( It.Get() != itk::NumericTraits<PixelType>::Zero )
+ {
+ // We extract both the 2-D location of the point
+ // and the pixel value of that point.
+
+ PointSetType::PointType point;
+ reader->GetOutput()->TransformIndexToPhysicalPoint( It.GetIndex(), point );
+
+ unsigned long i = pointSet->GetNumberOfPoints();
+ pointSet->SetPoint( i, point );
+
+ PointSetType::PixelType V( DataDimension );
+ V[0] = static_cast<RealType>( It.Get() );
+ pointSet->SetPointData( i, V );
+ }
+ }
+
+ // Instantiate the B-spline filter and set the desired parameters.
+ typedef ij::BSplineScatteredDataPointSetToImageFilter
+ <PointSetType, VectorImageType> FilterType;
+ FilterType::Pointer filter = FilterType::New();
+ filter->SetSplineOrder( 3 );
+ FilterType::ArrayType ncps;
+ ncps.Fill( 4 );
+ filter->SetNumberOfControlPoints( ncps );
+ filter->SetNumberOfLevels( 10 );
+
+ // Define the parametric domain.
+ filter->SetOrigin( reader->GetOutput()->GetOrigin() );
+ filter->SetSpacing( reader->GetOutput()->GetSpacing() );
+ filter->SetSize( reader->GetOutput()->GetLargestPossibleRegion().GetSize() );
+
+ filter->SetInput( pointSet );
+
+ try
+ {
+ filter->Update();
+ }
+ catch (...)
+ {
+ std::cerr << "Test 1: itkBSplineScatteredDataImageFilter exception thrown"
+ << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ // Write the output to an image.
+ typedef otb::Image<RealType, ParametricDimension> RealImageType;
+ RealImageType::Pointer image = RealImageType::New();
+ image->SetRegions( reader->GetOutput()->GetLargestPossibleRegion() );
+ image->Allocate();
+ itk::ImageRegionIteratorWithIndex<RealImageType>
+ Itt( image, image->GetLargestPossibleRegion() );
+
+ for ( Itt.GoToBegin(); !Itt.IsAtEnd(); ++Itt )
+ {
+ Itt.Set( filter->GetOutput()->GetPixel( Itt.GetIndex() )[0] );
+ }
+
+ typedef otb::ImageFileWriter<RealImageType> WriterType;
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetInput( image );
+ writer->SetFileName( argv[2] );
+ writer->Update();
+
+ return EXIT_SUCCESS;
+};
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Testing/Utilities/ossimIntegrationTest.cxx b/Testing/Utilities/ossimIntegrationTest.cxx
index 3dd4ee80f6..890fb74a56 100644
--- a/Testing/Utilities/ossimIntegrationTest.cxx
+++ b/Testing/Utilities/ossimIntegrationTest.cxx
@@ -50,7 +50,7 @@ void printOutputTypes();
ossimProjection* newUtmView(const ossimGpt& centerGround,
const ossimDpt& metersPerPixel);
-int main(int argc, char* argv[])
+int ossimIntegrationTest(int argc, char* argv[])
{
ossimInit::instance()->initialize(argc, argv);
@@ -163,7 +163,7 @@ ossimProjection* newUtmView(const ossimGpt& centerGround,
ossimUtmProjection* utm = new ossimUtmProjection;
// we will make it a square pixel in meters
- double averageGsd = (metersPerPixel.x + metersPerPixel.y)*.5;
+ //double averageGsd = (metersPerPixel.x + metersPerPixel.y)*.5;
utm->setZone(centerGround);
utm->setMetersPerPixel(ossimDpt(metersPerPixel));
diff --git a/Testing/Utilities/otbUtilitiesTests.cxx b/Testing/Utilities/otbUtilitiesTests.cxx
new file mode 100644
index 0000000000..1cb56926fc
--- /dev/null
+++ b/Testing/Utilities/otbUtilitiesTests.cxx
@@ -0,0 +1,31 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt 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.
+
+=========================================================================*/
+
+// this file defines the otbCommonTest for the test driver
+// and all it expects is that you have a function called RegisterTests
+#if defined(_MSC_VER)
+#pragma warning ( disable : 4786 )
+#endif
+#include <iostream>
+#include "otbTestMain.h"
+
+void RegisterTests()
+{
+REGISTER_TEST(ossimIntegrationTest);
+REGISTER_TEST(ijBSplineScatteredDataPointSetToImageFilterTest);
+}
diff --git a/Utilities/CMakeLists.txt b/Utilities/CMakeLists.txt
index 7633047034..4d130a0b76 100755
--- a/Utilities/CMakeLists.txt
+++ b/Utilities/CMakeLists.txt
@@ -4,7 +4,7 @@ IF(NOT OTB_USE_EXTERNAL_ITK)
SUBDIRS( ITK )
ENDIF(NOT OTB_USE_EXTERNAL_ITK)
-SUBDIRS(BGL otbsvm dxflib)
+SUBDIRS(BGL otbsvm dxflib InsightJournal)
OPTION(OTB_OSSIM_COMPILATION "EXPERIMENTAL: Compile OSSIM library embedded within OTB" OFF)
MARK_AS_ADVANCED(OTB_OSSIM_COMPILATION)
diff --git a/Utilities/InsightJournal/CMakeLists.txt b/Utilities/InsightJournal/CMakeLists.txt
new file mode 100755
index 0000000000..a7b82a2774
--- /dev/null
+++ b/Utilities/InsightJournal/CMakeLists.txt
@@ -0,0 +1,18 @@
+
+# Sources of non-templated classes.
+FILE(GLOB InsightJournal_SRCS "*.cxx" )
+
+ADD_LIBRARY(InsightJournal ${InsightJournal_SRCS})
+TARGET_LINK_LIBRARIES (InsightJournal ITKCommon)
+
+INSTALL(TARGETS InsightJournal
+RUNTIME DESTINATION ${OTB_INSTALL_BIN_DIR} COMPONENT RuntimeLibraries
+LIBRARY DESTINATION ${OTB_INSTALL_LIB_DIR} COMPONENT RuntimeLibraries
+ARCHIVE DESTINATION ${OTB_INSTALL_LIB_DIR} COMPONENT Development)
+
+FILE(GLOB __files1 "${CMAKE_CURRENT_SOURCE_DIR}/*.h")
+FILE(GLOB __files2 "${CMAKE_CURRENT_SOURCE_DIR}/*.txx")
+
+INSTALL(FILES ${__files1} ${__files2}
+ DESTINATION ${OTB_INSTALL_INCLUDE_DIR}/Utilities/InsightJournal
+ COMPONENT Development)
diff --git a/Utilities/InsightJournal/foo.cxx b/Utilities/InsightJournal/foo.cxx
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/Utilities/InsightJournal/ijBSplineKernelFunction.h b/Utilities/InsightJournal/ijBSplineKernelFunction.h
new file mode 100755
index 0000000000..121e0f34f1
--- /dev/null
+++ b/Utilities/InsightJournal/ijBSplineKernelFunction.h
@@ -0,0 +1,175 @@
+/*=========================================================================
+
+ 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 < 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
new file mode 100755
index 0000000000..7f82b08574
--- /dev/null
+++ b/Utilities/InsightJournal/ijBSplineKernelFunction.txx
@@ -0,0 +1,185 @@
+/*=========================================================================
+
+ 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, b;
+ 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
new file mode 100755
index 0000000000..6634d4b5ec
--- /dev/null
+++ b/Utilities/InsightJournal/ijBSplineScatteredDataPointSetToImageFilter.h
@@ -0,0 +1,228 @@
+/*=========================================================================
+
+ 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 "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 ij::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 vnl_matrix<RealType> GradientType;
+ typedef itk::Image<RealType,
+ itkGetStaticConstMacro( ImageDimension )> RealImageType;
+ typedef itk::Image<PointDataType,
+ itkGetStaticConstMacro( ImageDimension )> PointDataImageType;
+
+ 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();
+
+ /** Interpolation kernel type (default spline order = 3) */
+ typedef BSplineKernelFunction<3> KernelType;
+
+ 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
new file mode 100755
index 0000000000..bc2982ae88
--- /dev/null
+++ b/Utilities/InsightJournal/ijBSplineScatteredDataPointSetToImageFilter.txx
@@ -0,0 +1,854 @@
+#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 "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()
+{
+
+ Superclass::GenerateData();
+
+ this->m_InputPointData->Initialize();
+ this->m_OutputPointData->Initialize();
+
+ if ( this->m_UsePointWeights &&
+ ( this->m_PointWeights->Size() != this->GetInput()->GetNumberOfPoints() ) )
+ {
+ itkExceptionMacro( << "The number of weight points and input points must be equal." );
+ }
+
+ 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 = itk::NumericTraits<RealType>::Zero;
+ // RealType norm_p;
+
+ 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() )
+ {
+ avg_p += ( ItIn.Value()-ItOut.Value() ).GetNorm();
+ }
+ ++ItIn;
+ ++ItOut;
+ }
+ itkDebugMacro( << "The average difference norm of the point set is "
+ << avg_p / static_cast<RealType>( this->GetInput()->GetNumberOfPoints() ) );
+
+ 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();
+ }
+}
+
+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] >= 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 < 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] >= 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;
+
+ 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 )
+ {
+ typename ImageType::IndexType idx = It.GetIndex();
+ PointType point;
+ for ( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ point[i] = static_cast<RealType>( idx[i] )
+ /( static_cast<RealType>( this->m_Size[i]-1 ) );
+ }
+ PointDataType data;
+ this->Evaluate( point, data );
+ 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]*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] >= 1.0 )
+ {
+ params[i] = 0.99999;
+ }
+ if ( params[i] < 0.0 )
+ {
+ params[i] = 0.0;
+ }
+ 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();
+
+ itk::ImageRegionIteratorWithIndex<RealImageType>
+ Itw( w, w->GetLargestPossibleRegion() );
+
+ PointDataType val;
+ data.Fill( 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 = 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]*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] >= 1.0 )
+ {
+ params[i] = 0.99999;
+ }
+ if ( params[i] < 0.0 )
+ {
+ params[i] = 0.0;
+ }
+ 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();
+
+ itk::ImageRegionIteratorWithIndex<RealImageType>
+ Itw( w, w->GetLargestPossibleRegion() );
+
+ PointDataType val;
+ gradient.set_size( val.Size(), ImageDimension );
+ gradient.fill( 0.0 );
+
+ for ( unsigned int j = 0; j < gradient.cols(); 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( k, j ) += 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
new file mode 100755
index 0000000000..ce11e8d72e
--- /dev/null
+++ b/Utilities/InsightJournal/ijPointSetToImageFilter.h
@@ -0,0 +1,154 @@
+/*=========================================================================
+
+ 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
new file mode 100755
index 0000000000..5eeb977401
--- /dev/null
+++ b/Utilities/InsightJournal/ijPointSetToImageFilter.txx
@@ -0,0 +1,282 @@
+/*=========================================================================
+
+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
diff --git a/otbIncludeDirectories.cmake b/otbIncludeDirectories.cmake
index 9ff14dfbd1..2b272f2ae4 100644
--- a/otbIncludeDirectories.cmake
+++ b/otbIncludeDirectories.cmake
@@ -27,6 +27,7 @@ SET(OTB_INCLUDE_DIRS_BUILD_TREE ${OTB_INCLUDE_DIRS_BUILD_TREE}
${OTB_SOURCE_DIR}/Utilities/OSSIM/include
${OTB_SOURCE_DIR}/Utilities/OSSIM/include/ossim
${OTB_SOURCE_DIR}/Utilities/dxflib
+ ${OTB_SOURCE_DIR}/Utilities/InsightJournal
)
#-----------------------------------------------------------------------------
@@ -101,6 +102,7 @@ SET(OTB_INCLUDE_DIRS_INSTALL_TREE ${OTB_INCLUDE_DIRS_INSTALL_TREE}
${OTB_INSTALL_INCLUDE_DIR}/Utilities/BGL
${OTB_INSTALL_INCLUDE_DIR}/Utilities/BGL/boost
${OTB_INSTALL_INCLUDE_DIR}/Utilities/otbsvm
+ ${OTB_INSTALL_INCLUDE_DIR}/Utilities/InsightJournal
)
SET(OTB_INCLUDE_DIRS_INSTALL_TREE ${OTB_INCLUDE_DIRS_INSTALL_TREE}
--
GitLab