Merge branch 'pantex_texture_extraction' into 'develop'

Pantex texture extraction application

Closes #2022

See merge request orfeotoolbox/otb!704

Modules/Applications/AppTextures/app/CMakeLists.txt

@@ -27,3 +27,8 @@ otb_create_application(
   NAME           SFSTextureExtraction
   SOURCES        otbSFSTextureExtraction.cxx
   LINK_LIBRARIES ${${otb-module}_LIBRARIES})
+
+otb_create_application(
+  NAME           PantexTextureExtraction
+  SOURCES        otbPantexTextureExtraction.cxx
+  LINK_LIBRARIES ${${otb-module}_LIBRARIES})

Modules/Applications/AppTextures/app/otbPantexTextureExtraction.cxx

+/*
+ * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "otbWrapperApplication.h"
+#include "otbWrapperApplicationFactory.h"
+
+#include "otbMultiToMonoChannelExtractROI.h"
+#include "otbScalarImageToPanTexTextureFilter.h"
+#include "otbStreamingMinMaxImageFilter.h"
+
+namespace otb
+{
+namespace Wrapper
+{
+
+class PantexTextureExtraction : public Application
+{
+public:
+  
+  /** @name Standard class typedefs
+   * @{
+   */
+  using Self = PantexTextureExtraction;
+  using Superclass = Application;
+  using Pointer = itk::SmartPointer<Self>;
+  using ConstPointer = itk::SmartPointer<const Self>;
+  /** @} */
+
+  /** @name Standard macro
+   * @{
+   */
+  itkNewMacro(Self);
+  itkTypeMacro(PantexTextureExtraction, otb::Application);
+  /** @} */
+
+  using ExtractorFilterType =
+    otb::MultiToMonoChannelExtractROI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType>;
+  using PanTexTextureFilterType = otb::ScalarImageToPanTexTextureFilter<FloatImageType, FloatImageType>;
+  using MinMaxImageFilterType = otb::StreamingMinMaxImageFilter<FloatImageType>;
+
+private:
+  void DoInit() override
+  {
+    SetName("PantexTextureExtraction");
+    SetDescription("Computes Pantex textural features on the selected channel of the input image");
+
+    SetDocLongDescription(
+        "This application computes a texture-derived built-up presence index (PanTex) from textural"
+    "characteristics of scalar images. This is a contrast textural measure based on co-occurance."
+    );
+
+    SetDocLimitations("None");
+
+    SetDocAuthors("OTB-Team");
+    SetDocSeeAlso(" Pesari, M., A. Gerhardinger, F. Kayitakire. 2008.  A robust built-up area precense"
+      " index by anisotropic rotation-invariant textural measure."
+      " IEEE Journal of selected topics in applied earth observations and remote sensing.Vol1, NO3.");
+
+    AddDocTag(Tags::FeatureExtraction);
+    AddDocTag("Textures");
+    
+    AddParameter(ParameterType_InputImage, "in", "Input Image");
+    SetParameterDescription("in", "The input image to compute the features on.");
+
+    AddParameter(ParameterType_Int, "channel", "Selected Channel");
+    SetParameterDescription("channel", "The selected channel index");
+    SetDefaultParameterInt("channel", 1);
+    SetMinimumParameterIntValue("channel", 1);
+
+    AddParameter(ParameterType_OutputImage, "out", "Output Image");
+    SetParameterDescription("out", "Output image containing the selected texture features.");
+
+    AddParameter(ParameterType_Float, "min", "Image minimum");
+    SetParameterDescription("min", "Input image minimum. If this parameter is not set, the application will compute "
+                            "the minumum of the image.");
+    MandatoryOff("min");
+
+    AddParameter(ParameterType_Float, "max", "Image maximum");
+    SetParameterDescription("max", "Input image maximum. If this parameter is not set, the application will compute "
+                            "the maximum of the image.");
+    MandatoryOff("max");
+
+    AddParameter(ParameterType_Int, "sradx", "Window radius (x direction)");
+    SetParameterDescription("sradx", "Radius of the window on which textures are computed (x direction)");
+    SetMinimumParameterIntValue("sradx", 0);
+    SetDefaultParameterInt("sradx", 4);
+
+    AddParameter(ParameterType_Int, "srady", "Window radius (y direction)");
+    SetParameterDescription("srady", "Radius of the window on which textures are computed (y direction)");
+    SetMinimumParameterIntValue("srady", 0);
+    SetDefaultParameterInt("srady", 4);
+
+    AddParameter(ParameterType_Int, "nbin", "Number of bins per axis for histogram generation");
+    SetParameterDescription("nbin", "Number of bins per axis for histogram generation "
+                                    "(number of gray levels considered in the computation of co-occurance).");
+    SetDefaultParameterInt("nbin", 8);
+
+    AddRAMParameter();
+
+    // Doc example parameter settings
+    SetDocExampleParameterValue("in", "qb_RoadExtract.tif");
+    SetDocExampleParameterValue("channel", "2");
+    SetDocExampleParameterValue("min", "0");
+    SetDocExampleParameterValue("max", "255");
+    SetDocExampleParameterValue("nbin", "8");
+    SetDocExampleParameterValue("srady", "4");
+    SetDocExampleParameterValue("sradx", "4");
+
+    SetOfficialDocLink();
+  }
+  
+  
+  void DoUpdateParameters() override
+  {
+    // Nothing to do here : all parameters are independent
+  }
+
+  void DoExecute() override
+  {
+    auto inImage = GetParameterImage("in");
+    inImage->UpdateOutputInformation();
+
+    if ((unsigned int)GetParameterInt("channel") > inImage->GetNumberOfComponentsPerPixel())
+    {
+      itkExceptionMacro(<< "The specified channel index is invalid.");
+    }
+
+    auto extractorFilter = ExtractorFilterType::New();
+    extractorFilter->SetInput(inImage);
+    extractorFilter->SetStartX(inImage->GetLargestPossibleRegion().GetIndex(0));
+    extractorFilter->SetStartY(inImage->GetLargestPossibleRegion().GetIndex(1));
+    extractorFilter->SetSizeX(inImage->GetLargestPossibleRegion().GetSize(0));
+    extractorFilter->SetSizeY(inImage->GetLargestPossibleRegion().GetSize(1));
+    extractorFilter->SetChannel(GetParameterInt("channel"));
+
+    auto textureFilter = PanTexTextureFilterType::New();
+
+    textureFilter->SetNumberOfBinsPerAxis(GetParameterInt("nbin"));
+    textureFilter->SetRadius( {(unsigned int) GetParameterInt("sradx"),
+                               (unsigned int) GetParameterInt("srady")} );
+
+    // Compute min and max only if one the corresponding parameter has not been set
+    if (!HasValue("min") || !HasValue("max"))
+    {
+      auto minMaxFilter = MinMaxImageFilterType::New();
+      minMaxFilter->SetInput(extractorFilter->GetOutput());
+      minMaxFilter->Update();
+      
+      if (!HasValue("min"))
+        otbAppLogINFO(<< "Computed Minimum: " << minMaxFilter->GetMinimum());
+
+      if (!HasValue("max"))
+        otbAppLogINFO(<< "Computed Maximum: " << minMaxFilter->GetMaximum());
+
+      textureFilter->SetInputImageMinimum(HasValue("min") ? GetParameterFloat("min") 
+                                            : minMaxFilter->GetMinimum());
+      textureFilter->SetInputImageMaximum(HasValue("max") ? GetParameterFloat("max") 
+                                            : minMaxFilter->GetMaximum());
+
+    }
+    else
+    {
+      textureFilter->SetInputImageMinimum(GetParameterFloat("min"));
+      textureFilter->SetInputImageMaximum(GetParameterFloat("max"));
+    }
+    
+    textureFilter->SetInput(extractorFilter->GetOutput());
+    SetParameterOutputImage("out", textureFilter->GetOutput());
+
+    RegisterPipeline();
+  }
+};
+
+}
+}
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::PantexTextureExtraction)

