Merge branch 'app_domain_transform' into develop

Modules/Applications/AppDomainTransform/CMakeLists.txt

+project(OTBAppDomainTransform)
+otb_module_impl()

Modules/Applications/AppDomainTransform/app/CMakeLists.txt

+otb_create_application(
+  NAME           DomainTransform
+  SOURCES        otbDomainTransform.cxx
+  LINK_LIBRARIES ${${otb-module}_LIBRARIES})

Modules/Applications/AppDomainTransform/app/otbDomainTransform.cxx

+/*=========================================================================
+
+  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.
+
+  =========================================================================*/
+
+#include "otbWrapperApplication.h"
+#include "otbWrapperApplicationFactory.h"
+#include "otbWaveletImageFilter.h"
+#include "otbWaveletInverseImageFilter.h"
+#include "otbWaveletGenerator.h"
+
+#include <itkConfigure.h>
+#include <itkForwardFFTImageFilter.h>
+#include <itkInverseFFTImageFilter.h>
+#include <itkUnaryFunctorImageFilter.h>
+#include <itkFFTShiftImageFilter.h>
+
+#include "otbComplexToVectorImageCastFilter.h"
+
+namespace otb
+{
+namespace Wrapper
+{
+template< class TInput, class TOutput>
+class ToComplexPixel
+{
+public:
+  ToComplexPixel ( ) { };
+
+  ~ToComplexPixel( ) { };
+
+  bool operator!=( const ToComplexPixel & ) const
+    {
+    return false;
+    }
+  bool operator==( const ToComplexPixel & other ) const
+    {
+    return !(*this != other);
+    }
+  inline TOutput operator( )( const TInput & A ) const
+    {
+      return TOutput( static_cast<typename TOutput::value_type> ( A[0] ),
+		      static_cast<typename TOutput::value_type> ( A[1] ) );
+    }
+};
+
+class DomainTransform : public Application
+{
+public:
+  /** Standard class typedefs. */
+  typedef DomainTransform     Self;
+  typedef Application                   Superclass;
+  typedef itk::SmartPointer<Self>       Pointer;
+
+  typedef itk::SmartPointer<const Self> ConstPointer;
+  typedef float InputPixelType;
+  typedef float OutputPixelType;
+
+  /** Standard macro */
+  itkNewMacro(Self);
+
+  itkTypeMacro(Self, otb::Application);
+
+private:
+  DomainTransform() {}
+
+  ~DomainTransform() ITK_OVERRIDE
+    {
+    // This is a trick to make sure fftw will cleanup its threads when application
+    // shared lib is released.
+    #if defined(ITK_USE_FFTWF)
+      fftwf_cleanup_threads();
+      fftwf_cleanup();
+    #endif
+    #if defined(ITK_USE_FFTWD)
+      fftw_cleanup_threads();
+      fftw_cleanup();
+    #endif
+    }
+
+  void DoInit() ITK_OVERRIDE
+    {
+    SetName("DomainTransform");
+    SetDescription("Domain Transform application for wavelet and fourier");
+
+    // Documentation
+    SetDocName("DomainTransform");
+    SetDocLongDescription("Domain Transform application for wavelet and fourier");
+    SetDocLimitations("This application is not streamed, check your system resources when processing large images");
+    SetDocAuthors("OTB-Team");
+    SetDocSeeAlso("otbWaveletImageFilter, otbWaveletInverseImageFilter, otbWaveletTransform");
+    AddDocTag(Tags::Filter);
+
+    AddParameter(ParameterType_InputImage, "in",  "Input Image");
+    SetParameterDescription("in", "This will take an input image to be transformed"
+      " image. For FFT inverse transform, it expects a complex image as two-band"
+      " image in which first band represent real part and second band represent"
+      " imaginary part.");
+
+    AddRAMParameter();
+
+    AddParameter(ParameterType_Choice, "mode", "mode");
+    SetParameterDescription("mode", "This parameter allows one to select between"
+      " fft(fourier) and wavelet");
+    AddChoice("mode.fft", "FFT transform");
+    SetParameterDescription("mode.fft", "FFT transform");
+    AddParameter(ParameterType_Empty, "mode.fft.shift", "false");
+    SetParameterDescription("mode.fft.shift", "Shift transform of fft filter");
+    AddChoice("mode.wavelet", "wavelet");
+    SetParameterDescription("mode.wavelet", "Wavelet transform");
+    AddParameter(ParameterType_Choice, "mode.wavelet.form", "Select wavelet form");
+
+    AddChoice("mode.wavelet.form.haar", "HAAR");
+    AddChoice("mode.wavelet.form.db4", "DAUBECHIES4");
+    AddChoice("mode.wavelet.form.db6", "DAUBECHIES6");
+    AddChoice("mode.wavelet.form.db8", "DAUBECHIES8");
+    AddChoice("mode.wavelet.form.db12", "DAUBECHIES12");
+    AddChoice("mode.wavelet.form.db20", "DAUBECHIES20");
+    AddChoice("mode.wavelet.form.sb24", "SPLINE_BIORTHOGONAL_2_4");
+    AddChoice("mode.wavelet.form.sb44", "SPLINE_BIORTHOGONAL_4_4");
+    AddChoice("mode.wavelet.form.sym8", "SYMLET8");
+
+    //Default value
+    SetParameterString("mode", "wavelet");
+    SetParameterString("mode.wavelet.form", "haar");
+
+    AddParameter(ParameterType_Int,"mode.wavelet.nlevels","Number of decomposition levels");
+    SetParameterDescription("mode.wavelet.nlevels","Number of decomposition levels");
+    SetDefaultParameterInt("mode.wavelet.nlevels",2);
+    SetMinimumParameterIntValue("mode.wavelet.nlevels",2);
+
+    AddParameter(ParameterType_Choice,"dir", "Direction");
+
+    AddChoice("dir.fwd", "Forward");
+    AddChoice("dir.inv", "Inverse");
+
+    AddParameter(ParameterType_OutputImage, "out", "Output Image");
+    SetParameterDescription("out", "This parameter holds the output file name to"
+      " which transformed image will be written. This has a slightly different"
+      " behaviour depending on transform type. \n For Wavelet, output is a"
+      " single band image for both forward and inverse transform. \n For FFT"
+      " forward transform, output is two band image where first band represents"
+      " real part and second band represents imaginary part of a complex image.");
+
+    SetDocExampleParameterValue("in", "input.tif");
+    SetDocExampleParameterValue("mode.wavelet.form", "haar");
+    SetDocExampleParameterValue("out", "output_wavelet_haar.tif");
+    }
+
+  void DoUpdateParameters() ITK_OVERRIDE
+    {
+
+    }
+
+  void DoExecute() ITK_OVERRIDE
+    {
+    int dir = GetParameterInt("dir");
+    int mode = GetParameterInt("mode");
+
+    if( dir != 0 && dir != 1)
+      {
+      itkExceptionMacro(<< "-dir is '"
+         << dir << "'."
+         << "It must be either 'fwd' or 'inv'");
+      }
+
+    if( mode != 0 && mode != 1)
+      {
+      itkExceptionMacro( << "mode is '"
+         << mode << "'."
+         << "It must must be either 'fft' or 'wavelet'");
+      }
+
+    if ( mode == 1)
+      {
+      int wavelet_type = GetParameterInt("mode.wavelet.form");
+      unsigned int nlevels = GetParameterInt("mode.wavelet.nlevels");
+      switch (wavelet_type)
+        {
+        case 0:
+          {
+          DoWaveletTransform<otb::Wavelet::HAAR>(dir, nlevels);
+          break;
+          }
+        case 1:
+          {
+          DoWaveletTransform<otb::Wavelet::DB4>(dir, nlevels);
+          break;
+          }
+        case 2:
+          {
+          DoWaveletTransform<otb::Wavelet::DB4>(dir, nlevels);
+          break;
+          }
+        case 3:
+          {
+          DoWaveletTransform<otb::Wavelet::DB6>(dir, nlevels);
+          break;
+          }
+        case 4:
+          {
+          DoWaveletTransform<otb::Wavelet::DB8>(dir, nlevels);
+          break;
+          }
+        case 5:
+          {
+          DoWaveletTransform<otb::Wavelet::DB12>(dir, nlevels);
+          break;
+          }
+        case 6:
+          {
+          DoWaveletTransform<otb::Wavelet::DB20>(dir, nlevels);
+          break;
+          }
+        case 7:
+          {
+          DoWaveletTransform<otb::Wavelet::SPLINE_BIORTHOGONAL_2_4 >(dir, nlevels);
+          break;
+          }
+        case 8:
+          {
+          DoWaveletTransform<otb::Wavelet::SPLINE_BIORTHOGONAL_4_4>(dir, nlevels);
+          break;
+          }
+        case 9:
+          {
+          DoWaveletTransform<otb::Wavelet::SYMLET8>(dir, nlevels);
+          break;
+          }
+        default:
+          {
+          itkExceptionMacro( << "Invalid wavelet type: '" <<  wavelet_type << "'");
+          break;
+          }
+        }
+      }
+    else
+      {
+      // fft ttransform
+      bool shift = IsParameterEnabled( "mode.fft.shift");
+      typedef otb::Image< std::complex<OutputPixelType> >          ComplexOutputImageType;
+      
+      if (dir == 0 )
+        {
+        //forward fft
+        typedef otb::Image<InputPixelType>          TInputImage;
+        typedef TInputImage::Pointer TInputImagePointer;
+
+        //get input paramter as otb::Image<InputPixelType>
+        TInputImagePointer inImage = GetParameterImage<TInputImage>("in");
+
+        //typedef itk::::ForwardFFTImageFilter over otbImage< InputPixelType >
+
+        typedef itk::ForwardFFTImageFilter < TInputImage, ComplexOutputImageType > FFTFilter;
+        FFTFilter::Pointer fwdFilter = FFTFilter::New();
+        fwdFilter->SetInput( inImage );
+
+	
+        //typedef VectorImage for output of UnaryFunctorImageFilter
+        typedef otb::VectorImage<OutputPixelType>          TOutputImage;
+
+	typedef otb::ComplexToVectorImageCastFilter<
+	  ComplexOutputImageType,
+	  TOutputImage > ComplexToVectorImageCastFilter;
+	ComplexToVectorImageCastFilter::Pointer unaryFunctorImageFilter = ComplexToVectorImageCastFilter::New();
+	
+        if( shift)
+          {
+          otbAppLogINFO( << "Applying Shift image filter" );
+          typedef itk::FFTShiftImageFilter<
+            typename FFTFilter::OutputImageType,
+            typename FFTFilter::OutputImageType > FFTShiftFilterType;
+
+          FFTShiftFilterType::Pointer
+            fftShiftFilter = FFTShiftFilterType::New();
+
+          fftShiftFilter->SetInput( fwdFilter->GetOutput() );
+          fftShiftFilter->Update();
+          unaryFunctorImageFilter->SetInput(fftShiftFilter->GetOutput() );
+          }
+        else
+          {
+          unaryFunctorImageFilter->SetInput(fwdFilter->GetOutput());
+          }
+
+        unaryFunctorImageFilter->Update();
+
+        //set output image
+        SetParameterOutputImage<TOutputImage>("out",
+          unaryFunctorImageFilter->GetOutput());
+        }
+      else
+        {
+        //inverse fft
+        typedef otb::VectorImage<InputPixelType>          TInputImage;
+        typedef TInputImage::Pointer TInputImagePointer;
+
+        TInputImagePointer inImage = GetParameterImage("in");
+
+        // typedef TComplexImage for InverseFFTImageFilter input
+        // This a image type of std::complex<InputPixelType>
+        typedef otb::Image<
+          std::complex<InputPixelType>, 2 > TComplexImage;
+        //typedef TOutputImage for InverseFFTImageFilter output
+        typedef otb::Image< OutputPixelType >  TOutputImage;
+	
+        // a unary functor to convert vectorimage to complex image
+        typedef itk::UnaryFunctorImageFilter
+          <TInputImage,
+           TComplexImage,
+           ToComplexPixel
+           <TInputImage::PixelType,
+            TComplexImage::PixelType> > UnaryFunctorImageFilter;
+
+        UnaryFunctorImageFilter::Pointer
+          unaryFunctorImageFilter = UnaryFunctorImageFilter::New();
+
+        if( shift)
+          {
+          typedef itk::FFTShiftImageFilter<
+            TInputImage,
+            TInputImage > FFTShiftFilterType;
+
+          FFTShiftFilterType::Pointer
+            fftShiftFilter = FFTShiftFilterType::New();
+          fftShiftFilter->SetInput( inImage );
+
+          fftShiftFilter->Update();
+
+          unaryFunctorImageFilter->SetInput(fftShiftFilter->GetOutput() );
+          }
+        else
+          {
+          unaryFunctorImageFilter->SetInput(inImage);
+          }
+
+        unaryFunctorImageFilter->Update();
+
+        //typedef itk::::InverseFFTImageFilter over TComplexImage
+        typedef itk::InverseFFTImageFilter
+          < TComplexImage,
+            TOutputImage > FFTFilter;
+        FFTFilter::Pointer invFilter = FFTFilter::New();
+        invFilter->SetInput( unaryFunctorImageFilter->GetOutput() );
+        invFilter->Update();
+
+        //set output image
+        SetParameterOutputImage<TOutputImage>("out", invFilter->GetOutput());
+        }
+      }
+    }
+
+  template<otb::Wavelet::Wavelet TWaveletOperator>
+  void DoWaveletTransform(const int dir,
+                          const unsigned int nlevels,
+                          const std::string inkey = "in",
+                          const std::string outkey = "out")
+    {
+    typedef otb::Image< InputPixelType >  TInputImage;
+    typedef otb::Image< OutputPixelType >  TOutputImage;
+    typedef typename TInputImage::Pointer TInputImagePointer;
+
+    TInputImagePointer inImage = GetParameterImage<TInputImage>(inkey);
+
+    if( dir == 0)
+      {
+      typedef otb::WaveletImageFilter<
+        TInputImage,
+        TOutputImage,
+        TWaveletOperator> TWaveletImageFilter;
+
+      typedef typename
+        TWaveletImageFilter::Pointer
+        TWaveletImageFilterPointer;
+      TWaveletImageFilterPointer waveletImageFilter =
+        TWaveletImageFilter::New();
+
+      waveletImageFilter->SetInput(inImage);
+      waveletImageFilter->SetNumberOfDecompositions(nlevels);
+      waveletImageFilter->Update();
+      SetParameterOutputImage<TOutputImage>(outkey, waveletImageFilter->GetOutput() );
+      }
+    else
+      {
+      typedef otb::WaveletInverseImageFilter<
+        TInputImage,
+        TOutputImage,
+        TWaveletOperator > TWaveletImageFilter;
+
+      typedef typename
+        TWaveletImageFilter::Pointer
+        TWaveletImageFilterPointer;
+
+      TWaveletImageFilterPointer waveletImageFilter =
+        TWaveletImageFilter::New();
+
+      waveletImageFilter->SetInput(inImage);
+      waveletImageFilter->SetNumberOfDecompositions(nlevels);
+      waveletImageFilter->Update();
+
+      SetParameterOutputImage<TOutputImage>( outkey, waveletImageFilter->GetOutput() );
+      }
+    }
+};
+
+} // end of namespace Wrapper
+} // end of namespace otb
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::DomainTransform)

