diff --git a/Code/BasicFilters/otbPrintableImageFilter.h b/Code/BasicFilters/otbPrintableImageFilter.h
index 80d2687c227b93a5ffde71579a53d4105d771aff..9025dc2fb8813d4ea36868627a82427f60e2d4ee 100644
--- a/Code/BasicFilters/otbPrintableImageFilter.h
+++ b/Code/BasicFilters/otbPrintableImageFilter.h
@@ -22,77 +22,231 @@
#ifndef __otbPrintableImageFilter_h
#define __otbPrintableImageFilter_h
+#include "otbImage.h"
#include "itkImageToImageFilter.h"
#include "otbVectorRescaleIntensityImageFilter.h"
#include "otbMultiChannelExtractROI.h"
+#include "itkBinaryFunctorImageFilter.h"
namespace otb
{
+namespace Functor
+{
/**
- * \class PrintableImageFilter
- * \brief This class is a helper class to turn a vector image to a generic 8 bytes RGB image.
- *
- * It is useful for publications for instance.
- **/
+ * \class MaskFunctor
+ * \brief Output is a InputPixel if MaskPixel is m_Background a defined other value (m_ObjkeectColor) else.
+ */
+
+template<class TInputPixel, class TMaskPixel, class TOutputPixel>
+class ITK_EXPORT MaskFunctor
+{
+ public :
+ MaskFunctor()
+ {
+ m_BackgroundValue = 0;
+ m_ObjectColor.SetSize(3);
+ m_ObjectColor.Fill(255);
+ };
+ ~MaskFunctor(){};
+
+ typedef TInputPixel InputPixelType;
+ typedef TMaskPixel MaskPixelType;
+ typedef TOutputPixel OutputPixelType;
+ typedef typename OutputPixelType::ValueType OutputInternalPixelType;
+
+ MaskPixelType GetBackgroundValue()
+ {
+ return m_BackgroundValue;
+ }
+ void SetBackgroundValue( MaskPixelType val )
+ {
+ m_BackgroundValue = val;
+ }
+
+ OutputPixelType GetObjectColor()
+ {
+ return m_ObjectColor;
+ }
+ void SetObjectColor( OutputPixelType val )
+ {
+ m_ObjectColor = val;
+ }
+
+ inline OutputPixelType operator()(InputPixelType inPix, MaskPixelType maskPix) const
+ {
+ OutputPixelType outPix;
+ if( maskPix==m_BackgroundValue )
+ {
+ outPix.SetSize( inPix.Size() );
+ for(unsigned int i=0; i<outPix.Size(); i++)
+ {
+ outPix[i] = static_cast<OutputInternalPixelType>(inPix[i]);
+ }
+ }
+ else
+ {
+ outPix = m_ObjectColor;
+ }
+
+ return outPix;
+ }
+
+ protected:
+ MaskPixelType m_BackgroundValue;
+ OutputPixelType m_ObjectColor;
+};
+}
-template <class TInputImage>
+
+/**
+ * \class PrintableImageFilter
+ * \brief This class is a helper class to turn a vector image to a generic 8 bytes RGB image.
+ * A mask can be used to highlight some object represebted by the same value.
+ * The mask is a binary image. Background MaskValue is used to precise which
+ * value of the mask are objects (default 0).
+ * Output object color can be set using m_ObjectColor (default white).
+ * The output is a 3 channels image, each channel is a channel of the input image.
+ * They can be selected usin m_ChannelList or SetChannel(int ch ) method.
+ *
+ * It is useful for publications for instance.
+ *
+ * \sa itkImageToImageFilter
+ **/
+
+template <class TInputImage , class TMaskImage = otb::Image<unsigned char, 2> >
class ITK_EXPORT PrintableImageFilter :
- public itk::ImageToImageFilter<TInputImage, otb::VectorImage<unsigned char,2> >
+public itk::ImageToImageFilter<TInputImage, otb::VectorImage<unsigned char, 2> >
{
-public:
+ public:
typedef PrintableImageFilter Self;
typedef itk::ImageToImageFilter
- <TInputImage, otb::VectorImage<unsigned char,2> > Superclass;
+ <TInputImage, otb::VectorImage<unsigned char,2> > Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
-
+
typedef TInputImage InputImageType;
typedef typename InputImageType::PixelType InputPixelType;
- typedef unsigned char OutputPixelType;
- typedef otb::VectorImage<OutputPixelType,2> OutputImageType;
+ typedef typename InputImageType::InternalPixelType InputInternalPixelType;
+ typedef unsigned char OutputInternalPixelType;
+ typedef VectorImage<OutputInternalPixelType,2> OutputImageType;
+ typedef OutputImageType::PixelType OutputPixelType;
+
+ typedef TMaskImage MaskImageType;
+ typedef typename MaskImageType::Pointer MaskImagePointerType;
+ typedef typename MaskImageType::PixelType MaskPixelType;
- typedef otb::VectorRescaleIntensityImageFilter
- <InputImageType,OutputImageType> VectorRescalerType;
- typedef otb::MultiChannelExtractROI
- <OutputPixelType,OutputPixelType> ChannelExtractorType;
- typedef typename ChannelExtractorType::ChannelsType ChannelsType;
+ typedef VectorRescaleIntensityImageFilter<InputImageType,OutputImageType> VectorRescalerType;
+
+ typedef otb::MultiChannelExtractROI<InputInternalPixelType,
+ InputInternalPixelType> ChannelExtractorType;
+ typedef typename ChannelExtractorType::ChannelsType ChannelsType;
+ typedef Functor::MaskFunctor<InputPixelType,MaskPixelType,OutputPixelType> FunctorType;
+ typedef itk::BinaryFunctorImageFilter<OutputImageType, MaskImageType,
+ OutputImageType, FunctorType > FunctorFilterType;
+ typedef typename FunctorFilterType::Pointer FunctorFilterPointerType;
/** Method for creation through object factory */
itkNewMacro(Self);
-
+
/** Run-time type information */
itkTypeMacro(PrintableImageFilter,
itk::ImageToImageFilter);
-
+
/** Display */
void PrintSelf( std::ostream& os, itk::Indent indent ) const;
-
+
void SetChannel( unsigned int channel);
const ChannelsType GetChannels(void) const;
-
+
otbSetObjectMemberMacro(Rescaler,AutomaticInputMinMaxComputation,bool);
otbGetObjectMemberMacro(Rescaler,AutomaticInputMinMaxComputation,bool);
otbSetObjectMemberMacro(Rescaler,InputMinimum,InputPixelType);
otbGetObjectMemberMacro(Rescaler,InputMinimum,InputPixelType);
otbSetObjectMemberMacro(Rescaler,InputMaximum,InputPixelType);
otbGetObjectMemberMacro(Rescaler,InputMaximum,InputPixelType);
+
+
+ /**
+ * If set, only pixels within the mask will be classified.
+ * \param mask The input mask.
+ */
+ void SetInputMask(const MaskImageType * mask);
+
+ /**
+ * Get the input mask.
+ * \return The mask.
+ */
+ MaskImageType * GetInputMask(void);
+
+ itkSetMacro(UseMask, bool);
+ itkGetMacro(UseMask, bool);
+
+ ChannelsType const GetChannelList()
+ {
+ return m_ChannelList;
+ }
+ /* Set the selected channle index (order is important) */
+ void SetChannelList( ChannelsType chList )
+ {
+ if(chList.size() != 3)
+ {
+ itkExceptionMacro(<<"Invalid channel list, size is "<<chList.size()<<" instead of 3");
+ }
+ m_ChannelList = chList;
+ this->Modified();
+ }
+
+ /** Output Mask Object color. */
+ void SetObjectColor( OutputPixelType val )
+ {
+ if(val.GetSize() != 3)
+ {
+ itkExceptionMacro(<<"Invalid object color, size is "<<val.Size()<<" instead of 3");
+ }
+ m_ObjectColor = val;
+ m_MaskFilter->GetFunctor().SetObjectColor( val );
+ this->Modified();
+ }
+ itkGetMacro(ObjectColor, OutputPixelType);
+
+ void SetBackgroundMaskValue( MaskPixelType val )
+ {
+ m_BackgroundMaskValue = val;
+ m_MaskFilter->GetFunctor().SetBackgroundValue( val );
+ this->Modified();
+ }
+ itkGetMacro(BackgroundMaskValue, MaskPixelType);
-protected:
+ protected:
PrintableImageFilter();
+ void BeforeGenerateData();
void GenerateData();
-
-private:
-
+
+ private:
+
PrintableImageFilter(Self&); // intentionally not implemented
void operator=(const Self&); // intentionally not implemented
-
- typename VectorRescalerType::Pointer m_Rescaler;
+
+ typename VectorRescalerType::Pointer m_Rescaler;
typename ChannelExtractorType::Pointer m_Extractor;
+ // Foreground mask value
+ FunctorFilterPointerType m_MaskFilter;
+ // Objects (of the mask) will be displayer with the choosen color.
+ OutputPixelType m_ObjectColor;
+ // Use mask
+ bool m_UseMask;
+ // Used channel for output Image
+ ChannelsType m_ChannelList;
+ // Foreground mask value
+ //MaskPixelType m_ForegroundMaskValue;
+ // Background mask value
+ MaskPixelType m_BackgroundMaskValue;
+
};
-
+
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
diff --git a/Code/BasicFilters/otbPrintableImageFilter.txx b/Code/BasicFilters/otbPrintableImageFilter.txx
index d094287d5a6210194ef46ee79ad719f5bcfaafdf..71d758810ba2a3f901870b4ce959e408e8b9b66b 100644
--- a/Code/BasicFilters/otbPrintableImageFilter.txx
+++ b/Code/BasicFilters/otbPrintableImageFilter.txx
@@ -23,72 +23,157 @@
#include "otbPrintableImageFilter.h"
-
namespace otb
{
-template <class TInputImage>
-PrintableImageFilter<TInputImage>
+template <class TInputImage, class TMaskImage>
+PrintableImageFilter<TInputImage, TMaskImage>
::PrintableImageFilter()
{
m_Rescaler = VectorRescalerType::New();
m_Extractor = ChannelExtractorType::New();
+
+ m_Rescaler->SetInput( m_Extractor->GetOutput() );
+ m_MaskFilter = FunctorFilterType::New();
- m_Extractor->SetInput( m_Rescaler->GetOutput() );
+ m_ChannelList = ChannelsType(3,2);
+ m_ChannelList[1] = 3;
+ m_ChannelList[2] = 4;
+ m_ObjectColor.SetSize( 3 );
+ m_ObjectColor.Fill(255);
+ m_BackgroundMaskValue = 0;
}
-template <class TInputImage>
+template <class TInputImage, class TMaskImage>
void
-PrintableImageFilter<TInputImage>
+PrintableImageFilter<TInputImage, TMaskImage>
::SetChannel(unsigned int channel)
{
+ m_ChannelList.push_back(channel);
m_Extractor->SetChannel(channel);
this->Modified();
}
-template <class TInputImage>
-const typename PrintableImageFilter<TInputImage>::ChannelsType
-PrintableImageFilter<TInputImage>
+template <class TInputImage, class TMaskImage>
+const typename PrintableImageFilter<TInputImage,TMaskImage>::ChannelsType
+PrintableImageFilter<TInputImage, TMaskImage>
::GetChannels() const
{
- return m_Extractor->GetChannels();
+ return m_ChannelList;
}
-
-template <class TInputImage>
+template <class TInputImage, class TMaskImage>
void
-PrintableImageFilter<TInputImage>
-::GenerateData()
+PrintableImageFilter<TInputImage, TMaskImage>
+::SetInputMask(const MaskImageType * mask)
{
+ this->itk::ProcessObject::SetNthInput(1,const_cast<MaskImageType *>(mask));
+ m_UseMask = true;
+}
+
- if (m_Extractor->GetNbChannels() == 0)
+template <class TInputImage, class TMaskImage>
+typename PrintableImageFilter<TInputImage,TMaskImage>::MaskImageType *
+PrintableImageFilter<TInputImage, TMaskImage>
+::GetInputMask()
+{
+ if (this->GetNumberOfInputs()<2)
{
- m_Extractor->SetChannel(2);
- m_Extractor->SetChannel(3);
- m_Extractor->SetChannel(4);
+ return 0;
}
+ return static_cast<MaskImageType *>(this->itk::ProcessObject::GetInput(1));
+}
+
+
+template <class TInputImage, class TMaskImage>
+void
+PrintableImageFilter<TInputImage, TMaskImage>
+::BeforeGenerateData()
+{
+ if(m_Extractor->GetNbChannels() == 0 && m_ChannelList.size()>3)
+ {
+ itkExceptionMacro(<<"Invalid channel list (must be three channels instead of "<<m_ChannelList.size());
+ }
+ if(m_Extractor->GetNbChannels() == 0 && m_Extractor->GetNbChannels()>3)
+ {
+ itkExceptionMacro(<<"Invalid channel list (must be three channels instead of "<<m_Extractor->GetNbChannels());
+ }
+ if( m_ObjectColor.GetSize() != 3 )
+ {
+ itkExceptionMacro(<<"Invalid m_ObjectColor pixel size");
+ }
+
+ if( m_UseMask == true )
+ {
+ MaskImagePointerType inputMaskPtr = this->GetInputMask();
+ typename InputImageType::ConstPointer inputPtr = this->GetInput();
+
+ if( !inputMaskPtr )
+ {
+ itkExceptionMacro(<<"No mask detected");
+ }
+ if( inputMaskPtr->GetLargestPossibleRegion().GetSize() != inputPtr->GetLargestPossibleRegion().GetSize())
+ {
+ itkExceptionMacro(<<"Input size ("<<inputPtr->GetLargestPossibleRegion().GetSize()<<") and Mask size ("<<inputMaskPtr->GetLargestPossibleRegion().GetSize()<<") must be the same");
+ }
+ }
+}
+
+
+/*
+ * If no mask used, just rescale the input between 0 and 255.
+ * Else, sur M%askFunctor.
+ */
+template <class TInputImage, class TMaskImage>
+void
+PrintableImageFilter<TInputImage, TMaskImage>
+::GenerateData()
+{
+ this->BeforeGenerateData();
+
+ // let this loop to be compliant with previous version of the class where m_ChannelList didn't exist...
+ if(m_Extractor->GetNbChannels() == 0)
+ {
+ for(unsigned int i=0; i<m_ChannelList.size(); i++)
+ m_Extractor->SetChannel( m_ChannelList[i]);
+ }
+
+
+ m_Extractor->SetInput(this->GetInput());
+ m_Extractor->UpdateOutputInformation();
typename TInputImage::PixelType minimum,maximum;
- minimum.SetSize(this->GetInput()->GetNumberOfComponentsPerPixel());
- maximum.SetSize(this->GetInput()->GetNumberOfComponentsPerPixel());
+ minimum.SetSize(m_Extractor->GetNbChannels());
+ maximum.SetSize(m_Extractor->GetNbChannels());
minimum.Fill(0);
maximum.Fill(255);
- m_Rescaler->SetInput(this->GetInput());
m_Rescaler->SetOutputMinimum(minimum);
m_Rescaler->SetOutputMaximum(maximum);
m_Rescaler->SetClampThreshold(0.01);
- m_Extractor->GraftOutput( this->GetOutput() );
- m_Extractor->Update();
- this->GraftOutput( m_Extractor->GetOutput() );
+ if( m_UseMask == false )
+ {
+ m_Rescaler->GraftOutput( this->GetOutput() );
+ m_Rescaler->Update();
+ this->GraftOutput( m_Rescaler->GetOutput() );
+ }
+ else
+ {
+ m_MaskFilter->SetInput1(m_Rescaler->GetOutput());
+ m_MaskFilter->SetInput2(this->GetInputMask());
+ m_MaskFilter->GraftOutput( this->GetOutput() );
+ m_MaskFilter->Update();
+ this->GraftOutput( m_MaskFilter->GetOutput() );
+ }
+
}
-template <class TInputImage>
+template <class TInputImage, class TMaskImage>
void
-PrintableImageFilter<TInputImage>
+PrintableImageFilter<TInputImage, TMaskImage>
::PrintSelf( std::ostream& os, itk::Indent indent ) const
{
Superclass::PrintSelf(os,indent);
diff --git a/Code/BasicFilters/otbVectorRescaleIntensityImageFilter.txx b/Code/BasicFilters/otbVectorRescaleIntensityImageFilter.txx
index 4ddccd1c2c138b2a52f14f518a76d0bd8140e80b..b4cc883fad05b84ecd113299f57c640a1e029f09 100644
--- a/Code/BasicFilters/otbVectorRescaleIntensityImageFilter.txx
+++ b/Code/BasicFilters/otbVectorRescaleIntensityImageFilter.txx
@@ -76,6 +76,11 @@ void
VectorRescaleIntensityImageFilter<TInputImage, TOutputImage>
::BeforeThreadedGenerateData()
{
+ if(m_ClampThreshold==0.)
+ {
+ itkExceptionMacro(<<"Invalid Clamp Threshold must be greater than 0.0");
+ }
+
if (m_AutomaticInputMinMaxComputation)
{
diff --git a/Code/DisparityMap/otbBSplinesInterpolateTransformDeformationFieldGenerator.h b/Code/DisparityMap/otbBSplinesInterpolateTransformDeformationFieldGenerator.h
index 27b326a0c58cc38457cd5a930a4bbaf75b88cdcd..a47087fabcfaedd7ed919d09945b19b9afd75af6 100644
--- a/Code/DisparityMap/otbBSplinesInterpolateTransformDeformationFieldGenerator.h
+++ b/Code/DisparityMap/otbBSplinesInterpolateTransformDeformationFieldGenerator.h
@@ -19,7 +19,7 @@ PURPOSE. See the above copyright notices for more information.
#define __otbBSplinesInterpolateTransformDeformationFieldGenerator_h
#include "otbPointSetWithTransformToDeformationFieldGenerator.h"
-#include "ijBSplineScatteredDataPointSetToImageFilter.h"
+#include "itkBSplineScatteredDataPointSetToImageFilter.h"
#include "otbImage.h"
#include "otbMath.h"
#include <complex>
@@ -29,14 +29,14 @@ namespace otb
/** \class BSplinesInterpolateTransformDeformationFieldGenerator
* \brief This class generate the deformation field by using spline interpolation on the parameters of the transform.
*
- * Spline interpolation of non regularly scattered data is provided by a filter from the insight journal,
- * ij::BSplineScatteredDataPointSetToImageFilter. It allows interpolation using any spline order and implements a multi-level approach.
+ * Spline interpolation of non regularly scattered data is provided
+ * by the itk::BSplineScatteredDataPointSetToImageFilter. It allows interpolation using any spline order and implements a multi-level approach.
*
* This filter is used for each parameter. One can also specify the indices of the angular parameters. Angular parameters are first
* converted to complex exponential, the interpolated and converted back to the angular space. This is done to avoid interpolating angular discontinuities,
* which is a non-sense.
*
- * \sa ij::BSplineScatteredDataPointSetToImageFilter
+ * \sa itk::BSplineScatteredDataPointSetToImageFilter
*/
template <class TPointSet, class TDeformationField>
class ITK_EXPORT BSplinesInterpolateTransformDeformationFieldGenerator
@@ -72,7 +72,7 @@ public:
typedef itk::Vector<ValueType,2> PointSetDataType;
typedef otb::Image<PointSetDataType,DeformationFieldType::ImageDimension> InternalImageType;
typedef itk::PointSet<PointSetDataType,PointSetType::PointDimension> InternalPointSetType;
- typedef ij::BSplineScatteredDataPointSetToImageFilter<InternalPointSetType,InternalImageType>
+ typedef itk::BSplineScatteredDataPointSetToImageFilter<InternalPointSetType,InternalImageType>
SPlineInterpolateFilterType;
typedef typename SPlineInterpolateFilterType::Pointer SPlineInterpolateFilterPointerType;
diff --git a/Examples/IO/LidarToImageExample.cxx b/Examples/IO/LidarToImageExample.cxx
index 45466b20031ba780a020c6cb10ed67d1a6365e2d..067ff6c69272fbe6f9d52df72aa3a024ec5c0314 100644
--- a/Examples/IO/LidarToImageExample.cxx
+++ b/Examples/IO/LidarToImageExample.cxx
@@ -27,8 +27,8 @@
// Software Guide : BeginCommandLineArgs
// INPUTS: {TO_core_last_zoom.las}
-// OUTPUTS: {lidar-image-6.hdr}, {lidar-image-6-pretty.png}
-// 1.0 5 6
+// OUTPUTS: {lidar-image-8.hdr}, {lidar-image-8-pretty.png}
+// 1.0 5 8
// Software Guide : EndCommandLineArgs
@@ -261,8 +261,8 @@ int main( int argc, char* argv[] )
// \begin{figure}
// \center
// \includegraphics[width=0.40\textwidth]{lidar-image-4-pretty.eps}
- // \includegraphics[width=0.40\textwidth]{lidar-image-6-pretty.eps}
- // \itkcaption[Image from lidar data]{Image obtained with 4 level spline interpolation (left) and 6 levels (right)}
+ // \includegraphics[width=0.40\textwidth]{lidar-image-8-pretty.eps}
+ // \itkcaption[Image from lidar data]{Image obtained with 4 level spline interpolation (left) and 8 levels (right)}
// \label{fig:LIDARTOIMAGEEXAMPLE}
// \end{figure}
//
diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt
index 4fa8f0ae3a107e3649c810a18342528ff8d16898..e8075f3290f96bee6cf8ca1ccedad88488eada60 100644
--- a/Testing/Code/BasicFilters/CMakeLists.txt
+++ b/Testing/Code/BasicFilters/CMakeLists.txt
@@ -526,7 +526,16 @@ ADD_TEST(bfTvPrintableImageFilter ${BASICFILTERS_TESTS5}
${TEMP}/bfTvPrintableImageFilter.png
)
-
+ADD_TEST(bfTvPrintableImageFilterWithMask ${BASICFILTERS_TESTS5}
+ --compare-image ${NOTOL}
+ ${BASELINE}/bfTvPrintableImageFilterWithMask.png
+ ${TEMP}/bfTvPrintableImageFilterWithMask.png
+ otbPrintableImageFilterWithMask
+ ${INPUTDATA}/poupees.tif
+ ${INPUTDATA}/carre.png
+ ${TEMP}/bfTvPrintableImageFilterWithMask.png
+ 0 # background
+ )
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbBasicFiltersTests6 ~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1597,6 +1606,7 @@ otbUnaryImageFunctorWithVectorImageFilterNew.cxx
otbUnaryImageFunctorWithVectorImageFilter.cxx
otbPrintableImageFilterNew.cxx
otbPrintableImageFilter.cxx
+otbPrintableImageFilterWithMask.cxx
)
# A enrichir
diff --git a/Testing/Code/BasicFilters/otbBasicFiltersTests5.cxx b/Testing/Code/BasicFilters/otbBasicFiltersTests5.cxx
index d9d4a8d3d269e1f411c785af676478c539d27dde..31fdd53aa3c75a8e084a193770424ce3ff365b9f 100644
--- a/Testing/Code/BasicFilters/otbBasicFiltersTests5.cxx
+++ b/Testing/Code/BasicFilters/otbBasicFiltersTests5.cxx
@@ -27,16 +27,17 @@
void RegisterTests()
{
- REGISTER_TEST(otbStreamingResampleImageFilterNew);
- REGISTER_TEST(otbStreamingResampleImageFilter);
- REGISTER_TEST(otbStreamingResampleImageFilterCompareWithITK);
- REGISTER_TEST(otbResampleSLCImage);
- REGISTER_TEST(otbStreamingStatisticsVectorImageFilterNew);
- REGISTER_TEST(otbStreamingStatisticsVectorImageFilter);
- REGISTER_TEST(otbMatrixTransposeMatrixImageFilterNew);
- REGISTER_TEST(otbMatrixTransposeMatrixImageFilter);
- REGISTER_TEST(otbUnaryImageFunctorWithVectorImageFilterNew);
- REGISTER_TEST(otbUnaryImageFunctorWithVectorImageFilter);
- REGISTER_TEST(otbPrintableImageFilterNew);
- REGISTER_TEST(otbPrintableImageFilter);
+REGISTER_TEST(otbStreamingResampleImageFilterNew);
+REGISTER_TEST(otbStreamingResampleImageFilter);
+REGISTER_TEST(otbStreamingResampleImageFilterCompareWithITK);
+REGISTER_TEST(otbResampleSLCImage);
+REGISTER_TEST(otbStreamingStatisticsVectorImageFilterNew);
+REGISTER_TEST(otbStreamingStatisticsVectorImageFilter);
+REGISTER_TEST(otbMatrixTransposeMatrixImageFilterNew);
+REGISTER_TEST(otbMatrixTransposeMatrixImageFilter);
+REGISTER_TEST(otbUnaryImageFunctorWithVectorImageFilterNew);
+REGISTER_TEST(otbUnaryImageFunctorWithVectorImageFilter);
+REGISTER_TEST(otbPrintableImageFilterNew);
+REGISTER_TEST(otbPrintableImageFilter);
+REGISTER_TEST(otbPrintableImageFilterWithMask);
}
diff --git a/Testing/Code/BasicFilters/otbPrintableImageFilterWithMask.cxx b/Testing/Code/BasicFilters/otbPrintableImageFilterWithMask.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..132a85311233e6efefed390177c93f07f79019d7
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbPrintableImageFilterWithMask.cxx
@@ -0,0 +1,87 @@
+/*=========================================================================
+
+ 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.
+
+=========================================================================*/
+#if defined(_MSC_VER)
+#pragma warning ( disable : 4786 )
+#endif
+
+#include "otbImageFileReader.h"
+#include "otbStreamingImageFileWriter.h"
+#include "otbVectorImage.h"
+
+#include "otbPrintableImageFilter.h"
+
+int otbPrintableImageFilterWithMask( int argc, char * argv[] )
+{
+ const char * inputFilename = argv[1];
+ const char * masktFilename = argv[2];
+ const char * outputFilename = argv[3];
+
+
+ typedef double InputPixelType;
+ const unsigned int Dimension = 2;
+
+ typedef otb::VectorImage< InputPixelType, Dimension > InputImageType;
+ typedef otb::PrintableImageFilter< InputImageType> FilterType;
+ typedef FilterType::OutputImageType OutputImageType;
+ typedef OutputImageType::PixelType OutputPixelType;
+ typedef FilterType::MaskImageType MaskImageType;
+ typedef FilterType::MaskPixelType MaskPixelType;
+
+ typedef otb::ImageFileReader< InputImageType > InputReaderType;
+ typedef otb::ImageFileReader< MaskImageType > MaskReaderType;
+ typedef otb::StreamingImageFileWriter< OutputImageType > WriterType;
+
+
+ FilterType::Pointer printableImageFilter = FilterType::New();
+ InputReaderType::Pointer inputReader = InputReaderType::New();
+ MaskReaderType::Pointer maskReader = MaskReaderType::New();
+ WriterType::Pointer writer = WriterType::New();
+
+ inputReader->SetFileName( inputFilename );
+ maskReader->SetFileName( masktFilename );
+ maskReader->GenerateOutputInformation();
+ writer->SetFileName( outputFilename );
+
+ printableImageFilter->SetInput( inputReader->GetOutput() );
+ printableImageFilter->SetInputMask( maskReader->GetOutput()/*rescaler->GetOutput()*/ );
+
+
+ FilterType::ChannelsType chList;
+ chList.push_back(3);
+ chList.push_back(2);
+ chList.push_back(1);
+ printableImageFilter->SetChannelList(chList);
+
+ OutputPixelType objectColor;
+ objectColor.SetSize(3);
+ objectColor[0] = 0;
+ objectColor[1] = 255;
+ objectColor[2] = 0;
+ printableImageFilter->SetObjectColor( objectColor );
+ printableImageFilter->SetBackgroundMaskValue(static_cast<MaskPixelType>(atof(argv[4])));
+
+
+ writer->SetInput( printableImageFilter->GetOutput() );
+
+ writer->Update();
+
+
+ return EXIT_SUCCESS;
+}
+
+
diff --git a/Testing/Code/MultiScale/otbWaveletFilterBank.cxx b/Testing/Code/MultiScale/otbWaveletFilterBank.cxx
index d35b3aaac1a5d2941062df3ff158707977beb1c9..72bed2e0ef7c8da6b5ab9fe09192c68e7b47f82e 100644
--- a/Testing/Code/MultiScale/otbWaveletFilterBank.cxx
+++ b/Testing/Code/MultiScale/otbWaveletFilterBank.cxx
@@ -24,7 +24,6 @@
#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
-#include "itkRescaleIntensityImageFilter.h"
#include "otbWaveletOperator.h"
#include "otbWaveletFilterBank.h"
@@ -72,18 +71,10 @@ int otbWaveletFilterBank( int argc, char * argv[] )
invFilter->SetSubsampleImageFactor( filter->GetSubsampleImageFactor() );
/* Saving output */
- typedef unsigned char OutputPixelType;
- typedef otb::Image< OutputPixelType, Dimension > OutputImageType;
- typedef itk::RescaleIntensityImageFilter< ImageType, OutputImageType > RescalerType;
- RescalerType::Pointer rescaler = RescalerType::New();
- rescaler->SetOutputMinimum(0);
- rescaler->SetOutputMaximum(255);
- rescaler->SetInput( invFilter->GetOutput() );
-
- typedef otb::ImageFileWriter< OutputImageType > WriterType;
+ typedef otb::ImageFileWriter< ImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( outputFileName );
- writer->SetInput( rescaler->GetOutput() );
+ writer->SetInput( invFilter->GetOutput() );
writer->Update();
return EXIT_SUCCESS;
diff --git a/Testing/Code/MultiScale/otbWaveletPacketTransform.cxx b/Testing/Code/MultiScale/otbWaveletPacketTransform.cxx
index 67306a3f5e1e21c9f3842ce34a738593b5ed5ffb..8e50f7dbba7b0721ca2681c5c18649a535580b9f 100644
--- a/Testing/Code/MultiScale/otbWaveletPacketTransform.cxx
+++ b/Testing/Code/MultiScale/otbWaveletPacketTransform.cxx
@@ -24,7 +24,6 @@
#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
-#include "itkRescaleIntensityImageFilter.h"
#include "otbWaveletOperator.h"
#include "otbWaveletFilterBank.h"
@@ -83,18 +82,10 @@ int otbWaveletPacketTransform( int argc, char * argv[] )
invFilter->Update();
/* Writing the output */
- typedef unsigned char OutputPixelType;
- typedef otb::Image< OutputPixelType, Dimension > OutputImageType;
- typedef itk::RescaleIntensityImageFilter< ImageType, OutputImageType > RescalerType;
- RescalerType::Pointer rescaler = RescalerType::New();
- rescaler->SetOutputMinimum(0);
- rescaler->SetOutputMaximum(255);
- rescaler->SetInput( invFilter->GetOutput() );
-
- typedef otb::ImageFileWriter< OutputImageType > WriterType;
+ typedef otb::ImageFileWriter< ImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( outputFileName );
- writer->SetInput( rescaler->GetOutput() );
+ writer->SetInput( invFilter->GetOutput() );
writer->Update();
return EXIT_SUCCESS;
diff --git a/Testing/Code/MultiScale/otbWaveletTransform.cxx b/Testing/Code/MultiScale/otbWaveletTransform.cxx
index 118ec0a336e4216e4c97ae3710a088bc2c8c7fc1..840972003a39b3b2f68ebec202ac893130cbf9ce 100644
--- a/Testing/Code/MultiScale/otbWaveletTransform.cxx
+++ b/Testing/Code/MultiScale/otbWaveletTransform.cxx
@@ -77,18 +77,10 @@ int otbWaveletTransform( int argc, char * argv[] )
invFilter->Update();
/* Writing the output */
- typedef unsigned char OutputPixelType;
- typedef otb::Image< OutputPixelType, Dimension > OutputImageType;
- typedef itk::RescaleIntensityImageFilter< ImageType, OutputImageType > RescalerType;
- RescalerType::Pointer rescaler = RescalerType::New();
- rescaler->SetOutputMinimum(0);
- rescaler->SetOutputMaximum(255);
- rescaler->SetInput( invFilter->GetOutput() );
-
- typedef otb::ImageFileWriter<OutputImageType> WriterType;
+ typedef otb::ImageFileWriter<ImageType> WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( outputFileName );
- writer->SetInput( rescaler->GetOutput() );
+ writer->SetInput( invFilter->GetOutput() );
writer->Update();
return EXIT_SUCCESS;
diff --git a/Testing/Utilities/CMakeLists.txt b/Testing/Utilities/CMakeLists.txt
index 1bf939b29858c5446b03979d520ed8eb9f901ef7..99c9b3b7b0c11eee33c341e31f726acb7fdbbeda 100644
--- a/Testing/Utilities/CMakeLists.txt
+++ b/Testing/Utilities/CMakeLists.txt
@@ -51,16 +51,6 @@ ADD_TEST(utTuOssimElevManagerTest ${UTILITIES_TESTS}
ossimElevManagerTest
)
-
-ADD_TEST(utTvIjBSplineScatteredDataPointSetToImageFilterTest ${UTILITIES_TESTS}
---compare-image ${EPSILON}
- ${BASELINE}/ijBSplineScatteredDataPointSetToImageFilterTestOutput.hdr
- ${TEMP}/ijBSplineScatteredDataPointSetToImageFilterTestOutput.hdr
- ijBSplineScatteredDataPointSetToImageFilterTest
- ${INPUTDATA}/brain.png
- ${TEMP}/ijBSplineScatteredDataPointSetToImageFilterTestOutput.hdr
-)
-
# ------- lib otbsvm ------------------------------
ADD_TEST(utTuSvmKernelFunctorTest ${UTILITIES_TESTS}
@@ -442,7 +432,6 @@ SET(UtilitiesTests_SRCS
ossimIntegrationTest.cxx
ossimKeywordlistTest.cxx
ossimElevManagerTest.cxx
-ijBSplineScatteredDataPointSetToImageFilterTest.cxx
svmGenericKernelFunctor.cxx
svmTest.cxx
svmGenericKernelTest.cxx
@@ -492,7 +481,6 @@ ENDIF(NOT OTB_USE_EXTERNAL_EXPAT)
# ------- Select sources files suppress warning -----------------------------------
SET(UtilitiesTests_DisableWarning_SRCS
- ijBSplineScatteredDataPointSetToImageFilterTest.cxx
siftfast.cpp
openthreadsWorkCrew.cpp
expatchardata.cxx
diff --git a/Testing/Utilities/ijBSplineScatteredDataPointSetToImageFilterTest.cxx b/Testing/Utilities/ijBSplineScatteredDataPointSetToImageFilterTest.cxx
deleted file mode 100644
index af15484cf5a300af4b9762dbb459079e7f5ea7f9..0000000000000000000000000000000000000000
--- a/Testing/Utilities/ijBSplineScatteredDataPointSetToImageFilterTest.cxx
+++ /dev/null
@@ -1,122 +0,0 @@
-#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 double 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/otbUtilitiesTests.cxx b/Testing/Utilities/otbUtilitiesTests.cxx
index bdb10195cad7de53e0a938a488db37ebe4c2e6bc..fc88d8c0bc50510d68071c847145e1aa691749c9 100644
--- a/Testing/Utilities/otbUtilitiesTests.cxx
+++ b/Testing/Utilities/otbUtilitiesTests.cxx
@@ -29,7 +29,6 @@ void RegisterTests()
REGISTER_TEST(ossimIntegrationTest);
REGISTER_TEST(ossimKeywordlistTest);
REGISTER_TEST(ossimElevManagerTest);
-REGISTER_TEST(ijBSplineScatteredDataPointSetToImageFilterTest);
REGISTER_TEST(svmGenericKernelFunctor);
REGISTER_TEST(svmTest);
REGISTER_TEST(svmGenericKernelTest);
diff --git a/Utilities/ITK/Code/BasicFilters/itkBSplineDecompositionImageFilter.h b/Utilities/ITK/Code/BasicFilters/itkBSplineDecompositionImageFilter.h
index 1a86df8bca9ef255da050191806cd2af97f08894..970426d417c6e4f7c0b85826d1107b68dbe27915 100644
--- a/Utilities/ITK/Code/BasicFilters/itkBSplineDecompositionImageFilter.h
+++ b/Utilities/ITK/Code/BasicFilters/itkBSplineDecompositionImageFilter.h
@@ -12,8 +12,8 @@
Portions of this code are covered under the VTK copyright.
See VTKCopyright.txt or http://www.kitware.com/VTKCopyright.htm for details.
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ 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.
=========================================================================*/
@@ -34,7 +34,7 @@ namespace itk
* \brief Calculates the B-Spline coefficients of an image. Spline order may be from 0 to 5.
*
* This class defines N-Dimension B-Spline transformation.
- * It is based on:
+ * It is based on:
* [1] M. Unser,
* "Splines: A Perfect Fit for Signal and Image Processing,"
* IEEE Signal Processing Magazine, vol. 16, no. 6, pp. 22-38,
@@ -48,7 +48,7 @@ namespace itk
* IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 834-848,
* February 1993.
* And code obtained from bigwww.epfl.ch by Philippe Thevenaz
- *
+ *
* Limitations: Spline order must be between 0 and 5.
* Spline order must be set before setting the image.
* Uses mirror boundary conditions.
@@ -60,10 +60,10 @@ namespace itk
* ***TODO: Is this an ImageFilter? or does it belong to another group?
* \ingroup ImageFilters
* \ingroup SingleThreaded
- * \ingroup CannotBeStreamed
+ * \ingroup CannotBeStreamed
*/
template <class TInputImage, class TOutputImage>
-class ITK_EXPORT BSplineDecompositionImageFilter :
+class ITK_EXPORT BSplineDecompositionImageFilter :
public ImageToImageFilter<TInputImage,TOutputImage>
{
public:
@@ -75,7 +75,7 @@ public:
/** Run-time type information (and related methods). */
itkTypeMacro(BSplineDecompositionImageFilter, ImageToImageFilter);
-
+
/** New macro for creation of through a Smart Pointer */
itkNewMacro( Self );
@@ -102,13 +102,14 @@ public:
/** Begin concept checking */
itkConceptMacro(DimensionCheck,
(Concept::SameDimension<ImageDimension, OutputImageDimension>));
- itkConceptMacro(InputConvertibleToDoubleCheck,
- (Concept::Convertible<typename TInputImage::PixelType, double>));
- itkConceptMacro(OutputConvertibleToDoubleCheck,
- (Concept::Convertible<typename TOutputImage::PixelType, double>));
- itkConceptMacro(InputConvertibleToOutputCheck,
- (Concept::Convertible<typename TInputImage::PixelType,
- typename TOutputImage::PixelType>));
+// itkConceptMacro(InputConvertibleToDoubleCheck,
+// (Concept::Convertible<typename TInputImage::PixelType, double>));
+// itkConceptMacro(OutputConvertibleToDoubleCheck,
+// (Concept::Convertible<typename TOutputImage::PixelType, double>));
+// itkConceptMacro(InputConvertibleToOutputCheck,
+// (Concept::Convertible<typename TInputImage::PixelType,
+// typename TOutputImage::PixelType>));
+//FIXME probably need a fix in the OptBSplineInterpolateImageFilter
itkConceptMacro(DoubleConvertibleToOutputCheck,
(Concept::Convertible<double, typename TOutputImage::PixelType>));
/** End concept checking */
@@ -125,7 +126,7 @@ protected:
void GenerateInputRequestedRegion();
/** This filter must produce all of its output at once. */
- void EnlargeOutputRequestedRegion( DataObject *output );
+ void EnlargeOutputRequestedRegion( DataObject *output );
/** These are needed by the smoothing spline routine. */
std::vector<double> m_Scratch; // temp storage for processing of Coefficients
@@ -166,7 +167,7 @@ private:
/** Copies a vector of data from m_Scratch to the Coefficients image. */
void CopyScratchToCoefficients( OutputLinearIterator & );
-
+
};
diff --git a/Utilities/ITK/Code/Review/itkBSplineScatteredDataPointSetToImageFilter.txx b/Utilities/ITK/Code/Review/itkBSplineScatteredDataPointSetToImageFilter.txx
index 9ab890b7455e39f320a98478386023c0a885f441..f73d49fb2bf3beb1a4297c6646a882bab777e679 100644
--- a/Utilities/ITK/Code/Review/itkBSplineScatteredDataPointSetToImageFilter.txx
+++ b/Utilities/ITK/Code/Review/itkBSplineScatteredDataPointSetToImageFilter.txx
@@ -659,14 +659,20 @@ BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
}
RealType wc = this->m_PointWeights->GetElement( It.Index() );
RealType t = Itw.Get();
- omega->SetPixel( idx, omega->GetPixel( idx ) + wc*t*t );
+ if(idx[0]>=0 && idx[1]>=0 && idx[0]< omega->GetLargestPossibleRegion().GetSize()[0] && idx[1]< omega->GetLargestPossibleRegion().GetSize()[1])
+ omega->SetPixel( idx, omega->GetPixel( idx ) + wc*t*t );
+
PointDataType data = It.Value();
data *= ( t / w2_sum );
Itp.Set( data );
data *= ( t * t * wc );
+
+ if(idx[0]>=0 && idx[1]>=0 && idx[0]< delta->GetLargestPossibleRegion().GetSize()[0] && idx[1]< delta->GetLargestPossibleRegion().GetSize()[1])
+ {
delta->SetPixel( idx, delta->GetPixel( idx ) + data );
delta->GetPixel( idx ) + data;
+ }
}
}
@@ -913,7 +919,9 @@ BSplineScatteredDataPointSetToImageFilter<TInputPointSet, TOutputImage>
idx[dimension] %=
lattice->GetLargestPossibleRegion().GetSize()[dimension];
}
- data += ( lattice->GetPixel( idx ) * B );
+
+ if(idx[dimension]>=0 && idx[dimension]<lattice->GetLargestPossibleRegion().GetSize()[dimension])
+ data += ( lattice->GetPixel( idx ) * B );
}
It.Set( data );
}
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