diff --git a/Code/Common/otbSubsampledImageRegionConstIterator.h b/Code/Common/otbSubsampledImageRegionConstIterator.h
new file mode 100644
index 0000000000000000000000000000000000000000..9bb299a360b2f9fe7103cc8f6e192a287cee0d06
--- /dev/null
+++ b/Code/Common/otbSubsampledImageRegionConstIterator.h
@@ -0,0 +1,350 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+ Copyright (c) Institut Telecom / Telecom Bretagne. All rights reserved.
+ See ITCopyright.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.
+
+=========================================================================*/
+
+#ifndef __otbSubsampledImageRegionConstIterator_h
+#define __otbSubsampledImageRegionConstIterator_h
+
+#include "itkImageRegionConstIterator.h"
+
+namespace otb {
+
+/** SubsampledImageRegionConstIterator
+ * \brief Regular subsample iterator over an image
+ *
+ * This iterator is a itk::ImageRegionConstIterator that perform a subsampled
+ * scan over an image. It runs one pixel over X (in row, line, slice... dimensions),
+ * if X is the (integer) value of the SubsampleFactor.
+ *
+ * \ingroup ImageIterator
+ * \sa StationaryFilterBank
+ */
+template < class TImage >
+class ITK_EXPORT SubsampledImageRegionConstIterator
+ : public itk::ImageRegionConstIterator< TImage >
+{
+public:
+ /** Standard typedef. */
+ typedef SubsampledImageRegionConstIterator Self;
+ typedef itk::ImageRegionConstIterator<TImage> Superclass;
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(SubsampledImageRegionConstIterator, ImageRegionConstIterator);
+
+ /** Dimension of the image the iterator walks. This constant is needed so
+ * functions that are templated over image iterator type (as opposed to
+ * being templated over pixel type and dimension) can have compile time
+ * access to the dimension of the image that the iterator walks. */
+ itkStaticConstMacro(ImageIteratorDimension, unsigned int,
+ Superclass::ImageIteratorDimension);
+
+ /** Index typedef support. While this was already typdef'ed in the superclass
+ * it needs to be redone here for this subclass to compile properly with gcc. */
+ typedef typename Superclass::IndexType IndexType;
+
+ /** Size typedef support. While this was already typdef'ed in the superclass
+ * it needs to be redone here for this subclass to compile properly with gcc. */
+ typedef typename Superclass::SizeType SizeType;
+
+ /** Region typedef support. */
+ typedef typename Superclass::RegionType RegionType;
+
+
+ /** Image typedef support. While this was already typdef'ed in the superclass
+ * it needs to be redone here for this subclass to compile properly with gcc. */
+ typedef typename Superclass::ImageType ImageType;
+
+ /** PixelContainer typedef support. Used to refer to the container for
+ * the pixel data. While this was already typdef'ed in the superclass
+ * it needs to be redone here for this subclass to compile properly with gcc. */
+ typedef typename Superclass::PixelContainer PixelContainer;
+ typedef typename Superclass::PixelContainerPointer PixelContainerPointer;
+
+ /** Internal Pixel Type */
+ typedef typename Superclass::InternalPixelType InternalPixelType;
+
+ /** External Pixel Type */
+ typedef typename Superclass::PixelType PixelType;
+
+ /** Accessor type that convert data between internal and external
+ * representations. */
+ typedef typename Superclass::AccessorType AccessorType;
+
+ /** Index value used for pixel location */
+ typedef typename Superclass::IndexValueType IndexValueType;
+
+ /** Default constructor. Needed since we provide a cast constructor. */
+ SubsampledImageRegionConstIterator() : itk::ImageRegionConstIterator<TImage> ()
+ {
+ m_SubsampleFactor = 1;
+ m_SubSampledEndOffset = this->m_EndOffset;
+
+ const IndexType& startIndex = this->m_Region.GetIndex();
+ const SizeType& size = this->m_Region.GetSize();
+
+ for ( unsigned int i = 0; i < ImageIteratorDimension; i++ )
+ {
+ m_LastUsableIndex[i] = startIndex[i]
+ + static_cast<IndexValueType>( m_SubsampleFactor * ((size[i]-1) / m_SubsampleFactor) );
+ }
+ }
+
+ /** Constructor establishes an iterator to walk a particular image and a
+ * particular region of that image. */
+ SubsampledImageRegionConstIterator ( const ImageType *ptr,
+ const RegionType ®ion)
+ : itk::ImageRegionConstIterator< TImage > ( ptr, region )
+ {
+ m_SubsampleFactor = 1;
+ m_SubSampledEndOffset = this->m_EndOffset;
+
+ const IndexType& startIndex = this->m_Region.GetIndex();
+ const SizeType& size = this->m_Region.GetSize();
+
+ for ( unsigned int i = 0; i < ImageIteratorDimension; i++ )
+ {
+ m_LastUsableIndex[i] = startIndex[i]
+ + static_cast<IndexValueType>( m_SubsampleFactor * ((size[i]-1) / m_SubsampleFactor) );
+ }
+ }
+
+ /** Constructor that can be used to cast from an ImageIterator to an
+ * SubsampledImageRegionConstIterator. Many routines return an ImageIterator
+ * but for a particular task, you may want an SubsampledImageRegionConstIterator.
+ * Rather than provide overloaded APIs that return different types of Iterators, itk
+ * returns ImageIterators and uses constructors to cast from an
+ * ImageIterator to a SubsampledImageRegionConstIterator. */
+ SubsampledImageRegionConstIterator( const itk::ImageIterator<TImage> &it )
+ : itk::ImageRegionConstIterator< TImage >( it )
+ {
+ m_SubsampleFactor = 1;
+ m_SubSampledEndOffset = this->m_EndOffset;
+
+ const IndexType& startIndex = this->m_Region.GetIndex();
+ const SizeType& size = this->m_Region.GetSize();
+
+ for ( unsigned int i = 0; i < ImageIteratorDimension; i++ )
+ {
+ m_LastUsableIndex[i] = startIndex[i]
+ + static_cast<IndexValueType>( m_SubsampleFactor * ((size[i]-1) / m_SubsampleFactor) );
+ }
+ }
+
+ /** Constructor that can be used to cast from an ImageConstIterator to an
+ * SubsampledImageRegionConstIterator. Many routines return an ImageIterator
+ * but for a particular task, you may want an SubsampledImageRegionConstIterator.
+ * Rather than provide overloaded APIs that return different types of Iterators, itk
+ * returns ImageIterators and uses constructors to cast from an
+ * ImageIterator to a SubsampledImageRegionConstIterator. */
+ SubsampledImageRegionConstIterator( const itk::ImageConstIterator<TImage> &it)
+ : itk::ImageRegionConstIterator< TImage >( it )
+ {
+ m_SubsampleFactor = 1;
+
+ const IndexType& startIndex = this->m_Region.GetIndex();
+ const SizeType& size = this->m_Region.GetSize();
+
+ for ( unsigned int i = 0; i < ImageIteratorDimension; i++ )
+ {
+ m_LastUsableIndex[i] = startIndex[i]
+ + static_cast<IndexValueType>( m_SubsampleFactor * ((size[i]-1) / m_SubsampleFactor) );
+ }
+
+ }
+
+ /** Set / Get the subsample factor */
+ void SetSubsampleFactor ( unsigned int factor )
+ {
+ this->m_SubsampleFactor = factor;
+
+ // Evaluate the last possible pixel.
+ const IndexType& startIndex = this->m_Region.GetIndex();
+ const SizeType& size = this->m_Region.GetSize();
+
+ for ( unsigned int i = 0; i < ImageIteratorDimension; i++ )
+ {
+ m_LastUsableIndex[i] = startIndex[i]
+ + static_cast<IndexValueType>( m_SubsampleFactor * ( (size[i]-1) / m_SubsampleFactor ) );
+ }
+
+ m_SubSampledEndOffset = this->m_Image->ComputeOffset( m_LastUsableIndex ) + 1;
+
+ otbGenericMsgDebugMacro(<<"m_EndOffset=" << this->m_EndOffset);
+ otbGenericMsgDebugMacro(<<"m_SubSampledEndOffset="<< m_SubSampledEndOffset);
+
+ for ( unsigned int i = 0; i < ImageIteratorDimension; i++ )
+ {
+ otbGenericMsgDebugMacro(<<"m_LastUsableIndex[" << i << "] = " << m_LastUsableIndex[i]);
+ }
+
+}
+
+ unsigned int GetSubsampleFactor () const
+ {
+ return this->m_SubsampleFactor;
+ }
+
+ /** Move an iterator to the beginning of the region. "Begin" is
+ * defined as the first pixel in the region. */
+ void GoToBegin()
+ {
+ Superclass::GoToBegin();
+
+ const SizeType& size = this->m_Region.GetSize();
+
+ this->m_SpanBeginOffset = this->m_Offset;
+ this->m_SpanEndOffset = this->m_Offset
+ + static_cast<IndexValueType>( m_SubsampleFactor * ((size[0]-1) / m_SubsampleFactor) )
+ + 1;
+ }
+
+ /** Move an iterator to the end of the region. "End" is defined as
+ * one pixel past the last pixel of the region. */
+ void GoToEnd()
+ {
+ Superclass::GoToEnd();
+ this->m_Offset = m_SubSampledEndOffset - 1;
+ this->m_SpanEndOffset = this->m_Offset + 1;
+ this->m_SpanBeginOffset = this->m_Offset - m_LastUsableIndex[0];
+ }
+
+ /** Is the iterator at the beginning of the region?
+ * "Begin" is defined here as before the first pixel
+ * of the region.
+ * On the contrary to Superclasses, now you cannot get
+ * *it even if IsAtBegin() is true!!!
+ */
+ bool IsAtBegin(void) const
+ {
+ return ( this->m_Offset <= this->m_BeginOffset );
+ }
+
+ /** Is the iterator at the end of the region?
+ * "End" is defined as past the last candidate pixel
+ * of the region (under the subsample point of view).
+ *
+ * For instance, if in 1D, the Size is 10 with a subsample
+ * factor of 4, the last candidate offset would be 8.
+ */
+ bool IsAtEnd(void) const
+ {
+ return ( this->m_Offset >= m_SubSampledEndOffset );
+ }
+
+ /** Set the index. No bounds checking is performed. This is overridden
+ * from the parent because we have an extra ivar.
+ * \sa GetIndex */
+ void SetIndex(const IndexType &ind)
+ {
+ Superclass::SetIndex( ind );
+ }
+
+ /** Increment (prefix) the fastest moving dimension of the iterator's index.
+ * This operator will constrain the iterator within the region (i.e. the
+ * iterator will automatically wrap from the end of the row of the region
+ * to the beginning of the next row of the region) up until the iterator
+ * tries to moves past the last pixel of the region. Here, the iterator
+ * will be set to be one pixel past the end of the region.
+ * \sa operator++(int) */
+ Self & operator++()
+ {
+ // On the contrary to itk::ImageRegionConstIterator, m_Offset to
+ // not incremented before the test
+ if( this->m_Offset + m_SubsampleFactor >= this->m_SpanEndOffset )
+ {
+ this->Increment();
+ }
+ else
+ {
+ // Now, the increment is performed
+ this->m_Offset += m_SubsampleFactor;
+ }
+ return *this;
+ }
+
+ /** Decrement (prefix) the fastest moving dimension of the iterator's index.
+ * This operator will constrain the iterator within the region (i.e. the
+ * iterator will automatically wrap from the beginning of the row of the region
+ * to the end of the next row of the region) up until the iterator
+ * tries to moves past the first pixel of the region. Here, the iterator
+ * will be set to be one pixel past the beginning of the region.
+ * \sa operator--(int) */
+ Self & operator--()
+ {
+ // On the contrary to itk::ImageRegionConstIterator, m_Offset
+ // is not decremented here (it may become negative!)
+ if ( this->m_Offset < this->m_SpanBeginOffset + m_SubsampleFactor )
+ {
+ this->Decrement();
+ }
+ else
+ {
+ // Now we can
+ this->m_Offset -= m_SubsampleFactor;
+ }
+ return *this;
+ }
+
+ /** In order to help copy into a new Image, give the new region parameters
+ */
+ RegionType GetNewRegion () const
+ {
+ IndexType startIndex = this->m_Region.GetIndex();
+ SizeType size = this->m_Region.GetSize();
+
+ for ( unsigned int i = 0; i < ImageIteratorDimension; i++ )
+ {
+ startIndex[i] /= m_SubsampleFactor;
+ --size[i];
+ size[i] /= m_SubsampleFactor;
+ ++size[i];
+ }
+
+ RegionType newRegion;
+ newRegion.SetIndex( startIndex );
+ newRegion.SetSize( size );
+
+ for ( unsigned int i = 0; i < ImageIteratorDimension; i++ )
+ {
+ otbGenericMsgDebugMacro(<<"NewRegionIndex[" << i << "] = " << startIndex[i]);
+ otbGenericMsgDebugMacro(<<"NewRegionSize [" << i << "] = " << size[i]);
+ }
+
+ return newRegion;
+ }
+
+protected:
+ unsigned int m_SubsampleFactor;
+ unsigned long m_SubSampledEndOffset;
+ IndexType m_LastUsableIndex;
+
+private:
+ void Increment(); // jump in a direction other than the fastest moving
+ void Decrement(); // go back in a direction other than the fastest moving
+};
+
+
+} // end of namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSubsampledImageRegionConstIterator.txx"
+#endif
+
+#endif
+
diff --git a/Code/Common/otbSubsampledImageRegionConstIterator.txx b/Code/Common/otbSubsampledImageRegionConstIterator.txx
new file mode 100644
index 0000000000000000000000000000000000000000..022e25c3da378cbea0310ff4a9b3ef2adeb32a52
--- /dev/null
+++ b/Code/Common/otbSubsampledImageRegionConstIterator.txx
@@ -0,0 +1,116 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+ Copyright (c) Institut Telecom / Telecom Bretagne. All rights reserved.
+ See ITCopyright.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.
+
+=========================================================================*/
+
+#ifndef __otbSubsampledImageRegionConstIterator_txx
+#define __otbSubsampledImageRegionConstIterator_txx
+
+#include "otbSubsampledImageRegionConstIterator.h"
+
+namespace otb {
+
+template < class TImage >
+void
+SubsampledImageRegionConstIterator< TImage >
+::Increment ()
+{
+ // Get the index of the last pixel on the span (row)
+ IndexType ind = this->m_Image->ComputeIndex( static_cast<typename Superclass::OffsetValueType>(this->m_Offset) );
+
+ const IndexType& startIndex = this->m_Region.GetIndex();
+ const SizeType& size = this->m_Region.GetSize();
+
+ // Increment along a row, then wrap at the end of the region row.
+ unsigned int dim;
+
+ // Check to see if we are past the last pixel in the region
+ // Note that ++ind[0] moves to the next pixel along the row.
+ ind[0] += static_cast<IndexValueType>(m_SubsampleFactor);
+ bool done = (ind[0] > m_LastUsableIndex[0] );
+ for (unsigned int i=1; done && i < ImageIteratorDimension; i++)
+ {
+ done = (ind[i] >= m_LastUsableIndex[i]);
+ }
+
+ // if the iterator is outside the region (but not past region end) then
+ // we need to wrap around the region
+ dim = 0;
+ if (!done)
+ {
+ while ( ( dim+1 < ImageIteratorDimension )
+ && (ind[dim] > m_LastUsableIndex[dim]) )
+ {
+ ind[dim] = startIndex[dim];
+ ind[++dim] += static_cast<IndexValueType>(m_SubsampleFactor);
+ }
+ }
+ this->m_Offset = this->m_Image->ComputeOffset( ind );
+ this->m_SpanEndOffset = this->m_Offset
+ + static_cast<IndexValueType>( m_SubsampleFactor * ((size[0]-1) / m_SubsampleFactor) )
+ + 1;
+ this->m_SpanBeginOffset = this->m_Offset;
+}
+
+template < class TImage >
+void
+SubsampledImageRegionConstIterator< TImage >
+::Decrement ()
+{
+ // Get the index of the first pixel on the span (row)
+ typename itk::ImageIterator<TImage>::IndexType
+ ind = this->m_Image->ComputeIndex( static_cast<IndexValueType>(this->m_Offset) );
+
+ const typename itk::ImageIterator<TImage>::IndexType&
+ startIndex = this->m_Region.GetIndex();
+
+ // Deccrement along a row, then wrap at the beginning of the region row.
+ bool done;
+ unsigned int dim;
+
+ // Check to see if we are past the first pixel in the region
+ // Note that --ind[0] moves to the previous pixel along the row.
+ ind[0] -= static_cast< IndexValueType >( m_SubsampleFactor );
+ done = (ind[0] <= startIndex[0] - 1);
+ for (unsigned int i=1; done && i < ImageIteratorDimension; i++)
+ {
+ done = (ind[i] <= startIndex[i]);
+ }
+
+ // if the iterator is outside the region (but not past region begin) then
+ // we need to wrap around the region
+ dim = 0;
+ if (!done)
+ {
+ while ( (dim < ImageIteratorDimension - 1)
+ && (ind[dim] < startIndex[dim]) )
+ {
+ ind[dim] = m_LastUsableIndex[dim];
+ ind[++dim] -= static_cast<IndexValueType>( m_SubsampleFactor );
+ }
+ }
+ this->m_Offset = this->m_Image->ComputeOffset( ind );
+ this->m_SpanEndOffset = this->m_Offset + 1;
+ this->m_SpanBeginOffset = this->m_Offset - m_LastUsableIndex[0];
+}
+
+} // end of namespace otb
+
+
+#endif
+
diff --git a/Code/MultiScale/otbHaarOperator.h b/Code/MultiScale/otbHaarOperator.h
index cd96cf0c4e3f5af03b1b3d52c3b6c44946b7f302..31377b73fcfdffebb8e4694ba2b97f41b185aa5b 100644
--- a/Code/MultiScale/otbHaarOperator.h
+++ b/Code/MultiScale/otbHaarOperator.h
@@ -43,7 +43,8 @@ namespace otb {
*
* \ingroup Operators
*/
-template<class TPixel, unsigned int VDimension,
+template < InverseOrForwardTransformationEnum TDirectionOfTransformation,
+ class TPixel, unsigned int VDimension,
class TAllocator = itk::NeighborhoodAllocator< TPixel > >
class ITK_EXPORT LowPassHaarOperator
: public WaveletOperator<TPixel, VDimension, TAllocator>
@@ -55,16 +56,20 @@ public:
itkTypeMacro(LowPassHaarOperator, WaveletOperator);
+ typedef InverseOrForwardTransformationEnum DirectionOfTransformationEnumType;
+ itkStaticConstMacro(DirectionOfTransformation,DirectionOfTransformationEnumType,TDirectionOfTransformation);
+
LowPassHaarOperator()
- {
+ {
this->SetRadius(1);
- this->CreateToRadius(1);
+ this->CreateToRadius(1);
+ this->SetWavelet( "Haar" );
}
LowPassHaarOperator(const Self& other)
: WaveletOperator<TPixel, VDimension, TAllocator>(other)
{
- /* ... */
+ this->SetWavelet( "Haar" );
}
/**
@@ -105,6 +110,17 @@ protected:
// stands for z^1
coeff.push_back(0.0);
+#if 0
+ std::cerr << "Coeff H(" << this->GetWavelet();
+ if ( (int) DirectionOfTransformation == (int) FORWARD )
+ std::cerr << " Forward ) = ";
+ else
+ std::cerr << " Inverse ) = ";
+ for ( typename CoefficientVector::const_iterator iter = coeff.begin(); iter < coeff.end(); ++iter )
+ std::cerr << *iter << " ";
+ std::cerr << "\n";
+#endif
+
return Superclass::UpSamplingCoefficients( coeff );
}
@@ -134,7 +150,8 @@ protected:
*
* \ingroup Operators
*/
-template<class TPixel, unsigned int VDimension,
+template< InverseOrForwardTransformationEnum TDirectionOfTransformation,
+ class TPixel, unsigned int VDimension,
class TAllocator = itk::NeighborhoodAllocator< TPixel > >
class ITK_EXPORT HighPassHaarOperator
: public WaveletOperator<TPixel, VDimension, TAllocator>
@@ -146,16 +163,20 @@ public:
itkTypeMacro(HighPassHaarOperator, WaveletOperator);
+ typedef InverseOrForwardTransformationEnum DirectionOfTransformationEnumType;
+ itkStaticConstMacro(DirectionOfTransformation,DirectionOfTransformationEnumType,TDirectionOfTransformation);
+
HighPassHaarOperator()
- {
+ {
this->SetRadius(1);
- this->CreateToRadius(1);
+ this->CreateToRadius(1);
+ this->SetWavelet( "Haar" );
}
HighPassHaarOperator(const Self& other)
: WaveletOperator<TPixel, VDimension, TAllocator>(other)
{
- /* ... */
+ this->SetWavelet( "Haar" );
}
/**
@@ -185,13 +206,49 @@ protected:
/**
* Set operator coefficients.
+ * Due to the number of template parameter, we cannot specialized it with regard
+ * to the TDirectionOfTransformation value...
*/
CoefficientVector GenerateCoefficients()
{
CoefficientVector coeff;
- coeff.push_back(-0.5);
- coeff.push_back(0.5);
- coeff.push_back(0.0);
+
+ typedef LowPassHaarOperator< FORWARD, TPixel, VDimension, TAllocator >
+ LowPassOperatorType;
+ LowPassOperatorType lowPassOperator;
+ lowPassOperator.SetDirection(0);
+ lowPassOperator.CreateDirectional();
+
+ CoefficientVector lowPassCoeff;
+ lowPassCoeff.resize( lowPassOperator.GetSize()[0] );
+
+ for ( typename LowPassOperatorType::ConstIterator iter = lowPassOperator.Begin();
+ iter != lowPassOperator.End(); ++iter )
+ lowPassCoeff.push_back( *iter );
+
+ switch ( DirectionOfTransformation )
+ {
+ case FORWARD:
+ coeff = this->GetHighPassFilterFromLowPassFilter( lowPassCoeff );
+ break;
+ case INVERSE:
+ coeff = this->GetInverseHighPassFilterFromForwardLowPassFilter( lowPassCoeff );
+ break;
+ default:
+ itkExceptionMacro(<<"Wavelet operator has to be INVERSE or FORWARD only!!");
+ break;
+ }
+
+#if 1
+ std::cerr << "Coeff G(" << this->GetWavelet();
+ if ( (int) DirectionOfTransformation == (int) FORWARD )
+ std::cerr << " Forward ) = ";
+ else
+ std::cerr << " Inverse ) = ";
+ for ( typename CoefficientVector::const_iterator iter = coeff.begin(); iter < coeff.end(); ++iter )
+ std::cerr << *iter << " ";
+ std::cerr << "\n";
+#endif
return Superclass::UpSamplingCoefficients( coeff );
}
@@ -203,7 +260,6 @@ protected:
}
};
-
} // end of namespace otb
#endif
diff --git a/Code/MultiScale/otbSplineBiOrthogonalOperator.h b/Code/MultiScale/otbSplineBiOrthogonalOperator.h
new file mode 100644
index 0000000000000000000000000000000000000000..99dc9d354cb0202c64af48c459a0561c0dac1b32
--- /dev/null
+++ b/Code/MultiScale/otbSplineBiOrthogonalOperator.h
@@ -0,0 +1,178 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+ Copyright (c) Institut Telecom / Telecom Bretagne. All rights reserved.
+ See ITCopyright.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.
+
+=========================================================================*/
+#ifndef __otbSplieBiOrthogonalOperator_h
+#define __otbSplieBiOrthogonalOperator_h
+
+#include "itkMacro.h"
+#include "itkExceptionObject.h"
+
+#include "otbWaveletOperator.h"
+
+namespace otb {
+
+/**
+ * \class LowPassSplineBiOrthogonalOperator
+ *
+ * \brief A NeighborhoodOperator for performing a Spline Bi-Orthogonal-based filtering
+ * at a pixel location.
+ *
+ * It may be configure with many filter orders. AT present time, so-called "9/7" filter is
+ * implemented. It may be fixed with SetWavelet method
+ *
+ * The LowPassSplineBiOrthogonalOperator is a NeighborhoodOperator that should be applied a
+ * NeighborhoodIterator using the NeighborhoodInnerProduct method.
+ * The 9/7 Operator is defiend in 1D as
+ * \f$ H(z) = 0.026748757411 z^{-4} -0.016864118443 z^{-3} -0.078223266529 z^{-2}
+ * + 0.266864118443 z^{-1} + 0.602949018236 + 0.266864118443 z
+ * -0.078223266529 z^2 -0.016864118443 z^3 + 0.026748757411 z^4 \f$.
+ *
+ * In N dimensions, the operator is directional
+ *
+ * \sa WaveletOperator
+ * \sa NeighborhoodOperator
+ * \sa Neighborhood
+ * \sa ForwardDifferenceOperator
+ * \sa BackwardDifferenceOperator
+ *
+ * \ingroup Operators
+ */
+template < InverseOrForwardTransformationEnum TDirectionOfTransformation,
+ class TPixel, unsigned int VDimension,
+ class TAllocator = itk::NeighborhoodAllocator< TPixel > >
+class ITK_EXPORT LowPassSplineBiOrthogonalOperator
+ : public WaveletOperator<TPixel, VDimension, TAllocator>
+{
+public:
+ /** Standard typedefs */
+ typedef LowPassSplineBiOrthogonalOperator Self;
+ typedef WaveletOperator<TPixel, VDimension, TAllocator> Superclass;
+
+ itkTypeMacro(LowPassSplineBiOrthogonalOperator, WaveletOperator);
+
+ typedef InverseOrForwardTransformationEnum DirectionOfTransformationEnumType;
+ itkStaticConstMacro(DirectionOfTransformation,DirectionOfTransformationEnumType,TDirectionOfTransformation);
+
+ LowPassSplineBiOrthogonalOperator() ;
+ LowPassSplineBiOrthogonalOperator(const Self& other);
+
+ /**
+ * Assignment operator
+ */
+ Self &operator=(const Self& other);
+
+ /**
+ * Prints some debugging information
+ */
+ virtual void PrintSelf(std::ostream &os, itk::Indent i) const ;
+
+protected:
+ /**
+ * Typedef support for coefficient vector type. Necessary to
+ * work around compiler bug on VC++.
+ */
+ typedef typename Superclass::CoefficientVector CoefficientVector;
+ typedef typename Superclass::PixelType PixelType;
+
+ /**
+ * Set operator coefficients.
+ */
+ CoefficientVector GenerateCoefficients();
+
+ /** Arranges coefficients spatially in the memory buffer. */
+ void Fill(const CoefficientVector& coeff);
+};
+
+/**
+ * \class HighPassSplineBiOrthogonalOperator
+ *
+ * \brief A NeighborhoodOperator for performing a Spline Bi-Orthogonal-based filtering
+ * at a pixel location.
+ *
+ * HighPassSplineBiOrthogonalOperator is a NeighborhoodOperator that should be applied a
+ * NeighborhoodIterator using the NeighborhoodInnerProduct method.
+ * The 9/7 Operator is defiend in 1D as
+ * \f$ H(z) = 0.045635881557 z^{-3} -0.028771763114 z^{-2} -0.295635881557 z^{-1}
+ * + 0.557543526229 - 0.295635881557 z -0.028771763114 z^2 + 0.045635881557 z^3 \f$.
+ *
+ * In N dimensions, the operator is directional
+ *
+ * \sa WaveletOperator
+ * \sa NeighborhoodOperator
+ * \sa Neighborhood
+ * \sa ForwardDifferenceOperator
+ * \sa BackwardDifferenceOperator
+ *
+ * \ingroup Operators
+ */
+template < InverseOrForwardTransformationEnum TDirectionOfTransformation,
+class TPixel, unsigned int VDimension,
+ class TAllocator = itk::NeighborhoodAllocator< TPixel > >
+class ITK_EXPORT HighPassSplineBiOrthogonalOperator
+ : public WaveletOperator<TPixel, VDimension, TAllocator>
+{
+public:
+ /** Standard typedefs */
+ typedef HighPassSplineBiOrthogonalOperator Self;
+ typedef WaveletOperator<TPixel, VDimension, TAllocator> Superclass;
+
+ itkTypeMacro(HighPassSplineBiOrthogonalOperator, WaveletOperator);
+
+ typedef InverseOrForwardTransformationEnum DirectionOfTransformationEnumType;
+ itkStaticConstMacro(DirectionOfTransformation,DirectionOfTransformationEnumType,TDirectionOfTransformation);
+
+ HighPassSplineBiOrthogonalOperator () ;
+ HighPassSplineBiOrthogonalOperator(const Self& other);
+
+ /**
+ * Assignment operator
+ */
+ Self &operator=(const Self& other);
+
+ /**
+ * Prints some debugging information
+ */
+ virtual void PrintSelf(std::ostream &os, itk::Indent i) const ;
+
+protected:
+ /**
+ * Typedef support for coefficient vector type. Necessary to
+ * work around compiler bug on VC++.
+ */
+ typedef typename Superclass::CoefficientVector CoefficientVector;
+ typedef typename Superclass::PixelType PixelType;
+
+ /**
+ * Set operator coefficients.
+ */
+ CoefficientVector GenerateCoefficients();
+
+ /** Arranges coefficients spatially in the memory buffer. */
+ void Fill(const CoefficientVector& coeff);
+};
+
+} // end of namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSplineBiOrthogonalOperator.txx"
+#endif
+
+
+#endif
+
diff --git a/Code/MultiScale/otbStationaryFilterBank.h b/Code/MultiScale/otbStationaryFilterBank.h
index a45ed5b9b775d4ee4e0c31abfcc5c2e7b560f1bd..ac7141b235e213f2adf2165ae39bea200614c59d 100644
--- a/Code/MultiScale/otbStationaryFilterBank.h
+++ b/Code/MultiScale/otbStationaryFilterBank.h
@@ -38,7 +38,8 @@ namespace otb {
* (ie. convolution-like operation).
*
* the inner operator are supposed to be defined through 1D filters. Then, the
- * transformation yields \f$ 2^Dim \f$ output images.
+ * forward transformation yields \f$ 2^{\test{Dim}} \f$ output images, while the inverse
+ * transformation requires \f$ 2^{\text{Dim}} \f$ input image for one output.
*
* In case of 1D, GetOutput(0) -> LowPass
*
@@ -63,7 +64,12 @@ namespace otb {
*
* -> x_(n-1) x_(n-2) x_(n-3) x_1 x_0
*
+ * And conversely in the inverse transformation.
*
+ * TODO: At present version, there is not consideration on meta data information that can be transmited
+ * from the input(s) to the output(s)...
+ *
+ * TODO: In the case of multiresolution decomposition, take care of the size of outputRegionForThread...
*
* \sa LowPassHaarOperator
* \sa HighPassHaarOperator
@@ -72,7 +78,9 @@ namespace otb {
*
* \ingroup Streamed
*/
-template < class TInputImage, class TOutputImage, class TLowPassOperator, class THighPassOperator >
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
class ITK_EXPORT StationaryFilterBank
: public itk::ImageToImageFilter< TInputImage, TOutputImage >
{
@@ -81,68 +89,113 @@ public:
typedef StationaryFilterBank Self;
typedef itk::ImageToImageFilter< TInputImage, TOutputImage > Superclass;
typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
- /** Type macro */
- itkNewMacro(Self);
+ /** Type macro */
+ itkNewMacro(Self);
- /** Creation through object factory macro */
- itkTypeMacro(StationaryFilterBank,ImageToImageFilter);
+ /** Creation through object factory macro */
+ itkTypeMacro(StationaryFilterBank,ImageToImageFilter);
/** Template parameters typedefs */
typedef TInputImage InputImageType;
- typedef typename InputImageType::Pointer InputImagePointerType;
- typedef typename InputImageType::RegionType InputImageRegionType;
- typedef typename InputImageType::SizeType InputSizeType;
- typedef typename InputImageType::IndexType InputIndexType;
- typedef typename InputImageType::PixelType InputPixelType;
+ typedef typename InputImageType::Pointer InputImagePointerType;
+ typedef typename InputImageType::RegionType InputImageRegionType;
+ typedef typename InputImageType::SizeType InputSizeType;
+ typedef typename InputImageType::IndexType InputIndexType;
+ typedef typename InputImageType::PixelType InputPixelType;
typedef TOutputImage OutputImageType;
- typedef typename OutputImageType::Pointer OutputImagePointerType;
+ typedef typename OutputImageType::Pointer OutputImagePointerType;
typedef typename OutputImageType::RegionType OutputImageRegionType;
typedef typename OutputImageType::SizeType OutputSizeType;
typedef typename OutputImageType::IndexType OutputIndexType;
- typedef typename OutputImageType::PixelType OutputPixelType;
+ typedef typename OutputImageType::PixelType OutputPixelType;
typedef TLowPassOperator LowPassOperatorType;
typedef THighPassOperator HighPassOperatorType;
+ typedef InverseOrForwardTransformationEnum DirectionOfTransformationEnumType;
+ itkStaticConstMacro(DirectionOfTransformation,DirectionOfTransformationEnumType,TDirectionOfTransformation);
+
/** Inner product iterators */
- typedef itk::ConstNeighborhoodIterator< OutputImageType > NeighborhoodIteratorType;
- typedef itk::NeighborhoodInnerProduct< OutputImageType > InnerProductType;
- typedef itk::ImageRegionIterator< OutputImageType > IteratorType;
+ typedef itk::ConstNeighborhoodIterator< InputImageType > NeighborhoodIteratorType;
+ typedef itk::NeighborhoodInnerProduct< InputImageType > InnerProductType;
+ typedef itk::ImageRegionIterator< InputImageType > IteratorType;
typedef typename itk::NeighborhoodAlgorithm
- ::ImageBoundaryFacesCalculator< OutputImageType > FaceCalculatorType;
+ ::ImageBoundaryFacesCalculator< InputImageType > FaceCalculatorType;
typedef typename FaceCalculatorType::FaceListType FaceListType;
typedef typename FaceListType::iterator FaceListIterator;
/** Dimension */
- itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
- itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension);
+ itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
+ itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension);
/**
- * Set/Get the level of up sampling of the filter used in the A-trou algorithm
+ * Set/Get the level of up sampling of the filter used in the A-trou algorithm.
*/
- itkGetMacro(UpSampleFactor,unsigned int);
- itkSetMacro(UpSampleFactor,unsigned int);
+ itkGetMacro(UpSampleFilterFactor,unsigned int);
+ itkSetMacro(UpSampleFilterFactor,unsigned int);
+
+ /**
+ * Set/Get the level of down sampling of the image used in forward algorithm.
+ * (or upsampling in the inverse case)
+ *
+ * In this implementation, we are dealing with M-band decomposition then m_SubSampleImageFactor
+ * is most likely to be 1 or 2... but in any case integer and not real...
+ */
+ itkGetMacro(SubSampleImageFactor,unsigned int);
+ itkSetMacro(SubSampleImageFactor,unsigned int);
+
protected:
StationaryFilterBank();
virtual ~StationaryFilterBank() { }
+ /** GenerateOutputInformation
+ * Set the size of the output image depending on the decimation factor
+ * Copy informations from the input image if existing.
+ **/
+ virtual void GenerateOutputInformation();
+
+ /** BeforeThreadedGenerateData
+ * If SubSampleImageFactor neq 1, it is necessary to up sample input images in the INVERSE mode
+ */
+ virtual void BeforeThreadedGenerateData ();
+
+ /** AfterThreadedGenerateData
+ * If SubSampleImageFactor neq 1, it is necessary to down sample output images in the FORWARD mode
+ */
+ virtual void AfterThreadedGenerateData ();
+
/** Generate data redefinition */
virtual void ThreadedGenerateData ( const OutputImageRegionType& outputRegionForThread, int threadId );
- /** Specific case applied on the input image */
- virtual void ThreadedGenerateData ( itk::ProgressReporter & reporter,
- const OutputImageRegionType& outputRegionForThread, int threadId );
+ /** Specification-like in the case where DirectionOfTransformation stands for FORWARD or INVERSE */
+ virtual void ThreadedForwardGenerateData ( const OutputImageRegionType& outputRegionForThread, int threadId );
+ virtual void ThreadedInverseGenerateData ( const OutputImageRegionType& outputRegionForThread, int threadId );
- /** Iterative call to the filter bank at each dimension */
- virtual void ThreadedGenerateData ( unsigned int idx, unsigned int direction,
- itk::ProgressReporter & reporter,
- const OutputImageRegionType& outputRegionForThread, int threadId );
+ /** Specific case applied on the input image in forward transformation */
+ virtual void ThreadedForwardGenerateData ( itk::ProgressReporter & reporter,
+ const OutputImageRegionType& outputRegionForThread, int threadId );
+
+ /** Iterative call to the forward filter bank at each dimension */
+ virtual void ThreadedForwardGenerateData ( unsigned int idx, unsigned int direction,
+ itk::ProgressReporter & reporter,
+ const OutputImageRegionType& outputRegionForThread, int threadId );
+
+ /** Iterative call to the inverse filter bank at each dimension */
+ virtual void ThreadedInverseGenerateData ( unsigned int idx, unsigned int direction,
+ InputImagePointerType & outputImage,
+ itk::ProgressReporter & reporter,
+ const InputImageRegionType& inputRegionForThread, int threadId );
+
+ /** Specific case applied on the last reconstruction in inverse transformation */
+ virtual void ThreadedInverseGenerateData ( itk::ProgressReporter & reporter,
+ const OutputImageRegionType& outputRegionForThread, int threadId );
- unsigned int m_UpSampleFactor;
+ unsigned int m_UpSampleFilterFactor;
+ unsigned int m_SubSampleImageFactor;
private:
StationaryFilterBank( const Self & );
diff --git a/Code/MultiScale/otbStationaryFilterBank.txx b/Code/MultiScale/otbStationaryFilterBank.txx
index f1c80aa58632b29bb0c7d56ea0af8255e930377e..8f1a2800145bfca0a3ba3511683d73712bce2f17 100644
--- a/Code/MultiScale/otbStationaryFilterBank.txx
+++ b/Code/MultiScale/otbStationaryFilterBank.txx
@@ -23,44 +23,183 @@
#include "otbStationaryFilterBank.h"
+#include "otbMacro.h"
+
#include "itkNeighborhoodAlgorithm.h"
#include "itkZeroFluxNeumannBoundaryCondition.h"
#include "itkProgressReporter.h"
namespace otb {
-template < class TInputImage, class TOutputImage, class TLowPassOperator, class THighPassOperator >
-StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOperator >
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
::StationaryFilterBank ()
{
- this->SetNumberOfRequiredInputs(1);
+ switch ( DirectionOfTransformation )
+ {
+ case FORWARD:
+ {
+ this->SetNumberOfRequiredInputs(1);
+
+ unsigned int numOfOutputs = 1<<InputImageType::ImageDimension;
+
+ this->SetNumberOfOutputs( numOfOutputs );
+ for ( unsigned i = 1; i < numOfOutputs; i++ )
+ {
+ this->SetNthOutput(i,OutputImageType::New());
+ }
+
+ break;
+ }
+ case INVERSE:
+ {
+ this->SetNumberOfRequiredOutputs(1);
+ this->SetNumberOfInputs( 1 << InputImageType::ImageDimension );
+ break;
+ }
+ default:
+ {
+ itkExceptionMacro(<<"FilterBank transformation may be FORWARD or INVERSE only!!");
+ break;
+ }
+ }
+ m_UpSampleFilterFactor = 0;
+ m_SubSampleImageFactor = 1;
+}
- unsigned int numOfOutput = 1<<InputImageType::ImageDimension;
- this->SetNumberOfOutputs( numOfOutput );
- for ( unsigned i = 1; i < numOfOutput; i++ )
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
+void
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
+::GenerateOutputInformation()
+{
+ Superclass::GenerateOutputInformation();
+
+ switch ( DirectionOfTransformation )
+ {
+ case FORWARD:
+ {
+ /**
+ * If necessary, we do not say to the filter, at this
+ * step, that we will down sample the output images...
+ */
+ break;
+ }
+ case INVERSE:
+ {
+ otbGenericMsgDebugMacro( << " down sampling output regions by a factor of " << GetSubSampleImageFactor() );
+
+ this->GetOutput()->CopyInformation( this->GetInput() );
+ this->GetOutput()->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
+
+ if ( GetSubSampleImageFactor() > 1 )
+ {
+ for ( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ this->GetOutput()->GetRegion().GetIndex()[i] /= GetSubSampleImageFactor();
+ this->GetOutput()->GetRegion().GetSize()[i] /= GetSubSampleImageFactor();
+ }
+ }
+ break;
+ }
+ default:
+ {
+ itkExceptionMacro(<<"FilterBank transformation may be FORWARD or INVERSE only!!");
+ break;
+ }
+ }
+}
+
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
+void
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
+::BeforeThreadedGenerateData ()
+{
+ if ( static_cast<int>( DirectionOfTransformation ) == INVERSE )
{
- this->SetNthOutput(i,OutputImageType::New());
+ otbGenericMsgDebugMacro(<<"Subsample inputs by a factor of " << GetSubSampleImageFactor() );
+
+ // Sub sample inputs
}
+}
- m_UpSampleFactor = 0;
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
+void
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
+::AfterThreadedGenerateData ()
+{
+ if ( static_cast<int>( DirectionOfTransformation ) == FORWARD )
+ {
+ otbGenericMsgDebugMacro(<<"Downsample outputs by a factor of " << GetSubSampleImageFactor() );
+
+ // Down sample inputs
+ }
}
-template < class TInputImage, class TOutputImage, class TLowPassOperator, class THighPassOperator >
+
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
void
-StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOperator >
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
::ThreadedGenerateData ( const OutputImageRegionType& outputRegionForThread, int threadId )
{
+ switch ( DirectionOfTransformation )
+ {
+ case FORWARD:
+ return ThreadedForwardGenerateData( outputRegionForThread, threadId );
+ break;
+ case INVERSE:
+ return ThreadedInverseGenerateData( outputRegionForThread, threadId );
+ break;
+ default:
+ itkExceptionMacro(<<"FilterBank transformation may be FORWARD or INVERSE only!!");
+ break;
+ }
+}
+
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
+void
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
+::ThreadedForwardGenerateData ( const OutputImageRegionType& outputRegionForThread, int threadId )
+{
+
itk::ProgressReporter reporter ( this, threadId,
outputRegionForThread.GetNumberOfPixels() * this->GetNumberOfOutputs() );
- ThreadedGenerateData( reporter, outputRegionForThread, threadId );
+ ThreadedForwardGenerateData( reporter, outputRegionForThread, threadId );
}
-template < class TInputImage, class TOutputImage, class TLowPassOperator, class THighPassOperator >
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
void
-StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOperator >
-::ThreadedGenerateData
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
+::ThreadedForwardGenerateData
( itk::ProgressReporter & reporter,
const OutputImageRegionType& outputRegionForThread, int threadId )
{
@@ -79,8 +218,8 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
if ( outputHighPass == 0 )
{
itk::OStringStream msg;
- msg << "Sortie 1<<" << idx << " = " << (1<<idx) << " nulle\n";
- msg << "Nombre de sortie attendue " << this->GetNumberOfOutputs() << "\n";
+ msg << "Output number 1<<" << idx << " = " << (1<<idx) << " null\n";
+ msg << "Number of excpected outputs " << this->GetNumberOfOutputs() << "\n";
throw itk::ExceptionObject( __FILE__, __LINE__, msg.str().c_str(), ITK_LOCATION );
}
@@ -97,7 +236,7 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
// High pass part calculation
HighPassOperatorType highPassOperator;
highPassOperator.SetDirection(idx);
- highPassOperator.SetUpSampleFactor( this->GetUpSampleFactor() );
+ highPassOperator.SetUpSampleFactor( this->GetUpSampleFilterFactor() );
highPassOperator.CreateDirectional();
faceList = faceCalculator( input, outputRegionForThread, highPassOperator.GetRadius() );
@@ -117,7 +256,7 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
// Low pass part calculation
LowPassOperatorType lowPassOperator;
lowPassOperator.SetDirection(idx);
- lowPassOperator.SetUpSampleFactor( this->GetUpSampleFactor() );
+ lowPassOperator.SetUpSampleFactor( this->GetUpSampleFilterFactor() );
lowPassOperator.CreateDirectional();
faceList = faceCalculator( input, outputRegionForThread, lowPassOperator.GetRadius() );
@@ -136,15 +275,19 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
if ( idx > 0 )
{
- ThreadedGenerateData( 0, idx-1, reporter, outputRegionForThread, threadId );
- ThreadedGenerateData( 1<<idx, idx-1, reporter, outputRegionForThread, threadId );
+ ThreadedForwardGenerateData( 0, idx-1, reporter, outputRegionForThread, threadId );
+ ThreadedForwardGenerateData( 1<<idx, idx-1, reporter, outputRegionForThread, threadId );
}
}
-template < class TInputImage, class TOutputImage, class TLowPassOperator, class THighPassOperator >
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
void
-StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOperator >
-::ThreadedGenerateData
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
+::ThreadedForwardGenerateData
( unsigned int idx, unsigned int direction, itk::ProgressReporter & reporter,
const OutputImageRegionType& outputRegionForThread, int threadId )
{
@@ -162,7 +305,7 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
// High pass part calculation
HighPassOperatorType highPassOperator;
highPassOperator.SetDirection(direction);
- highPassOperator.SetUpSampleFactor( this->GetUpSampleFactor() );
+ highPassOperator.SetUpSampleFactor( this->GetUpSampleFilterFactor() );
highPassOperator.CreateDirectional();
faceList = faceCalculator( input, outputRegionForThread, highPassOperator.GetRadius() );
@@ -182,7 +325,7 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
// Low pass part calculation
LowPassOperatorType lowPassOperator;
lowPassOperator.SetDirection(direction);
- lowPassOperator.SetUpSampleFactor( this->GetUpSampleFactor() );
+ lowPassOperator.SetUpSampleFactor( this->GetUpSampleFilterFactor() );
lowPassOperator.CreateDirectional();
faceList = faceCalculator( input, outputRegionForThread, lowPassOperator.GetRadius() );
@@ -205,7 +348,7 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
for ( inIt.GoToBegin(), outIt.GoToBegin(); !inIt.IsAtEnd(); ++inIt, ++outIt )
{
- outIt.Set( inIt.Get() );
+ outIt.Set( static_cast<InputPixelType>( inIt.Get() ) );
}
if ( direction != 0 )
@@ -215,6 +358,206 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
}
}
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
+void
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
+::ThreadedInverseGenerateData ( const OutputImageRegionType& outputRegionForThread, int threadId )
+{
+ itk::ProgressReporter reporter ( this, threadId,
+ outputRegionForThread.GetNumberOfPixels() * this->GetNumberOfInputs() );
+
+ ThreadedInverseGenerateData( reporter, outputRegionForThread, threadId );
+}
+
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
+void
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
+::ThreadedInverseGenerateData
+( itk::ProgressReporter & reporter,
+ const OutputImageRegionType& outputRegionForThread, int threadId )
+{
+ InputImageRegionType inputRegionForThread;
+ this->CallCopyOutputRegionToInputRegion( inputRegionForThread, outputRegionForThread );
+
+ unsigned int idx = InputImageType::ImageDimension-1;
+
+ InputImagePointerType imgLowPassPointer = InputImageType::New();
+ InputImageType * imgLowPass = this->GetInput(0);
+ if ( idx != 0 )
+ {
+ imgLowPassPointer->SetRegions( inputRegionForThread );
+ imgLowPassPointer->Allocate();
+
+ ThreadedInverseGenerateData( 0, idx-1, imgLowPassPointer, reporter, inputRegionForThread, threadId );
+
+ imgLowPass = imgLowPassPointer.GetPointer();
+ }
+
+ InputImagePointerType imgHighPassPointer = InputImageType::New();
+ InputImageType * imgHighPass = this->GetInput(1);
+ if ( idx != 0 )
+ {
+ imgHighPassPointer->SetRegions( inputRegionForThread );
+ imgHighPassPointer->Allocate();
+
+ ThreadedInverseGenerateData( 1<<idx, idx-1, imgHighPassPointer, reporter, inputRegionForThread, threadId );
+
+ imgHighPass = imgHighPassPointer.GetPointer();
+ }
+
+ // Low pass part calculation
+ LowPassOperatorType lowPassOperator;
+ lowPassOperator.SetDirection(idx);
+ lowPassOperator.SetUpSampleFactor( this->GetUpSampleFilterFactor() );
+ lowPassOperator.CreateDirectional();
+
+ // High pass part calculation
+ HighPassOperatorType highPassOperator;
+ highPassOperator.SetDirection(idx);
+ highPassOperator.SetUpSampleFactor( this->GetUpSampleFilterFactor() );
+ highPassOperator.CreateDirectional();
+
+ // Faces iterations
+ typedef NeighborhoodIteratorType::RadiusType radiusMax;
+ for ( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ radiusMax[i] = lowPassOperator.GetRadius()[i];
+ if ( radius[i] < highPassOperator.GetRadius()[i] )
+ radius[i] = highPassOperator.GetRadius()[i];
+ }
+
+ FaceCalculatorType faceCalculator;
+ FaceListType faceList;
+ faceList = faceCalculator( imgHighPass, inputRegionForThread, radiusMax );
+
+ InnerProductType innerProduct;
+
+ for ( FaceListIterator faceIt = faceList.begin(); faceIt != faceList.end(); ++faceIt )
+ {
+ itk::ImageRegionIterator< OutputImageType > out ( this->GetOutput(), *faceIt );
+
+ NeighborhoodIteratorType lowIter ( lowPassOperator.GetRadius(), imgLowPass, *faceIt );
+ NeighborhoodIteratorType highIter ( highPassOperator.GetRadius(), imgHighPass, *faceIt );
+
+ lowIter.GoToBegin();
+ highIter.GoToBegin();
+ out.GoToBegin();
+
+ while ( !out.IsAtEnd() )
+ {
+ out.Set( static_cast< OutputPixelType >(
+ innerProduct( lowIter, lowPassOperator )
+ + innerProduct( highIter, highPassOperator ) ) );
+
+ ++lowIter;
+ ++highIter;
+ ++out;
+
+ reporter.CompletedPixel();
+ }
+ }
+}
+
+template < class TInputImage, class TOutputImage,
+ class TLowPassOperator, class THighPassOperator,
+ InverseOrForwardTransformationEnum TDirectionOfTransformation >
+void
+StationaryFilterBank< TInputImage, TOutputImage,
+ TLowPassOperator, THighPassOperator,
+ TDirectionOfTransformation >
+::ThreadedInverseGenerateData
+( unsigned int idx, unsigned int direction,
+ InputImagePointerType & outputImage,
+ itk::ProgressReporter & reporter,
+ const InputImageRegionType& inputRegionForThread, int threadId )
+{
+ InputImagePointerType imgLowPassPointer = InputImageType::New();
+ InputImageType * imgLowPass = this->GetInput(idx);
+ if ( direction != 0 )
+ {
+ imgLowPassPointer->SetRegions( inputRegionForThread );
+ imgLowPassPointer->Allocate();
+
+ ThreadedGenerateData( idx, direction-1, imgLowPassPointer, reporter, outputRegionForThread, threadId );
+
+ imgLowPass = imgLowPassPointer.GetPointer();
+ }
+
+ InputImagePointerType imgHighPassPointer = InputImageType::New();
+ InputImageType * imgHighPass = this->GetInput(idx + (1<<direction));
+ if ( direction != 0 )
+ {
+ imgHighPassPointer->SetReions( inputRegionForThread );
+ imgHighPassPointer->Allocate();
+
+ ThreadedInverseGenerateData( idx + (1<<direction), direction-1, imgHighPassPointer, reporter, inputRegionForThread, threadId );
+
+ imgHighPass = imgHighPassPointer.GetPointer();
+ }
+
+ // Low pass part calculation
+ LowPassOperatorType lowPassOperator;
+ lowPassOperator.SetDirection(idx);
+ lowPassOperator.SetUpSampleFactor( this->GetUpSampleFilterFactor() );
+ lowPassOperator.CreateDirectional();
+
+ // High pass part calculation
+ HighPassOperatorType highPassOperator;
+ highPassOperator.SetDirection(idx);
+ highPassOperator.SetUpSampleFactor( this->GetUpSampleFilterFactor() );
+ highPassOperator.CreateDirectional();
+
+ // Faces iterations
+ typedef NeighborhoodIteratorType::RadiusType radiusMax;
+ for ( unsigned int i = 0; i < ImageDimension; i++ )
+ {
+ radiusMax[i] = lowPassOperator.GetRadius()[i];
+ if ( radius[i] < highPassOperator.GetRadius()[i] )
+ radius[i] = highPassOperator.GetRadius()[i];
+ }
+
+ FaceCalculatorType faceCalculator;
+ FaceListType faceList;
+ faceList = faceCalculator( imgHighPass, inputRegionForThread, radiusMax );
+
+ InnerProductType innerProduct;
+
+ for ( FaceListIterator faceIt = faceList.begin(); faceIt != faceList.end(); ++faceIt )
+ {
+ IteratorType out ( this->GetOutput(), *faceIt );
+
+ NeighborhoodIteratorType lowIter ( lowPassOperator.GetRadius(), imgLowPass, *faceIt );
+ NeighborhoodIteratorType highIter ( highPassOperator.GetRadius(), imgHighPass, *faceIt );
+
+ lowIter.GoToBegin();
+ highIter.GoToBegin();
+ out.GoToBegin();
+
+ while ( !out.IsAtEnd() )
+ {
+ out.Set( innerProduct( lowIter, lowPassOperator )
+ + innerProduct( highIter, highPassOperator ) );
+
+ ++lowIter;
+ ++highIter;
+ ++out;
+
+ reporter.CompletedPixel();
+ }
+ }
+}
+
+
+
+
} // end of namespace
#endif
diff --git a/Code/MultiScale/otbWaveletForwardTransform.h b/Code/MultiScale/otbWaveletForwardTransform.h
index a94474b841ef417babbb8453dbe497d59ed37210..18554318e892de5bd3a59b4f2eadcc9fa097f9a9 100644
--- a/Code/MultiScale/otbWaveletForwardTransform.h
+++ b/Code/MultiScale/otbWaveletForwardTransform.h
@@ -54,13 +54,13 @@ public:
typedef WaveletForwardTransform Self;
typedef itk::ImageToImageFilter< TInputImage, TOutputImage > Superclass;
typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
- /** Type macro */
- itkNewMacro(Self);
+ /** Type macro */
+ itkNewMacro(Self);
- /** Creation through object factory macro */
- itkTypeMacro(WaveletForwardTransform,ImageToImageFilter);
+ /** Creation through object factory macro */
+ itkTypeMacro(WaveletForwardTransform,ImageToImageFilter);
typedef TInputImage InputImageType;
typedef typename InputImageType::Pointer InputImagePointer;
diff --git a/Code/MultiScale/otbWaveletOperator.h b/Code/MultiScale/otbWaveletOperator.h
index b5263022cb0a19463627421c73e90a0f7648e317..5fd869196756be64e28f22d6284c1a6e9cf54eea 100644
--- a/Code/MultiScale/otbWaveletOperator.h
+++ b/Code/MultiScale/otbWaveletOperator.h
@@ -23,6 +23,9 @@
#include "itkExceptionObject.h"
#include "itkNeighborhoodOperator.h"
+// This include is needed to define InverseOrForwardTransformationEnum only...
+#include "otbGenericMapProjection.h"
+
namespace otb {
/**
@@ -65,15 +68,16 @@ public:
itkTypeMacro(WaveletOperator,NeighborhoodOperator);
/** Construction */
- WaveletOperator()
- {
- m_UpSampleFactor = 0;
- }
+ WaveletOperator() :
+ m_UpSampleFactor( 0 ), m_Wavelet( "Unknown" )
+ { }
+
/** Construction by copy */
WaveletOperator( const Self & other )
: itk::NeighborhoodOperator<TPixel, VDimension, TAllocator> (other)
{
m_UpSampleFactor = other.GetUpSampleFactor();
+ m_Wavelet = other.m_Wavelet;
}
virtual ~WaveletOperator() {}
@@ -82,6 +86,7 @@ public:
{
Superclass::operator=(other);
m_UpSampleFactor = other.GetUpSampleFactor();
+ m_Wavelet = other.m_Wavelet;
return *this;
}
/**
@@ -91,6 +96,7 @@ public:
{
Superclass::PrintSelf(os, i.GetNextIndent());
os << i << "Up-Sampling factor " << this->m_UpSampleFactor << "\n";
+ os << i << "Wavelet kind : " << GetWavelet() << std::endl;
}
/**
@@ -105,6 +111,19 @@ public:
this->m_UpSampleFactor = upSampleFactor;
}
+ /**
+ * Set/Get the name of the wavelet when necessary
+ */
+ virtual void SetWavelet ( const std::string & str ) {
+ this->m_Wavelet = str;
+ }
+
+ virtual void SetWavelet ( const char * str ) {
+ this->m_Wavelet = str; }
+
+ virtual const char * GetWavelet () const {
+ return this->m_Wavelet.c_str(); }
+
protected:
/**
* Typedef support for coefficient vector type. Necessary to
@@ -141,7 +160,85 @@ protected:
return coeff;
}
+ /**
+ * Performs the definition of synthesis filter from analysis one.
+ * Input is the forward low pass filter coefficients.
+ * It performs \f$ {\tilde G}(z) = -H(-z) \f$.
+ */
+ CoefficientVector GetInverseHighPassFilterFromForwardLowPassFilter ( CoefficientVector & coeff )
+ {
+ unsigned long radius = static_cast<unsigned long>( coeff.size() );
+ this->SetRadius( radius );
+
+ unsigned long medianPosition = radius/2;
+
+ CoefficientVector highPassedCoeff;
+ highPassedCoeff = coeff;
+
+ for ( unsigned int i = 0; i < medianPosition; i+=2 )
+ {
+ highPassedCoeff[medianPosition + i] = - coeff[ medianPosition + i ];
+ highPassedCoeff[medianPosition - i] = - coeff[ medianPosition - i ];
+ }
+
+ coeff = highPassedCoeff;
+ return coeff;
+ }
+
+ /**
+ * Performs the definition of synthesis filter from analysis one.
+ * Input is the forward high pass filter coefficients.
+ * It performs \f$ {\tilde H}(z) = G(-z) \f$.
+ */
+ CoefficientVector GetInverseLowPassFilterFromForwardHighPassFilter ( CoefficientVector & coeff )
+ {
+ unsigned long radius = static_cast<unsigned long>( coeff.size() );
+ this->SetRadius( radius );
+
+ unsigned long medianPosition = radius/2;
+
+ CoefficientVector lowPassedCoeff;
+ lowPassedCoeff = coeff;
+
+ for ( unsigned int i = 1; i < medianPosition; i+=2 )
+ {
+ lowPassedCoeff[medianPosition + i] = - coeff[ medianPosition + i ];
+ lowPassedCoeff[medianPosition - i] = - coeff[ medianPosition - i ];
+ }
+
+ coeff = lowPassedCoeff;
+ return coeff;
+ }
+
+ /**
+ * Performs the definition of high pass filter from the low pass one
+ * in a Quadrature mirror filter bank framework. It is then valid to
+ * define high pass filters from orthogonal wavelets...
+ * Input is the forward low pass filter coefficients.
+ * It performs \f$ G(z) = H(-z) \f$.
+ */
+ CoefficientVector GetHighPassFilterFromLowPassFilter ( CoefficientVector & coeff )
+ {
+ unsigned long radius = static_cast<unsigned long>( coeff.size() );
+ this->SetRadius( radius );
+
+ unsigned long medianPosition = radius/2;
+
+ CoefficientVector highPassedCoeff;
+ highPassedCoeff = coeff;
+
+ for ( unsigned int i = 1; i < medianPosition; i+=2 )
+ {
+ highPassedCoeff[medianPosition + i] = - coeff[ medianPosition + i ];
+ highPassedCoeff[medianPosition - i] = - coeff[ medianPosition - i ];
+ }
+
+ coeff = highPassedCoeff;
+ return coeff;
+ }
+
unsigned int m_UpSampleFactor;
+ std::string m_Wavelet;
};
} // end of namespace otb