Modules/Applications/AppDomainTransform/otb-module.cmake

+set(DOCUMENTATION "Application to perform Wavelet and Fourier transform. This application is not streamed")
+
+otb_module(OTBAppDomainTransform
+DEPENDS
+  OTBWavelet
+  OTBApplicationEngine
+
+  TEST_DEPENDS
+    OTBTestKernel
+    OTBCommandLine
+
+  DESCRIPTION
+    "${DOCUMENTATION}"
+)
+

Modules/Applications/AppDomainTransform/test/CMakeLists.txt

+otb_module_test()
+otb_test_application(NAME apTvDomainTransform_wav_db20_fwd
+  APP  DomainTransform
+  OPTIONS -in ${INPUTDATA}/QB_Toulouse_Ortho_PAN.tif
+  -mode wavelet
+  -mode.wavelet.form db20
+  -mode.wavelet.nlevels 2
+  -dir fwd
+  -out ${TEMP}/apTvDomainTransform_wav_db20_fwd.tif
+  )
+
+otb_test_application(NAME apTvDomainTransform_wav_haar_inv
+  APP  DomainTransform
+  OPTIONS -in ${INPUTDATA}/QB_Toulouse_Ortho_PAN.tif
+  -mode wavelet
+  -mode.wavelet.form haar
+  -dir inv
+  -out ${TEMP}/apTvDomainTransform_wav_haar_inv.tif
+  )
+
+otb_test_application(NAME apTvDomainTransform_fft_shift_fwd
+  APP  DomainTransform
+  OPTIONS -in ${INPUTDATA}/QB_Toulouse_Ortho_PAN.tif
+  -mode fft
+  -mode.fft.shift true
+  -dir fwd
+  -out ${TEMP}/apTvDomainTransform_fft_shift_fwd.tif
+  )
+		   
+
+otb_test_application(NAME apTvDomainTransform_fft_inv
+  APP  DomainTransform
+  OPTIONS -in ${INPUTDATA}/QB_Toulouse_Ortho_PAN.tif
+  -mode fft
+  -dir inv
+  -out ${TEMP}/apTvDomainTransform_fft_shift_inv.tif
+  )
+		   

