ENH: New dimensionnalityReduction Application.

Applications/DimensionnalityReduction/CMakeLists.txt

@@ -2,5 +2,9 @@
 OTB_CREATE_APPLICATION(NAME MaximumAutocorrelationFactor
                        SOURCES otbMaximumAutocorrelationFactor.cxx
                        LINK_LIBRARIES OTBIO;OTBCommon;OTBBasicFilters)
+                       
+OTB_CREATE_APPLICATION(NAME DimensionnalityReduction
+                       SOURCES otbDimensionnalityReduction.cxx
+                       LINK_LIBRARIES OTBIO;OTBCommon;OTBBasicFilters)
 
 

Applications/DimensionnalityReduction/otbDimensionnalityReduction.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 "itkImageToImageFilter.h"
+
+#include "otbMaximumAutocorrelationFactorImageFilter.h"
+#include "otbPCAImageFilter.h"
+#include "otbNAPCAImageFilter.h"
+#include "otbLocalActivityVectorImageFilter.h"
+#include "otbMaximumAutocorrelationFactorImageFilter.h"
+#include "otbFastICAImageFilter.h"
+#include "otbFastICAInternalOptimizerVectorImageFilter.h"
+
+//#include "otbVirtualDimensionality.h"
+
+#include "otbStreamingMinMaxVectorImageFilter.h"
+#include "otbVectorRescaleIntensityImageFilter.h"
+
+namespace otb
+{
+namespace Wrapper
+{
+
+class DimensionnalityReduction: public Application
+{
+public:
+  /** Standard class typedefs. */
+  typedef DimensionnalityReduction Self;
+  typedef Application Superclass;
+  typedef itk::SmartPointer<Self> Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  // Dimensionnality reduction typedef
+  typedef otb::MaximumAutocorrelationFactorImageFilter<FloatVectorImageType, FloatVectorImageType> MAFFilterType;
+
+  typedef itk::ImageToImageFilter<FloatVectorImageType, FloatVectorImageType> DimensionnalityReductionFilter;
+
+  // Reduction dimensio filters
+  typedef otb::PCAImageFilter<FloatVectorImageType, FloatVectorImageType, otb::Transform::FORWARD> PCAForwardFilterType;
+  typedef otb::PCAImageFilter<FloatVectorImageType, FloatVectorImageType, otb::Transform::INVERSE> PCAInverseFilterType;
+  //typedef otb::PCAImageFilter< FloatVectorImageType, FloatVectorImageType, otb::Transform::FORWARD >
+
+  typedef otb::LocalActivityVectorImageFilter<FloatVectorImageType, FloatVectorImageType> NoiseFilterType;
+
+  typedef otb::NAPCAImageFilter<FloatVectorImageType, FloatVectorImageType, NoiseFilterType, otb::Transform::FORWARD>
+      NAPCAForwardFilterType;
+  typedef otb::NAPCAImageFilter<FloatVectorImageType, FloatVectorImageType, NoiseFilterType, otb::Transform::INVERSE>
+      NAPCAInverseFilterType;
+
+  typedef otb::MaximumAutocorrelationFactorImageFilter<FloatVectorImageType, FloatVectorImageType> MAFForwardFilterType;
+
+  typedef otb::FastICAImageFilter<FloatVectorImageType, FloatVectorImageType, otb::Transform::FORWARD>
+      ICAForwardFilterType;
+  typedef otb::FastICAImageFilter<FloatVectorImageType, FloatVectorImageType, otb::Transform::INVERSE>
+      ICAInverseFilterType;
+
+  //typedef otb::StreamingStatisticsVectorImageFilter<FloatVectorImageType> StreamingStatisticsVectorImageFilterType;
+  //typedef otb::VirtualDimensionality<double> VDFilterType;
+
+
+  // output rescale
+  typedef otb::StreamingMinMaxVectorImageFilter<FloatVectorImageType> MinMaxFilterType;
+  typedef otb::VectorRescaleIntensityImageFilter<FloatVectorImageType> RescaleImageFilterType;
+
+  /** Standard macro */
+  itkNewMacro(Self)
+  ;
+
+  itkTypeMacro(DimensionnalityReduction, otb::Wrapper::Application)
+  ;
+
+private:
+  void DoInit()
+  {
+    SetName("DimensionnalityReduction");
+    SetDescription("Perform Dimension reduction of the input image.");
+    SetDocName("Dimensionnality reduction application");
+    SetDocLongDescription("Performs dimensionnality reduction on input image. PCA,NA-PCA,MAF,ICA methods are available.");
+    SetDocLimitations(
+                      "Though the inverse transform can be computed, this application only provides the forward transform for now.");
+    SetDocAuthors("OTB-Team");
+    SetDocSeeAlso(
+                  "\"Kernel maximum autocorrelation factor and minimum noise fraction transformations,\" IEEE Transactions on Image Processing, vol. 20, no. 3, pp. 612-624, (2011)");
+
+    AddDocTag(Tags::DimensionReduction);
+    AddDocTag(Tags::Filter);
+
+    AddParameter(ParameterType_InputImage, "in", "Input Image");
+    SetParameterDescription("in", "The input image to apply dimensionnality reduction.");
+    AddParameter(ParameterType_OutputImage, "out", "Output Image");
+    SetParameterDescription("out", "output image. Components are ordered by decreasing eigenvalues.");
+    AddParameter(ParameterType_Group, "rescale", "Rescale Output.");
+
+    MandatoryOff("rescale");
+    //  AddChoice("rescale.no","No rescale");
+    // AddChoice("rescale.minmax","rescale to min max value");
+
+    AddParameter(ParameterType_Float, "rescale.outmin", "Output min value");
+    AddParameter(ParameterType_Float, "rescale.outmax", "Output max value");
+    SetDefaultParameterFloat("rescale.outmin", 0.0);
+    SetParameterDescription("rescale.outmin", "Minimum value of the output image.");
+    SetDefaultParameterFloat("rescale.outmax", 255.0);
+    SetParameterDescription("rescale.outmax", "Maximum value of the output image.");
+
+    AddParameter(ParameterType_OutputImage, "outinv", " Inverse Output Image");
+    SetParameterDescription("outinv", "reconstruct output image.");
+    MandatoryOff("outinv");
+
+    AddParameter(ParameterType_Choice, "method", "Algorithm");
+    SetParameterDescription("method", "Selection of the reduction dimension method.");
+
+    AddChoice("method.pca", "PCA");
+    SetParameterDescription("method.pca", "Principal Component Analysis.");
+    AddChoice("method.napca", "NA-PCA");
+    SetParameterDescription("method.napca", "Noise Adjusted Principal Component Analysis.");
+    AddParameter(ParameterType_Int, "method.napca.radiusx", "Set the x radius of the sliding window.");
+    SetMinimumParameterIntValue("method.napca.radiusx", 1);
+    SetDefaultParameterInt("method.napca.radiusx", 1);
+    AddParameter(ParameterType_Int, "method.napca.radiusy", "Set the y radius of the sliding window.");
+    SetMinimumParameterIntValue("method.napca.radiusy", 1);
+    SetDefaultParameterInt("method.napca.radiusy", 1);
+
+    AddChoice("method.maf", "MAF");
+    SetParameterDescription("method.maf", "Maximum Autocorrelation Factor.");
+    AddChoice("method.ica", "ICA");
+    SetParameterDescription("method.ica", "Independant Component Analysis.");
+    AddParameter(ParameterType_Int, "method.ica.iter", "number of iterations ");
+    SetMinimumParameterIntValue("method.ica.iter", 1);
+    SetDefaultParameterInt("method.ica.iter", 20);
+    MandatoryOff("method.ica.iter");
+
+    AddParameter(ParameterType_Float, "method.ica.mu", "Give the increment weight of W in [0, 1]");
+    SetMinimumParameterFloatValue("method.ica.mu", 0.);
+    SetMaximumParameterFloatValue("method.ica.mu", 1.);
+    SetDefaultParameterFloat("method.ica.mu", 1.);
+    MandatoryOff("method.ica.mu");
+
+    //AddChoice("method.vd","virual Dimension");
+    //SetParameterDescription("method.vd","Virtual Dimension.");
+    //MandatoryOff("method");
+
+    AddParameter(ParameterType_Int, "nbcomp", "Number of Components.");
+    SetParameterDescription("nbcomp", "");
+    SetDefaultParameterInt("nbcomp", 0);
+    MandatoryOff("nbcomp");
+    SetMinimumParameterIntValue("nbcomp", 0);
+
+    AddParameter(ParameterType_Empty, "normalize", "Normalize.");
+    SetParameterDescription("normalize", "center AND reduce data before Dimensionnality reduction.");
+    MandatoryOff("normalize");
+
+    // Doc example parameter settings
+    SetDocExampleParameterValue("in", "cupriteSubHsi.tif");
+    SetDocExampleParameterValue("out", "FilterOutput.tif");
+    SetDocExampleParameterValue("method", "pca");
+  }
+
+  void DoUpdateParameters()
+  {
+  }
+
+  void DoExecute()
+  {
+
+    // Get Parameters
+    int nbComp = GetParameterInt("nbcomp");
+    bool normalize = HasValue("normalize");
+    bool rescale = IsParameterEnabled("rescale");
+
+    bool invTransform = HasValue("outinv");
+    switch (GetParameterInt("method"))
+      {
+      // PCA Algorithm
+      case 0:
+        {
+
+        otbAppLogDEBUG( << "PCA Algorithm ");
+        PCAForwardFilterType::Pointer filter = PCAForwardFilterType::New();
+        m_ForwardFilter = filter;
+        PCAInverseFilterType::Pointer invFilter = PCAInverseFilterType::New();
+        m_InverseFilter = invFilter;
+
+        filter->SetInput(GetParameterFloatVectorImage("in"));
+        filter->SetNumberOfPrincipalComponentsRequired(nbComp);
+        filter->SetUseNormalization(normalize);
+        m_ForwardFilter->Update();
+        if (invTransform)
+          {
+          invFilter->SetInput(m_ForwardFilter->GetOutput());
+          if (normalize)
+            {
+            otbAppLogINFO( << "Normalization MeanValue :"<<filter->GetMeanValues()<<
+                "StdValue :" <<filter->GetStdDevValues() );
+            invFilter->SetMeanValues(filter->GetMeanValues());
+            invFilter->SetStdDevValues(filter->GetStdDevValues());
+            }
+
+          invFilter->SetTransformationMatrix(filter->GetTransformationMatrix());
+          }
+
+        break;
+        }
+      case 1:
+        {
+        otbAppLogDEBUG( << "NA-PCA Algorithm ");
+
+        // NA-PCA
+
+        unsigned int radiusX = static_cast<unsigned int> (GetParameterInt("method.napca.radiusx"));
+        unsigned int radiusY = static_cast<unsigned int> (GetParameterInt("method.napca.radiusy"));
+
+        // Noise filtering
+        NoiseFilterType::RadiusType radius = { { radiusX, radiusY } };
+
+        NAPCAForwardFilterType::Pointer filter = NAPCAForwardFilterType::New();
+        m_ForwardFilter = filter;
+        NAPCAInverseFilterType::Pointer invFilter = NAPCAInverseFilterType::New();
+        m_InverseFilter = invFilter;
+
+        filter->SetInput(GetParameterFloatVectorImage("in"));
+        filter->SetNumberOfPrincipalComponentsRequired(nbComp);
+        filter->SetUseNormalization(normalize);
+        filter->GetNoiseImageFilter()->SetRadius(radius);
+        m_ForwardFilter->Update();
+        if (invTransform)
+          {
+          otbAppLogDEBUG( << "Compute Inverse Transform");
+          invFilter->SetInput(m_ForwardFilter->GetOutput());
+          invFilter->SetMeanValues(filter->GetMeanValues());
+          if (normalize)
+            {
+            invFilter->SetStdDevValues(filter->GetStdDevValues());
+            }
+
+          invFilter->SetTransformationMatrix(filter->GetTransformationMatrix());
+
+          }
+
+        break;
+        }
+      case 2:
+        {
+        otbAppLogDEBUG( << "MAF Algorithm ");
+        MAFForwardFilterType::Pointer filter = MAFForwardFilterType::New();
+        m_ForwardFilter = filter;
+        filter->SetInput(GetParameterFloatVectorImage("in"));
+        m_ForwardFilter->Update();
+
+        otbAppLogINFO( << "V :"<<std::endl<<filter->GetV()<<"Auto-Correlation :"<<std::endl <<filter->GetAutoCorrelation() );
+
+        break;
+        }
+      case 3:
+        {
+        otbAppLogDEBUG( << "ICA Algorithm ");
+
+        unsigned int nbIterations = static_cast<unsigned int> (GetParameterInt("method.ica.iter"));
+        double mu = static_cast<double> (GetParameterFloat("method.ica.mu"));
+
+        ICAForwardFilterType::Pointer filter = ICAForwardFilterType::New();
+        m_ForwardFilter = filter;
+        ICAInverseFilterType::Pointer invFilter = ICAInverseFilterType::New();
+        m_InverseFilter = invFilter;
+
+        filter->SetInput(GetParameterFloatVectorImage("in"));
+        filter->SetNumberOfPrincipalComponentsRequired(nbComp);
+        filter->SetNumberOfIterations(nbIterations);
+        filter->SetMu(mu);
+        m_ForwardFilter->Update();
+        if (invTransform)
+          {
+          otbAppLogDEBUG( << "Compute Inverse Transform");
+          invFilter->SetInput(m_ForwardFilter->GetOutput());
+
+          if (normalize)
+            {
+            invFilter->SetMeanValues(filter->GetMeanValues());
+            invFilter->SetStdDevValues(filter->GetStdDevValues());
+            }
+          invFilter->SetPCATransformationMatrix(filter->GetPCATransformationMatrix());
+          invFilter->SetTransformationMatrix(filter->GetTransformationMatrix());
+          }
+        break;
+        }
+        /* case 4:
+         {
+         otbAppLogDEBUG( << "VD Algorithm");
+
+         break;
+         }*/
+
+      default:
+        {
+        otbAppLogFATAL(<<"non defined method "<<GetParameterInt("method")<<std::endl);
+        break;
+
+        }
+        return;
+      }
+
+    if (invTransform)
+      {
+      if (GetParameterInt("method") == 2) //MAF or VD
+        {
+        otbAppLogWARNING(<<"This application only provides the forward transform .");
+        }
+      else SetParameterOutputImage("outinv", m_InverseFilter->GetOutput());
+      }
+
+    if (!rescale)
+      {
+      SetParameterOutputImage("out", m_ForwardFilter->GetOutput());
+      }
+    else
+      {
+      otbAppLogDEBUG( << "Rescaling " )
+      otbAppLogDEBUG( << "Starting Min/Max computation" )
+
+      m_MinMaxFilter = MinMaxFilterType::New();
+      m_MinMaxFilter->SetInput(m_ForwardFilter->GetOutput());
+      m_MinMaxFilter->GetStreamer()->SetNumberOfLinesStrippedStreaming(50);
+
+      AddProcess(m_MinMaxFilter->GetStreamer(), "Min/Max computing");
+      m_MinMaxFilter->Update();
+
+      otbAppLogDEBUG( << "Min/Max computation done : min=" << m_MinMaxFilter->GetMinimum()
+          << " max=" << m_MinMaxFilter->GetMaximum() )
+
+      FloatVectorImageType::PixelType inMin, inMax;
+
+      m_RescaleFilter = RescaleImageFilterType::New();
+      m_RescaleFilter->SetInput(m_ForwardFilter->GetOutput());
+      m_RescaleFilter->SetInputMinimum(m_MinMaxFilter->GetMinimum());
+      m_RescaleFilter->SetInputMaximum(m_MinMaxFilter->GetMaximum());
+
+      FloatVectorImageType::PixelType outMin, outMax;
+      outMin.SetSize(m_ForwardFilter->GetOutput()->GetNumberOfComponentsPerPixel());
+      outMax.SetSize(m_ForwardFilter->GetOutput()->GetNumberOfComponentsPerPixel());
+      outMin.Fill(GetParameterFloat("rescale.outmin"));
+      outMax.Fill(GetParameterFloat("rescale.outmax"));
+
+      m_RescaleFilter->SetOutputMinimum(outMin);
+      m_RescaleFilter->SetOutputMaximum(outMax);
+      m_RescaleFilter->UpdateOutputInformation();
+
+      SetParameterOutputImage("out", m_RescaleFilter->GetOutput());
+      }
+
+  }
+
+  MinMaxFilterType::Pointer               m_MinMaxFilter;
+  RescaleImageFilterType::Pointer         m_RescaleFilter;
+  DimensionnalityReductionFilter::Pointer m_ForwardFilter;
+  DimensionnalityReductionFilter::Pointer m_InverseFilter;
+};
+
+}
+}
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::DimensionnalityReduction)