UPDATE : New Functor for DEMPolygonsAnalysis

app/CMakeLists.txt

@@ -53,3 +53,8 @@ OTB_CREATE_APPLICATION(NAME SARCoRegistration
                        SOURCES otbSARCoRegistration.cxx
                        LINK_LIBRARIES ${${otb-module}_LIBRARIES}
 )
+
+OTB_CREATE_APPLICATION(NAME SARInterpolation
+                       SOURCES otbSARInterpolation.cxx
+                       LINK_LIBRARIES ${${otb-module}_LIBRARIES}
+)

app/otbSARInterpolation.cxx

+/*
+ * Copyright (C) 2005-2018 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 "otbSARDEMPolygonsAnalysisImageFilter.h"
+#include "otbSARInterpolationFunctor.h"
+
+#include "otbWrapperOutputImageParameter.h"
+
+#include <iostream>
+#include <string>
+#include <fstream>
+
+namespace otb
+{
+namespace Wrapper
+{
+
+class SARInterpolation : public Application
+{
+public:
+  typedef SARInterpolation Self;
+  typedef itk::SmartPointer<Self> Pointer; 
+
+  itkNewMacro(Self);
+  itkTypeMacro(SARInterpolation, otb::Wrapper::Application);
+
+  // Pixels
+  typedef typename FloatVectorImageType::PixelType       ImageInPixelType;
+  typedef typename FloatImageType::PixelType             ImageOutPixelType;   
+
+  // Function
+  typedef otb::Function::SARPolygonsFunctor<ImageInPixelType, ImageOutPixelType> PolygonsFunctorType;
+  typedef otb::Function::SARInterpolationFunctor<ImageInPixelType, ImageInPixelType> InterFunctorType;
+
+  // Filters
+  typedef otb::SARDEMPolygonsAnalysisImageFilter < FloatVectorImageType, FloatVectorImageType, FloatImageType, ComplexFloatImageType, InterFunctorType > FilterType;
+
+private:
+  void DoInit()
+  {
+    SetName("SARInterpolation");
+    SetDescription("SAR Amplitude estimation thanks to the associated DEM.");
+
+    SetDocLongDescription("This application estimate an simulated amplitude image thankds to a DEM file.\n");
+
+    //Optional descriptors
+    SetDocLimitations("Only Sentinel 1 products are supported for now.");
+    SetDocAuthors("OTB-Team");
+    SetDocSeeAlso(" ");
+    AddDocTag(Tags::SAR);
+
+    //Parameter declarations
+    AddParameter(ParameterType_InputImage,  "indemproj",   "Input vector of DEM projected into SAR geometry");
+    SetParameterDescription("indemproj", "Input vector image for amplitude estimation.");
+
+    AddParameter(ParameterType_InputImage,  "indem",   "Input DEM");
+    SetParameterDescription("indem", "DEM to extract DEM geometry.");
+
+    AddParameter(ParameterType_ComplexInputImage,  "insar",   "Input SAR image");
+    SetParameterDescription("insar", "SAR Image to extract SAR geometry.");
+
+    AddParameter(ParameterType_Float,  "ingain",   "Multiplying gain");
+    SetParameterDescription("ingain", "Multiplying gain to obtain a mean radiometry at 100.");
+
+    AddParameter(ParameterType_Int,  "indirectiondemc", "Range direction for DEM scan");
+    SetParameterDescription("indirectiondemc", "Range direction for DEM scan.");
+    SetDefaultParameterInt("indirectiondemc", 1);
+    MandatoryOff("indirectiondemc");
+
+    AddParameter(ParameterType_Int,  "indirectiondeml", "Azimut direction for DEM scan");
+    SetParameterDescription("indirectiondeml", "Azimut direction for DEM scan.");
+    SetDefaultParameterInt("indirectiondeml", 1);
+    MandatoryOff("indirectiondeml");
+
+    AddParameter(ParameterType_OutputImage, "out", "Output amplitude Image");
+    SetParameterDescription("out","Output  amplitude Image.");
+
+    AddRAMParameter();
+
+    SetDocExampleParameterValue("indemproj","CLZY_S21E055.tiff");
+    SetDocExampleParameterValue("indem","S21E055.hgt");
+    SetDocExampleParameterValue("insar","s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002_SLC.tiff");
+    SetDocExampleParameterValue("out","s1a-s4-simu-amplitude.tiff");
+  }
+
+  void DoUpdateParameters()
+  {
+    // Nothing to do here : all parameters are independent
+  }
+
+  void DoExecute()
+  {  
+    // Get numeric parameters : gain and directions for DEM scan
+    double gain = GetParameterFloat("ingain");
+    int DEMScanDirectionC = GetParameterInt("indirectiondemc");
+    int DEMScanDirectionL = GetParameterInt("indirectiondeml");
+
+    otbAppLogINFO(<<"Gain : "<<gain);
+    otbAppLogINFO(<<"Direction (DEM scan) in range : "<<DEMScanDirectionC);
+    otbAppLogINFO(<<"Direction (DEM scan) in azimut : "<<DEMScanDirectionL);
+    
+    // Start the first pipelines (Read SAR and DEM image metedata)
+    ComplexFloatImageType::Pointer SARPtr = GetParameterComplexFloatImage("insar");
+    FloatImageType::Pointer DEMPtr = GetParameterFloatImage("indem");
+
+
+    // Function Type : SARInterpolationFunctor
+    InterFunctorType::Pointer functor =  InterFunctorType::New(1);
+    
+    // Interpolation Filter (use SARDEMPolygonsAnalysisImageFilter with SARInterpolationFunctor
+    // to estimate amplitude image)
+    FilterType::Pointer filterInterpolation = FilterType::New();
+    m_Ref.push_back(filterInterpolation.GetPointer());
+    filterInterpolation->SetSARImageKeyWorList(SARPtr->GetImageKeywordlist());
+    filterInterpolation->SetSARImagePtr(SARPtr);
+    filterInterpolation->SetDEMImagePtr(DEMPtr);
+    filterInterpolation->SetDEMInformation(gain, DEMScanDirectionC, DEMScanDirectionL);
+    filterInterpolation->SetFunctorPtr(functor);
+    filterInterpolation->initializeMarginAndRSTransform();
+
+    // Define the main pipeline (controlled with extended FileName)
+    std::string origin_FileName = GetParameterString("out");
+
+    // Check if FileName is extended (with the ? caracter)
+    // If not extended then override the FileName
+    if (origin_FileName.find("?") == std::string::npos && !origin_FileName.empty()) 
+      {
+	std::string extendedFileName = origin_FileName;
+
+	// Get the ram value (in MB)
+	int ram = GetParameterInt("ram");
+	
+	// Define with the ram value, the number of lines for the streaming
+	int nbColSAR = SARPtr->GetLargestPossibleRegion().GetSize()[0];
+	int nbLinesSAR = SARPtr->GetLargestPossibleRegion().GetSize()[1];
+	// To determine the number of lines : 
+	// nbColSAR * nbLinesStreaming * sizeof(OutputPixel) = RamValue/2 (/2 to be sure)
+	// OutputPixel are float => sizeof(OutputPixel) = 4 (Bytes)
+	long int ram_In_KBytes = (ram/2) * 1024;
+	long int nbLinesStreaming = (ram_In_KBytes / (nbColSAR)) * (1024/4);
+	
+	// Check the value of nbLinesStreaming
+	int nbLinesStreamingMax = 10000;
+	if (nbLinesStreamingMax > nbLinesSAR)
+	  {
+	    nbLinesStreamingMax = nbLinesSAR;
+	  }
+	if (nbLinesStreaming <= 0 || nbLinesStreaming >  nbLinesStreamingMax)
+	  {
+	    nbLinesStreaming =  nbLinesStreamingMax;
+	  } 
+	
+	// Construct the extendedPart
+	std::ostringstream os;
+	os << "?&streaming:type=stripped&streaming:sizevalue=" << nbLinesStreaming;
+	
+	// Add the extendedPart
+	std::string extendedPart = os.str();
+	extendedFileName.append(extendedPart);
+	
+	// Set the new FileName with extended options
+	SetParameterString("out", extendedFileName);
+      }
+
+    // Execute the main Pipeline
+    filterInterpolation->SetInput(GetParameterImage("indemproj"));
+    SetParameterOutputImage("out", filterInterpolation->GetOutput());
+  }
+  // Vector for filters 
+  std::vector<itk::ProcessObject::Pointer> m_Ref;
+
+};
+
+
+}
+
+}
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::SARInterpolation)

include/otbSARAmplitudeEstimationFunctor.h

@@ -66,6 +66,15 @@ public:
 			int nbCol_into_outputRegion, TOutputPixel * outValueAmplitude,
 			   bool isFirstMaille, double & tgElevationMaxForCurrentLine) 
   {
+    // Initialize
+    if (isFirstMaille)
+      {
+	for (unsigned int j = 0; j < nbCol_into_outputRegion; j++) 
+	  {
+	    outValueAmplitude[j] = static_cast<TOutputPixel>(0);
+	  }
+      }
+
     double aE, bE, cE, aS, bS, cS;
     double dy, dz, dc;
     int col1, col2;
@@ -159,7 +168,7 @@ public:
       }
        
     // Loop on colunm 
-    for (int k = std::max(col1, firstCol_into_outputRegion); k<= std::min(col2, nbCol_into_outputRegion+firstCol_into_outputRegion); k++)
+    for (int k = std::max(col1, firstCol_into_outputRegion); k<= std::min(col2, nbCol_into_outputRegion+firstCol_into_outputRegion-1); k++)
       {
 	a = (cS - k) / dc;                                                                                    
 	b = 1.0 - a;       	  

include/otbSARDEMPolygonsAnalysisImageFilter.txx

@@ -197,7 +197,7 @@ namespace otb
 	    itkExceptionMacro(<<"Output must be a Vector Image.");
 	    return;
 	  }
-
+	
 	// Set the number of Components
 	outputPtr->SetNumberOfComponentsPerPixel(numberOfEstimatedComponents);
       }
@@ -737,7 +737,7 @@ namespace otb
   ::ThreadedGenerateData(const ImageOutRegionType & outputRegionForThread,
 			 itk::ThreadIdType threadId)
   {
-    // Compute corresponding input region
+     // Compute corresponding input region
     ImageInRegionType inputRegionForThread = OutputRegionToInputRegion(outputRegionForThread);
     
     // Iterator on output (SAR geometry)
@@ -757,10 +757,7 @@ namespace otb
     ImageOutPixelType outValueTab[nbCol];
 
     // Initialize outValueTab
-    if (this->GetOutput()->GetNameOfClass() == "Image")
-      {
-	memset(outValueTab, 0., nbCol*sizeof(ImageOutPixelType));
-      }
+    m_FunctionOnPolygon->initalize(nbCol, outValueTab);
 
     // Allocate Vectors for polygon storage (vector of vectors to separate line and colunm of polygon into DEM)
     std::vector<std::vector<ImageInPixelType> * > Polygon_SideInUp_Vec;
@@ -930,8 +927,9 @@ namespace otb
 	Polygon_SideOutUp_Vec.clear();
 	Polygon_SideOutDown_Vec.clear();
 	
-	// Assignate the output with the contribution of the current polygon
+	// Assignate the output with the contribution of the current polygons
 	int counter = 0;
+	
 	while (!OutIt.IsAtEndOfLine()) 
 	  {
 	    // Affect OutIt
@@ -939,19 +937,13 @@ namespace otb
 	    ++counter;
 	    ++OutIt;
 	  } // End colunms (ouput)
-      
-	// reinitialize outValueTab
-	for (unsigned int j = 0; j < nbCol; j++) 
-	  {
-	    outValueTab[j] = static_cast<ImageOutPixelType>(0);
-	  }
-     
+
+
 	// Next Line
 	OutIt.NextLine();
      
       } // End lines (output)
- 
-  }
+   }
 
 
 

include/otbSARInterpolationFunctor.h

+/*
+ * Copyright (C) 2005-2018 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.
+ */
+
+#ifndef otbSARInterpolationFunctor_h
+#define otbSARInterpolationFunctor_h
+#include <cassert>
+
+#include "itkVariableLengthVector.h"
+#include "itkRGBPixel.h"
+#include "itkRGBAPixel.h"
+#include "itkObject.h"
+#include "itkObjectFactory.h"
+
+#include "otbSARPolygonsFunctor.h"
+
+namespace otb
+{
+namespace Function
+{
+/** \class SARInterpolationFunctor
+* \brief Base class for pixel representation functions.
+*
+*  \ingroup Visualization
+*
+* \sa AmplitudeFunctor PhaseFunctor
+ *
+ * \ingroup OTBCommon
+*/
+  template <class TInputPixel, class TOutputPixel>
+  class SARInterpolationFunctor : public SARPolygonsFunctor<TInputPixel, TOutputPixel>
+{
+public:
+  /** Standard class typedefs */
+  typedef SARInterpolationFunctor        Self;
+  typedef itk::Object                   Superclass;
+  typedef itk::SmartPointer<Self>       Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+  /** Method for creation through the object factory */
+  //itkNewMacro(Self);
+
+  /** Runtime information */
+  itkTypeMacro(SARInterpolationFunctor, itk::Object);
+
+  virtual void initalize(int nbElt, TOutputPixel * outValue)
+  {
+    int nbComponentsMax = this->GetNumberOfEstimatedComponents();
+
+    // Reserve Components into outValue and Initalize
+    for(int ind = 0; ind < nbElt; ind++)
+      {
+	outValue[ind].Reserve(nbComponentsMax);
+	
+	for (int iC = 0; iC < nbComponentsMax; iC++)
+	  {
+	    outValue[ind][iC] = 0;
+	  }	    
+      }
+  }
+
+  virtual void estimate(const int ind_LineSAR, TInputPixel * CLZY_InSideUpPtr, 
+			TInputPixel * CLZY_InSideDownPtr, 
+			TInputPixel * CLZY_OutSideUpPtr, 
+			TInputPixel * CLZY_OutSideDownPtr, 
+			int firstCol_into_outputRegion, 
+			int nbCol_into_outputRegion, TOutputPixel * outValue,
+			bool isFirstMaille, double & tgElevationMaxForCurrentLine) 
+  {
+    int nbComponentsMax = this->GetNumberOfEstimatedComponents();
+    // Initalize
+    if (isFirstMaille)
+      {
+	for(int ind = 0; ind < nbCol_into_outputRegion; ind++)
+	  {
+	    for (int iC = 0; iC < nbComponentsMax; iC++)
+	      {
+		outValue[ind][iC] = 0;
+	      }
+	    
+	  }
+      }
+
+
+    double aE, bE, cE, aS, bS, cS;
+    double dy, dz, dc;
+    int col1, col2;
+    // Elevation angle (for shadow)
+    double tgElevationMax = 0.;
+    double tgElevation = 0.;
+
+    double zE, yE, zS, yS;
+    double a, b, Z, Y;
+    double cotgIncidence;
+    double aux;
+
+    double a1, b1, c1, a2, b2, c2;
+    
+    TInputPixel CLZY_InSideUp = *CLZY_InSideUpPtr; 
+    TInputPixel CLZY_InSideDown = *CLZY_InSideDownPtr; 
+    TInputPixel CLZY_OutSideUp =  *CLZY_OutSideUpPtr;
+    TInputPixel CLZY_OutSideDown = *CLZY_OutSideDownPtr;
+
+    // Define the coef a and b
+    a1 =  CLZY_InSideUp[1] - ind_LineSAR;
+    a1 /= CLZY_InSideUp[1] - CLZY_InSideDown[1];
+
+    b1 = 1-a1;
+
+    c1 = b1*(CLZY_InSideUp[0]) + a1*(CLZY_InSideDown[0]);
+
+    a2 = CLZY_OutSideUp[1] - ind_LineSAR;
+    a2 /= (CLZY_OutSideUp[1] - CLZY_OutSideDown[1]);
+
+    b2 = 1-a2;
+
+    c2 = b2*(CLZY_OutSideUp[0]) + a2*(CLZY_OutSideDown[0]);
+     
+    // Caculate zE, yE, zS and yS
+    // Select input and output side for the current maille with c1 and c2 comparaison
+    if (c2 >= c1)
+      {
+	aE = a1;
+	bE = b1;
+	
+	zE = b1*(CLZY_InSideUp[2]) + 
+	  a1*(CLZY_InSideDown[2]);
+	yE = b1*(CLZY_InSideUp[3]) + 
+	  a1*(CLZY_InSideDown[3]);
+     
+	zS = b2*(CLZY_OutSideUp[2]) + 
+	  a2*(CLZY_OutSideDown[2]);
+	yS = b2*(CLZY_OutSideUp[3]) + 
+	  a2*(CLZY_OutSideDown[3]);
+
+	cE = c1;
+	cS = c2;
+      }
+    else
+      {
+	aS = a2;
+	bS = b2;
+
+	zE = b2*(CLZY_OutSideUp[2]) + 
+	  a2*(CLZY_OutSideDown[2]);
+	yE = b2*(CLZY_OutSideUp[3]) + 
+	  a2*(CLZY_OutSideDown[3]);
+
+	zS = b1*(CLZY_InSideUp[2]) + 
+	  a1*(CLZY_InSideDown[2]);
+	yS = b1*(CLZY_InSideUp[3]) + 
+	  a1*(CLZY_InSideDown[3]);
+	 
+	cE = c2;
+	cS = c1;
+      }
+      
+    dc = cS - cE;
+    dy = yS - yE;
+    dz = zS - zE;
+	      
+    // Colunm into // TODO: he current maille
+    col1 = (int) (std::min (cE, cS))+1;        // +1 ?                                                                                
+    col2 = (int) (std::max (cE, cS));
+     
+    // tgElevationMax per mailles (not sure)
+    if (isFirstMaille)
+      {
+	tgElevationMax = (zE > 0.0) ?  yE / zE : sign (1e20, yE);  
+	if (tgElevationMax < 0.) tgElevationMax = 0.;
+      }
+    else
+      {
+	tgElevationMax = tgElevationMaxForCurrentLine;
+      }
+       
+    // Loop on colunm 
+    for (int k = std::max(col1, firstCol_into_outputRegion); k<= std::min(col2, nbCol_into_outputRegion+firstCol_into_outputRegion-1); k++)
+      {
+	a = (cS - k) / dc;                                                                                    
+	b = 1.0 - a;       	  
+	Z = a * zE + b * zS;
+	Y = a * yE + b * yS;
+
+	 
+	tgElevation = (Z > 0.0) ? Y / Z : sign (1e20, Y);     
+	 
+	// If no into shadows
+	if (tgElevation > tgElevationMax)
+	  {                                                                         
+	    tgElevationMax = tgElevation;
+	     	     	     
+	    // Contribution for this maille
+	    int ind_for_out = k - firstCol_into_outputRegion;
+	    
+	    // Get the counter 
+	   int count = outValue[ind_for_out][0];
+	 
+	    // Nb Maille Max for each SAR pixel is m_MaxMaillePerPixel
+	    if (count < m_MaxMaillePerPixel)
+	      {
+		// Transfer the local Z and Y
+		outValue[ind_for_out][count*2 + 1] =  Z;
+		outValue[ind_for_out][count*2 + 2] =  Y;
+		
+		// Increment counter
+		outValue[ind_for_out][0] += 1;
+	      }
+	  }
+      }
+
+    // Store the tgElevationMax (for the next mailles)
+    tgElevationMaxForCurrentLine = tgElevationMax; 
+  }
+
+  /** Constructor */ // 
+  SARInterpolationFunctor(int maxMaillePerPixel)
+  {    
+    // Specific argument
+    m_MaxMaillePerPixel = maxMaillePerPixel;
+    
+    // Global argument (from PolygonsFunctor)
+    this->SetNumberOfExpectedComponents(4);
+    int estimatedComponents = 1 + m_MaxMaillePerPixel*2;
+    this->SetNumberOfEstimatedComponents(estimatedComponents);
+    this->SetSavePolygons(false);
+    this->SetOutputGeometry("SAR");
+  }
+
+  // Redefinition to use construction with maxMaillePerPixel
+static Pointer New(int maxMaillePerPixel) 
+  {
+  Pointer smartPtr = ::itk::ObjectFactory< Self >::Create();    
+    if ( smartPtr == nullptr )				     
+       {                                                
+         smartPtr = new Self(maxMaillePerPixel);                                
+        }                                                
+       smartPtr->UnRegister();                            
+       return smartPtr;
+  }
+
+  /** Destructor */
+  ~SARInterpolationFunctor() override {}
+  
+  
+private :
+  int m_MaxMaillePerPixel;
+
+  };
+
+}
+}
+
+#endif

include/otbSARPolygonsFunctor.h

@@ -52,6 +52,10 @@ public:
   /** Runtime information */
   itkTypeMacro(SARPolygonsFunctor, itk::Object);
 
+  void initalize(int nbElt, TOutputPixel * outValueAmplitude)
+  {
+  }
+  
   // Estimate method
   virtual void estimate(const int ind_LineSAR, TInputPixel * CLZY_InSideUpPtr, 
 			TInputPixel * CLZY_InSideDownPtr,