BUG: A trou algo for stationary wavelet paquets

Code/MultiScale/otbHaarOperator.h

@@ -55,7 +55,7 @@ public:
 
   itkTypeMacro(LowPassHaarOperator, WaveletOperator);
   
-  LowPassHaarOperator() : WaveletOperator ()
+  LowPassHaarOperator() 
 	{
     this->SetRadius(1);
 		this->CreateToRadius(1);

Code/MultiScale/otbStationaryFilterBank.h

@@ -120,6 +120,12 @@ public:
 	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
+   */
+  itkGetMacro(UpSampleFactor,unsigned int);
+  itkSetMacro(UpSampleFactor,unsigned int);
+
 protected:
   StationaryFilterBank();
   virtual ~StationaryFilterBank() { }
@@ -136,6 +142,7 @@ protected:
                                         itk::ProgressReporter & reporter,
                                         const OutputImageRegionType& outputRegionForThread, int threadId );
 
+  unsigned int m_UpSampleFactor;
 
 private:
   StationaryFilterBank( const Self & );

Code/MultiScale/otbStationaryFilterBank.txx

@@ -42,6 +42,8 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
   {
     this->SetNthOutput(i,OutputImageType::New());
   }
+
+  m_UpSampleFactor = 0;
 }
 
 template < class TInputImage, class TOutputImage, class TLowPassOperator, class THighPassOperator >
@@ -95,6 +97,7 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
   // High pass part calculation
   HighPassOperatorType highPassOperator;
   highPassOperator.SetDirection(idx);
+  highPassOperator.SetUpSampleFactor( this->GetUpSampleFactor() );
   highPassOperator.CreateDirectional();
 
   faceList = faceCalculator( input, outputRegionForThread, highPassOperator.GetRadius() );
@@ -114,6 +117,7 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
   // Low pass part calculation
   LowPassOperatorType lowPassOperator;
   lowPassOperator.SetDirection(idx);
+  lowPassOperator.SetUpSampleFactor( this->GetUpSampleFactor() );
   lowPassOperator.CreateDirectional();
 
   faceList = faceCalculator( input, outputRegionForThread, lowPassOperator.GetRadius() );
@@ -158,6 +162,7 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
   // High pass part calculation
   HighPassOperatorType highPassOperator;
   highPassOperator.SetDirection(direction);
+  highPassOperator.SetUpSampleFactor( this->GetUpSampleFactor() );
   highPassOperator.CreateDirectional();
 
   faceList = faceCalculator( input, outputRegionForThread, highPassOperator.GetRadius() );
@@ -177,6 +182,7 @@ StationaryFilterBank< TInputImage, TOutputImage, TLowPassOperator, THighPassOper
   // Low pass part calculation
   LowPassOperatorType lowPassOperator;
   lowPassOperator.SetDirection(direction);
+  lowPassOperator.SetUpSampleFactor( this->GetUpSampleFactor() );
   lowPassOperator.CreateDirectional();
 
   faceList = faceCalculator( input, outputRegionForThread, lowPassOperator.GetRadius() );

Code/MultiScale/otbWaveletForwardTransform.h

@@ -98,6 +98,12 @@ private:
   unsigned int m_NumberOfDecompositions;
   FilterListPointerType m_FilterList;
 
+  /*
+   * At this level, m_UseSubSampleImage is set to FALSE to perform
+   * multiscale analysis only
+   */
+  bool m_UseSubSampleImage;
+
 }; // end of class
 
 } // end of namespace

Code/MultiScale/otbWaveletForwardTransform.txx

@@ -34,6 +34,8 @@ WaveletForwardTransform< TInputImage, TOutputImage, TFilter >
 
   m_NumberOfDecompositions = 1;
   m_FilterList = FilterListType::New();
+
+  m_UseSubSampleImage = false;
 }
 
 template < class TInputImage, class TOutputImage, class TFilter >
@@ -86,6 +88,7 @@ WaveletForwardTransform< TInputImage, TOutputImage, TFilter >
   std::cerr << "outputIdxToGraft = " << outputIdxToGraft << "\n";
   */
 