Modules/Applications/AppTextures/otb-module.cmake

@@ -27,6 +27,7 @@ otb_module(OTBAppTextures
     OTBImageBase
     OTBApplicationEngine
     OTBImageManipulation
+    OTBStatistics
     OTBObjectList
 
   TEST_DEPENDS

Modules/Applications/AppTextures/test/CMakeLists.txt

@@ -43,3 +43,17 @@ otb_test_application(NAME  apTvFESFSTextureExtraction
                    			 ${BASELINE}/apTvFESFSTextureExtraction.tif
                  		     ${TEMP}/apTvFESFSTextureExtraction.tif)
 
+
+#----------- PantexTextureExtraction TESTS ----------------
+otb_test_application(NAME  apTvFEPantexTextureExtraction
+                     APP  PantexTextureExtraction
+                     OPTIONS -in ${INPUTDATA}/QB_Toulouse_Ortho_PAN.tif
+                             -channel 1
+                             -out ${TEMP}/apTvFEPantexTextureExtraction.tif
+                             -nbin 8
+                             -sradx 4
+                             -srady 4
+                     VALID   --compare-image ${NOTOL}
+                   			 ${BASELINE}/apTvFEPantexTextureExtraction.tif
+                 		     ${TEMP}/apTvFEPantexTextureExtraction.tif)
+

Modules/Feature/Textures/include/otbGreyLevelCooccurrenceIndexedList.h

@@ -41,7 +41,7 @@ namespace otb
 * iteration over the given input image. This class keep an internal itk::Array
 * as a lookup array with size as [nbbins x nbbins]. The lookup array stores
 * position CooccurrencePairType in the VectorType. It ensures us that all elements
-* in Vector are unqiue in terms of the index value in the pair. For any given
+* in Vector are unique in terms of the index value in the pair. For any given
 * pixel index, -1 value indicates zero existence of the index in the
 * VectorType. This avoid searching all elements in VectorType for each pixel
 * index added during neighborhood iterator. It is also used to decide wheather
@@ -136,9 +136,7 @@ public:
 
 protected:
   GreyLevelCooccurrenceIndexedList();
-  ~GreyLevelCooccurrenceIndexedList() override
-  {
-  }
+  ~GreyLevelCooccurrenceIndexedList() override = default;
 
   /** create a cooccurrence pair with given index and frequency = 1
     * value. Next occurrence of same index is checked via m_LookupArray and the

Modules/Feature/Textures/include/otbScalarImageToPanTexTextureFilter.h

@@ -113,7 +113,7 @@ protected:
   /** Constructor */
   ScalarImageToPanTexTextureFilter();
   /** Destructor */
-  ~ScalarImageToPanTexTextureFilter() override;
+  ~ScalarImageToPanTexTextureFilter() override = default;
   /** Generate the input requested region */
   void GenerateInputRequestedRegion() override;
   /** Parallel textures extraction */

Modules/Feature/Textures/include/otbScalarImageToPanTexTextureFilter.hxx

@@ -37,36 +37,9 @@ ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>::ScalarImageToPanTex
   // There are 1 output corresponding to the Pan Tex texture indice
   this->SetNumberOfRequiredOutputs(1);
 
-  // Fill the offset list for contrast computation
-  OffsetType off;
-  off[0] = 0;
-  off[1] = 1;
-  m_OffsetList.push_back(off); //(0, 1)
-  off[1] = 2;
-  m_OffsetList.push_back(off); //(0, 2)
-  off[0] = 1;
-  off[1] = -2;
-  m_OffsetList.push_back(off); //(1, -2)
-  off[1] = -1;
-  m_OffsetList.push_back(off); //(1, -1)
-  off[1] = 0;
-  m_OffsetList.push_back(off); //(1, 0)
-  off[1] = 1;
-  m_OffsetList.push_back(off); //(1, 1)
-  off[1] = 2;
-  m_OffsetList.push_back(off); //(1, 2)
-  off[0] = 2;
-  off[1] = -1;
-  m_OffsetList.push_back(off); //(2, -1)
-  off[1] = 0;
-  m_OffsetList.push_back(off); //(2, 0)
-  off[1] = 1;
-  m_OffsetList.push_back(off); //(2, 1)
-}
-
-template <class TInputImage, class TOutputImage>
-ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>::~ScalarImageToPanTexTextureFilter()
-{
+  // Ten offsets are selected for contrast computation (2 pixels displacement grid, and given that the
+  // co-occurance matrix is symmetric
+  m_OffsetList = { {0, 1}, {0, 2}, {1, -2}, {1, -1}, {1, 0}, {1, 1}, {1, 2}, {2, -1}, {2, 0}, {2, 1} };
 }
 
 template <class TInputImage, class TOutputImage>