Modules/Filtering/Wavelet/include/otbSubsampleImageFilter.h

@@ -102,7 +102,6 @@ protected:
   SubsampleImageFilter ()
     {
     m_SubsampleFactor.Fill(1);
-    //this->SetNumberOfThreads(1);
     }
   ~SubsampleImageFilter() ITK_OVERRIDE {}
 

Modules/Filtering/Wavelet/include/otbSubsampleImageFilter.txx

@@ -83,7 +83,7 @@ SubsampleImageFilter<TInputImage, TOutputImage, TDirectionOfTransformation>
 {
   Superclass::CallCopyOutputRegionToInputRegion(destRegion, srcRegion);
 
-  if (static_cast<int>(DirectionOfTransformation) == Wavelet::INVERSE)
+  if (DirectionOfTransformation == Wavelet::INVERSE)
     {
     typename OutputImageRegionType::IndexType srcIndex = srcRegion.GetIndex();
     typename OutputImageRegionType::SizeType  srcSize = srcRegion.GetSize();
@@ -112,7 +112,7 @@ SubsampleImageFilter<TInputImage, TOutputImage, TDirectionOfTransformation>
 {
   Superclass::CallCopyInputRegionToOutputRegion(destRegion, srcRegion);
 
-  if (static_cast<int>(DirectionOfTransformation) == Wavelet::INVERSE)
+  if (DirectionOfTransformation == Wavelet::INVERSE)
     {
     typename InputImageRegionType::IndexType srcIndex = srcRegion.GetIndex();
     typename InputImageRegionType::SizeType  srcSize = srcRegion.GetSize();
@@ -161,8 +161,7 @@ SubsampleImageFilter<TInputImage, TOutputImage, TDirectionOfTransformation>
   SubsampledImageRegionConstIterator<InputImageType> inputIter
     (this->GetInput(), inputRegionForThread);
 
-  if (static_cast<int>(itkGetStaticConstMacro(DirectionOfTransformation))
-        == static_cast<int>(Wavelet::FORWARD))
+  if (DirectionOfTransformation  == Wavelet::FORWARD)
     {
     inputIter.SetSubsampleFactor(GetSubsampleFactor());
     inputIter.GoToBegin();

Modules/Filtering/Wavelet/include/otbWaveletFilterBank.h

@@ -399,6 +399,11 @@ protected:
   WaveletFilterBank();
   ~WaveletFilterBank() ITK_OVERRIDE {}
 
+  void VerifyInputInformation() ITK_OVERRIDE
+  {
+
+  }
+
   /** GenerateOutputInformation
     * Set the size of the output image depending on the decimation factor
     * Copy information from the input image if existing.

Modules/Filtering/Wavelet/include/otbWaveletFilterBank.txx

@@ -27,9 +27,6 @@
 
 #include "itkPeriodicBoundaryCondition.h"
 
-// FIXME
-#define __myDebug__ 0
-
 namespace otb {
 
 /**
@@ -54,7 +51,6 @@ WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
   m_UpSampleFilterFactor = 0;
   m_SubsampleImageFactor = 1;
 
-  //this->SetNumberOfThreads(1);
 }
 
 template <class TInputImage, class TOutputImage, class TWaveletOperator>
@@ -66,11 +62,9 @@ WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
 
   if (GetSubsampleImageFactor() == 1) return;
 
-#if __myDebug__
   otbGenericMsgDebugMacro(<< " down sampling output regions by a factor of " << GetSubsampleImageFactor());
   otbGenericMsgDebugMacro(<< "initial region    " << this->GetInput()->GetLargestPossibleRegion().GetSize()[0]
                           << "," << this->GetInput()->GetLargestPossibleRegion().GetSize()[1]);
-#endif
 
   OutputImageRegionType newRegion;
   this->CallCopyInputRegionToOutputRegion(newRegion, this->GetInput()->GetLargestPossibleRegion());
@@ -80,9 +74,8 @@ WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
     this->GetOutput(i)->SetRegions(newRegion);
     }
 
-#if __myDebug__
   otbGenericMsgDebugMacro(<< "new region output " << newRegion.GetSize()[0] << "," << newRegion.GetSize()[1]);
-#endif
+
 }
 
 template <class TInputImage, class TOutputImage, class TWaveletOperator>
@@ -185,12 +178,13 @@ WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::FORWARD>
     this->CallCopyInputRegionToOutputRegion(InputImageDimension - 1 - direction,
                                             smallerRegion, largerRegion);
 
-    for (unsigned int i = 0; i < m_InternalImages[direction].size(); ++i)
+    const unsigned int d = InputImageDimension - 2 - direction;
+    for (unsigned int i = 0; i < m_InternalImages[d].size(); ++i)
       {
-      m_InternalImages[InputImageDimension - 2 - direction][i] = OutputImageType::New();
-      m_InternalImages[InputImageDimension - 2 - direction][i]->SetRegions(smallerRegion);
-      m_InternalImages[InputImageDimension - 2 - direction][i]->Allocate();
-      m_InternalImages[InputImageDimension - 2 - direction][i]->FillBuffer(0);
+      m_InternalImages[d][i] = OutputImageType::New();
+      m_InternalImages[d][i]->SetRegions(smallerRegion);
+      m_InternalImages[d][i]->Allocate();
+      m_InternalImages[d][i]->FillBuffer(0);
       }
 
     largerRegion = smallerRegion;
@@ -642,21 +636,18 @@ WaveletFilterBank<TInputImage, TOutputImage, TWaveletOperator, Wavelet::INVERSE>
       }
     }