+  unsigned int subSampleFactor = 0;
   while ( nbDecomp > 0 )
   { 
     /*
@@ -99,6 +102,9 @@ WaveletForwardTransform< TInputImage, TOutputImage, TFilter >
     while ( idx-- > 1 )
       GetFilter(nbDecomp)->GraftNthOutput( idx, this->GetOutput( outputIdxForGraft-- ) );
 
+    if ( m_UseSubSampleImage == false )
+      GetFilter(nbDecomp)->SetUpSampleFactor( ++subSampleFactor );
+
     GetFilter( nbDecomp )->Update();
 
     idx = this->GetFilter(nbDecomp)->GetNumberOfOutputs();
@@ -121,6 +127,9 @@ WaveletForwardTransform< TInputImage, TOutputImage, TFilter >
   while ( idx-- > 0 )
     GetFilter(nbDecomp)->GraftNthOutput( idx, this->GetOutput( outputIdxForGraft-- ) );
 
+  if ( m_UseSubSampleImage == false )
+      GetFilter(nbDecomp)->SetUpSampleFactor( ++subSampleFactor );
+
   GetFilter( nbDecomp )->Update();
 
   idx = this->GetFilter(nbDecomp)->GetNumberOfOutputs();
@@ -128,7 +137,6 @@ WaveletForwardTransform< TInputImage, TOutputImage, TFilter >
     this->GetOutput( outputIdxToGraft-- )->Graft( GetFilter(nbDecomp)->GetOutput( idx ) );
 }
 
-
 } // end of namespace otb
 
 #endif

Code/MultiScale/otbWaveletOperator.h

@@ -73,7 +73,7 @@ public:
   WaveletOperator( const Self & other ) 
     : itk::NeighborhoodOperator<TPixel, VDimension, TAllocator> (other)
   {
-    m_UpsampleFactor = other.GetUpSampleFactor();
+    m_UpSampleFactor = other.GetUpSampleFactor();
   }
   virtual ~WaveletOperator() {}
 
@@ -96,10 +96,16 @@ public:
   /**
    * Set/Get the level of up sampling of the filter
    */
-  itkGetMacro(UpSampleFactor,unsigned int);
-  itkSetMacro(UpSampleFactor,unsigned int);
+  unsigned int GetUpSampleFactor () const 
+  {
+    return this->m_UpSampleFactor;
+  }
+  void SetUpSampleFactor ( unsigned int upSampleFactor )
+  {
+    this->m_UpSampleFactor = upSampleFactor;
+  }
 
-    protected:
+protected:
   /**
    * Typedef support for coefficient vector type.  Necessary to
    * work around compiler bug on VC++.
@@ -113,13 +119,14 @@ public:
    */
   CoefficientVector UpSamplingCoefficients ( CoefficientVector & coeff ) 
   {
-    if ( m_UpsampleFactor == 0 )
+    if ( m_UpSampleFactor == 0 )
       return coeff;
 
-    for ( unsigned int up = 0; up < m_UpsampleFactor; ++up )
+    for ( unsigned int up = 0; up < m_UpSampleFactor; ++up )
     {
-      this->SetRadius( 2*this->GetRadius() );
-      this->CreateToRadius( this->GetRadius() );
+      unsigned long radius = static_cast<unsigned long>( coeff.size() ) >> 1;
+
+      this->SetRadius( 2*radius );
 
       CoefficientVector upSampledCoeff;
       upSampledCoeff.push_back( coeff[0] );
@@ -134,11 +141,11 @@ public:
     return coeff;
   }
 
-
-
   unsigned int m_UpSampleFactor;
 };
 
+} // end of namespace otb
+
 #endif
 
 

Code/MultiScale/otbWaveletPacketForwardTransform.txx

@@ -41,8 +41,6 @@ void
 WaveletPacketForwardTransform< TInputImage, TOutputImage, TFilter, TCost >
 ::GenerateData ()
 {
-  //this->AllocateOutputs();
-  
   /*
    * Start with a decomposition
    */
@@ -74,6 +72,7 @@ WaveletPacketForwardTransform< TInputImage, TOutputImage, TFilter, TCost >
     this->GetFilterList()->PushBack( FilterType::New() );
     FilterType * filter = this->GetFilterList()->GetNthElement( this->GetFilterList()->Size()-1 );
     filter->SetInput( outputIt.Get() );
+    filter->SetUpSampleFactor( depth );
     filter->Update();
 
     depth++;