@@ -93,9 +66,7 @@ void ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>::GenerateInputR
   InputRegionType inputRequestedRegion = outputRequestedRegion;
 
   // Apply the radius
-  SizeType maxOffsetSize;
-  maxOffsetSize[0] = 2;
-  maxOffsetSize[1] = 2;
+  SizeType maxOffsetSize = {2, 2};
   inputRequestedRegion.PadByRadius(m_Radius + maxOffsetSize);
 
   // Try to apply the requested region to the input image
@@ -123,8 +94,6 @@ void ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>::ThreadedGenera
   OutputImagePointerType outputPtr = this->GetOutput();
 
   itk::ImageRegionIteratorWithIndex<OutputImageType> outputIt(outputPtr, outputRegionForThread);
-
-  // Go to begin
   outputIt.GoToBegin();
 
   // Set-up progress reporting
@@ -154,13 +123,11 @@ void ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>::ThreadedGenera
         inputIndex[dim] = outputIt.GetIndex()[dim] - m_Radius[dim];
         inputSize[dim]  = 2 * m_Radius[dim] + 1;
       }
-      // Build the input  region
-      InputRegionType inputRegion;
-      inputRegion.SetIndex(inputIndex);
-      inputRegion.SetSize(inputSize);
+      
+      // Build the input region
+      InputRegionType inputRegion = {inputIndex, inputSize};
       inputRegion.Crop(inputPtr->GetRequestedRegion());
 
-
       SizeType neighborhoodRadius;
       /** calculate minimum offset and set it as neighborhood radius **/
       unsigned int minRadius = 0;
@@ -178,8 +145,8 @@ void ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>::ThreadedGenera
       GLCIList->Initialize(m_NumberOfBinsPerAxis, m_InputImageMinimum, m_InputImageMaximum);
 
       typedef itk::ConstNeighborhoodIterator<InputImageType> NeighborhoodIteratorType;
-      NeighborhoodIteratorType                               neighborIt;
-      neighborIt = NeighborhoodIteratorType(neighborhoodRadius, inputPtr, inputRegion);
+      NeighborhoodIteratorType neighborIt = NeighborhoodIteratorType(neighborhoodRadius, inputPtr, inputRegion);
+
       for (neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt)
       {
         const InputPixelType centerPixelIntensity = neighborIt.GetCenterPixel();
@@ -201,8 +168,8 @@ void ScalarImageToPanTexTextureFilter<TInputImage, TOutputImage>::ThreadedGenera
       constVectorIt  = glcVector.begin();
       while (constVectorIt != glcVector.end())
       {
-        CooccurrenceIndexType index     = (*constVectorIt).first;
-        RelativeFrequencyType frequency = (*constVectorIt).second / totalFrequency;
+        CooccurrenceIndexType index     = constVectorIt->first;
+        RelativeFrequencyType frequency = constVectorIt->second / totalFrequency;
         inertia += (index[0] - index[1]) * (index[0] - index[1]) * frequency;
         ++constVectorIt;
       }