ENH : New filter to grouped by interferogram estimation into ML geometry...

ENH : New filter to grouped by interferogram estimation into ML geometry (Compensated Complex image as input)

include/otbSARGroupedByMLImageFilter.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 otbSARGroupedByMLImageFilter_h
+#define otbSARGroupedByMLImageFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkSmartPointer.h"
+#include "itkPoint.h"
+
+#include "itkSimpleFastMutexLock.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageScanlineIterator.h"
+
+#include "otbImageKeywordlist.h"
+
+
+namespace otb
+{
+/** \class SARGroupedByMLImageFilter 
+ * \brief Creates an interferogram into ml geometry between two SAR images. 
+ * 
+ * This filter built the ML interferogram between a master and a slave image.
+ *
+ * The output is a vector image with amplitude, phase coherency and occurrences. 
+ * \ingroup DiapOTBModule
+ */
+
+  template <typename TImageSAR, typename TImageOut> 
+  class ITK_EXPORT SARGroupedByMLImageFilter :
+    public itk::ImageToImageFilter<TImageSAR,TImageOut>
+{
+public:
+
+  // Standard class typedefs
+  typedef SARGroupedByMLImageFilter                    Self;
+  typedef itk::ImageToImageFilter<TImageSAR,TImageOut>    Superclass;
+  typedef itk::SmartPointer<Self>                        Pointer;
+  typedef itk::SmartPointer<const Self>                  ConstPointer;
+
+  // Method for creation through object factory
+  itkNewMacro(Self);
+  // Run-time type information
+  itkTypeMacro(SARGroupedByMLImageFilter,ImageToImageFilter);
+
+  /** Typedef to image input type (master and slave reech images) :  otb::Image (Complex)  */
+  typedef TImageSAR                                  ImageSARType;
+  /** Typedef to describe the inout image pointer type. */
+  typedef typename ImageSARType::Pointer             ImageSARPointer;
+  typedef typename ImageSARType::ConstPointer        ImageSARConstPointer;
+  /** Typedef to describe the inout image region type. */
+  typedef typename ImageSARType::RegionType          ImageSARRegionType;
+  /** Typedef to describe the type of pixel and point for inout image. */
+  typedef typename ImageSARType::PixelType           ImageSARPixelType;
+  typedef typename ImageSARType::PointType           ImageSARPointType;
+  /** Typedef to describe the image index, size types and spacing for inout image. */
+  typedef typename ImageSARType::IndexType           ImageSARIndexType;
+  typedef typename ImageSARType::IndexValueType      ImageSARIndexValueType;
+  typedef typename ImageSARType::SizeType            ImageSARSizeType;
+  typedef typename ImageSARType::SizeValueType       ImageSARSizeValueType;
+  typedef typename ImageSARType::SpacingType         ImageSARSpacingType;
+  typedef typename ImageSARType::SpacingValueType    ImageSARSpacingValueType;
+  
+  /** Typedef to image output type : otb::VectorImage  */
+  typedef TImageOut                                  ImageOutType;
+  /** Typedef to describe the output image pointer type. */
+  typedef typename ImageOutType::Pointer             ImageOutPointer;
+  typedef typename ImageOutType::ConstPointer        ImageOutConstPointer;
+  /** Typedef to describe the output image region type. */
+  typedef typename ImageOutType::RegionType          ImageOutRegionType;
+  /** Typedef to describe the type of pixel and point for output image. */
+  typedef typename ImageOutType::PixelType           ImageOutPixelType;
+  typedef typename ImageOutType::PointType           ImageOutPointType;
+  /** Typedef to describe the image index, size types and spacing for output image. */
+  typedef typename ImageOutType::IndexType           ImageOutIndexType;
+  typedef typename ImageOutType::IndexValueType      ImageOutIndexValueType;
+  typedef typename ImageOutType::SizeType            ImageOutSizeType;
+  typedef typename ImageOutType::SizeValueType       ImageOutSizeValueType;
+  typedef typename ImageOutType::SpacingType         ImageOutSpacingType;
+  typedef typename ImageOutType::SpacingValueType    ImageOutSpacingValueType;
+
+  // Define Point2DType and Point3DType
+  using Point2DType = itk::Point<double,2>;
+  using Point3DType = itk::Point<double,3>;
+
+  // Typedef for iterators
+  typedef itk::ImageScanlineConstIterator< ImageSARType > InputSARIterator;
+  typedef itk::ImageScanlineIterator< ImageOutType > OutputIterator;
+
+  typedef itk::SimpleFastMutexLock         MutexType;
+
+  // Setter
+  itkSetMacro(MLRan, unsigned int);
+  itkSetMacro(MLAzi, unsigned int);
+  itkSetMacro(MarginRan, unsigned int);
+  itkSetMacro(MarginAzi, unsigned int);
+  itkSetMacro(Gain, float);
+
+  // Getter
+  itkGetMacro(MLRan, unsigned int);
+  itkGetMacro(MLAzi, unsigned int);
+  itkGetMacro(MarginRan, unsigned int);
+  itkGetMacro(MarginAzi, unsigned int);
+  itkGetMacro(Gain, float);
+
+
+protected:
+  // Constructor
+  SARGroupedByMLImageFilter();
+
+  // Destructor
+  virtual ~SARGroupedByMLImageFilter() ITK_OVERRIDE;
+
+  // Print
+  void PrintSelf(std::ostream & os, itk::Indent indent) const ITK_OVERRIDE;
+  
+  /** SARGroupedByMLImageFilter produces an image which are into output geometry. The differences between 
+   * output and input images can be the size of images and the dimensions. 
+   * As such, SARGroupedByMLImageFilter needs to provide an implementation for 
+   * GenerateOutputInformation() in order to inform the pipeline execution model. 
+   */ 
+  virtual void GenerateOutputInformation() ITK_OVERRIDE;
+
+  /** SARGroupedByMLImageFilter needs a input requested region that corresponds to the margin and ML factors 
+   * into the requested region of our output requested region. 
+   * As such, SARGroupedByMLImageFilter needs to provide an implementation for 
+   * GenerateInputRequestedRegion() in order to inform the pipeline execution model. 
+   * \sa ProcessObject::GenerateInputRequestedRegion() */
+  virtual void GenerateInputRequestedRegion() ITK_OVERRIDE;
+
+  /**
+   * OutputRegionToInputRegion returns the input SAR region 
+   */
+  ImageSARRegionType OutputRegionToInputRegion(const ImageOutRegionType& outputRegion, bool withMargin) const;
+
+
+  /** 
+   * SARGroupedByMLImageFilterr can be implemented as a multithreaded filter.
+   * Therefore, this implementation provides a ThreadedGenerateData() routine
+   * which is called for each processing thread. The main output image data is
+   * allocated automatically by the superclass prior to calling
+   * ThreadedGenerateData().  ThreadedGenerateData can only write to the
+   * portion of the output image specified by the parameter
+   * "outputRegionForThread"
+   *
+   * \sa ImageToImageFilter::ThreadedGenerateData(),
+   *     ImageToImageFilter::GenerateData() */
+  virtual void ThreadedGenerateData(const ImageOutRegionType& outputRegionForThread, 
+				    itk::ThreadIdType threadId) ITK_OVERRIDE;
+  
+
+
+ private:
+  SARGroupedByMLImageFilter(const Self&); // purposely not implemented
+  void operator=(const Self &); // purposely not 
+  
+  // ML factor in range 
+  unsigned int m_MLRan;
+  
+  // ML factor in azimut
+  unsigned int m_MLAzi;
+  
+  // gain for amplitude estimation
+  float m_Gain;
+
+  // SAR Dimensions
+  int m_nbLinesSAR;
+  int m_nbColSAR;
+
+  // Margin in order to spend our averaging
+  unsigned int m_MarginRan;
+  unsigned int m_MarginAzi;
+
+};
+
+} // End namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSARGroupedByMLImageFilter.txx"
+#endif
+
+
+
+#endif

include/otbSARGroupedByMLImageFilter.txx

+/*
+ * 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 otbSARGroupedByMLImageFilter_txx
+#define otbSARGroupedByMLImageFilter_txx
+
+#include "otbSARGroupedByMLImageFilter.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageScanlineIterator.h"
+#include "itkProgressReporter.h"
+#include "itkNumericTraitsPointPixel.h"
+#include "itkContinuousIndex.h"
+
+#include <cmath>
+#include <algorithm>
+
+
+namespace otb
+{
+  /** 
+   * Constructor with default initialization
+   */
+  template <class TImageSAR, class TImageOut> 
+  SARGroupedByMLImageFilter< TImageSAR, TImageOut >::SARGroupedByMLImageFilter()
+    :  m_MLRan(1), m_MLAzi(1), m_Gain(1), m_nbLinesSAR(0), m_nbColSAR(0), 
+       m_MarginRan(1), m_MarginAzi(1)
+  {
+  }
+    
+  /** 
+   * Destructor
+   */
+  template <class TImageSAR, class TImageOut> 
+  SARGroupedByMLImageFilter< TImageSAR, TImageOut >::~SARGroupedByMLImageFilter()
+  {
+  }
+
+  /**
+   * Print
+   */
+  template<class TImageSAR, class TImageOut>
+  void
+  SARGroupedByMLImageFilter< TImageSAR, TImageOut >
+  ::PrintSelf(std::ostream & os, itk::Indent indent) const
+  {
+    Superclass::PrintSelf(os, indent);
+
+    os << indent << "ML factors : " << m_MLRan << ", " << m_MLAzi << std::endl;
+    os << indent << "Gain for amplitude estimation : " << m_Gain << std::endl;
+  }
+
+
+  /**
+   * Method GenerateOutputInformaton()
+   **/
+  template<class TImageSAR, class TImageOut>
+  void
+  SARGroupedByMLImageFilter< TImageSAR, TImageOut >
+  ::GenerateOutputInformation()
+  {
+    // Call superclass implementation
+    Superclass::GenerateOutputInformation();
+
+    // Get pointers to the inputs and output
+    ImageSARConstPointer inputPtr = this->GetInput();
+    ImageOutPointer outputPtr = this->GetOutput();
+
+    // KeyWordList
+    ImageKeywordlist masterKWL = inputPtr->GetImageKeywordlist();
+
+    // Input Dimensions
+    m_nbLinesSAR = this->GetInput()->GetLargestPossibleRegion().GetSize()[1];
+    m_nbColSAR = this->GetInput()->GetLargestPossibleRegion().GetSize()[0];
+  
+    //////////////////////////// Main Output : GroupedByML into ML Geo ///////////////////////////
+    // Vector Image  :
+    // At Least 4 Components :  
+    //                _ Amplitude
+    //                _ Phase
+    //                _ Coherence
+    //                _ IsData
+    outputPtr->SetNumberOfComponentsPerPixel(4);
+
+    // The output is defined with the Input Image (Compensated Complex image into SAR geomtry) 
+    // and by ML factors 
+    // Origin, Spacing and Size (SAR master geometry) in function of ML factors
+    ImageOutSizeType outputSize;
+
+    outputSize[0] = std::max<ImageOutSizeValueType>(static_cast<ImageOutSizeValueType>(std::floor( (double) m_nbColSAR / (double) m_MLRan)),1);
+    outputSize[1] = std::max<ImageOutSizeValueType>(static_cast<ImageOutSizeValueType>(std::floor( (double) m_nbLinesSAR / (double) m_MLAzi)),1);
+
+    ImageOutPointType outOrigin;
+    outOrigin = inputPtr->GetOrigin();
+    ImageOutSpacingType outSP;
+    outSP = inputPtr->GetSpacing();
+    outSP[0] *= m_MLRan;
+    outSP[1] *= m_MLAzi;
+
+    // Define Output Largest Region
+    ImageOutRegionType outputLargestPossibleRegion = inputPtr->GetLargestPossibleRegion();
+    outputLargestPossibleRegion.SetSize(outputSize);
+    outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
+    outputPtr->SetOrigin(outOrigin);
+    outputPtr->SetSpacing(outSP);
+
+    std::cout << "outputLargestPossibleRegion : " << outputLargestPossibleRegion << std::endl;
+
+    // Add ML factors and bands meaning into keyWordList
+    ImageKeywordlist outputKWL = masterKWL;
+    outputKWL.AddKey("support_data.band.Amplitude", std::to_string(0));
+    outputKWL.AddKey("support_data.band.Phase", std::to_string(1));
+    outputKWL.AddKey("support_data.band.Coherency", std::to_string(2));
+    outputKWL.AddKey("support_data.band.isData", std::to_string(3));
+
+    outputKWL.AddKey("support_data.ml_ran", std::to_string(m_MLRan));
+    outputKWL.AddKey("support_data.ml_azi", std::to_string(m_MLAzi));
+
+    // Set new keyword list to output image
+    outputPtr->SetImageKeywordList(outputKWL);    
+    
+
+     ///////// Checks margin with ML factors (for averaging) /////////////
+     if (m_MLRan != 0 || m_MLAzi != 0)
+       {
+	 //if (m_MLRan < m_MarginRan || m_MLAzi < m_MarginAzi)
+	 // {
+	 //   itkExceptionMacro(<<"Margin for averaging can't be superior to ML Factors.");   
+	 // } 
+       }
+  }
+
+  /** 
+   * Method OutputRegionToInputRegion for GenerateInputRequestedRegion
+   */
+  template<class TImageSAR, class TImageOut>
+  typename SARGroupedByMLImageFilter< TImageSAR, TImageOut >::ImageSARRegionType 
+  SARGroupedByMLImageFilter< TImageSAR, TImageOut >
+  ::OutputRegionToInputRegion(const ImageOutRegionType& outputRegion, bool withMargin) const
+  {
+    // Compute the input requested region (size and start index)
+    // Use the image transformations to insure an input requested region
+    // that will provide the proper range
+    const ImageOutSizeType & outputRequestedRegionSize = outputRegion.GetSize();
+    ImageOutIndexType outputRequestedRegionIndex = outputRegion.GetIndex();
+
+    // Add a margin for averaging
+    int marge_MLRan = m_MarginRan;
+    int marge_MLAzi = m_MarginAzi;
+
+    if (!withMargin)
+      {
+	marge_MLRan = 0;
+	marge_MLAzi = 0;
+      }
+
+    // Multiply to output region by ML factor to get input SAR region
+    ImageSARIndexType inputRequestedRegionIndex; 
+    inputRequestedRegionIndex[0] = static_cast<ImageSARIndexValueType>(outputRequestedRegionIndex[0] * m_MLRan
+								       - marge_MLRan);
+    inputRequestedRegionIndex[1] = static_cast<ImageSARIndexValueType>(outputRequestedRegionIndex[1] * m_MLAzi 
+								       - marge_MLAzi);
+    ImageSARSizeType inputRequestedRegionSize;
+    inputRequestedRegionSize[0] = static_cast<ImageSARSizeValueType>(outputRequestedRegionSize[0] * m_MLRan + 
+								     2*marge_MLRan);
+    inputRequestedRegionSize[1] = static_cast<ImageSARSizeValueType>(outputRequestedRegionSize[1] * m_MLAzi + 
+								     2*marge_MLAzi);  
+   
+    // Check Index and Size
+    if (inputRequestedRegionIndex[0] < this->GetInput()->GetLargestPossibleRegion().GetIndex()[0])
+      {
+	inputRequestedRegionIndex[0] = this->GetInput()->GetLargestPossibleRegion().GetIndex()[0];
+      }
+    if (inputRequestedRegionIndex[1] < this->GetInput()->GetLargestPossibleRegion().GetIndex()[1])
+      {
+	inputRequestedRegionIndex[1] = this->GetInput()->GetLargestPossibleRegion().GetIndex()[1];
+      }
+    if ((inputRequestedRegionSize[0] + inputRequestedRegionIndex[0]) > 
+	this->GetInput()->GetLargestPossibleRegion().GetSize()[0])
+      {
+	inputRequestedRegionSize[0] = this->GetInput()->GetLargestPossibleRegion().GetSize()[0] - 
+	  inputRequestedRegionIndex[0];
+      }
+    if ((inputRequestedRegionSize[1] + inputRequestedRegionIndex[1]) > 
+	this->GetInput()->GetLargestPossibleRegion().GetSize()[1])
+      {
+	inputRequestedRegionSize[1] = this->GetInput()->GetLargestPossibleRegion().GetSize()[1] - 
+	  inputRequestedRegionIndex[1];
+      }
+
+    // Transform into a region
+    ImageSARRegionType inputRequestedRegion = outputRegion;
+    inputRequestedRegion.SetIndex(inputRequestedRegionIndex);
+    inputRequestedRegion.SetSize(inputRequestedRegionSize);
+    
+    
+    return inputRequestedRegion;  
+  }
+
+
+  /** 
+   * Method GenerateInputRequestedRegion
+   */
+  template<class TImageSAR, class TImageOut>
+  void
+  SARGroupedByMLImageFilter< TImageSAR, TImageOut >
+  ::GenerateInputRequestedRegion()
+  {
+    // call the superclass' implementation of this method
+    Superclass::GenerateInputRequestedRegion();
+    
+    // Get Output requested region
+    ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
+
+    ///////////// Find the region into SAR input images (same geometry for SAR after coregistration ////////////
+    ImageSARRegionType inputSARRequestedRegion = OutputRegionToInputRegion(outputRequestedRegion, true);
+
+    ImageSARPointer  inputPtr = const_cast< ImageSARType * >( this->GetInput() );
+    
+    inputPtr->SetRequestedRegion(inputSARRequestedRegion);
+  }
+  
+
+  /**
+   * Method ThreadedGenerateData
+   */
+  template<class TImageSAR, class TImageOut>
+  void
+  SARGroupedByMLImageFilter< TImageSAR, TImageOut >
+  ::ThreadedGenerateData(const ImageOutRegionType & outputRegionForThread,
+			 itk::ThreadIdType threadId)
+  {
+    // Compute corresponding input region for SAR image (Copensated Complex Image) 
+    ImageSARRegionType inputRegionForThread = OutputRegionToInputRegion(outputRegionForThread, true);
+  
+    //  Define/declare an iterator that will walk the input regions for this
+    // thread. 
+    InputSARIterator  inIt(this->GetInput(), inputRegionForThread);
+  
+    // Define/declare an iterator that will walk the output region for this
+    // thread.
+    OutputIterator outIt(this->GetOutput(), outputRegionForThread);
+
+    // Support progress methods/callbacks
+    itk::ProgressReporter progress( this, threadId, outputRegionForThread.GetNumberOfPixels() );
+
+    inIt.GoToBegin();
+    outIt.GoToBegin();
+
+    // Get the number of column of the region
+    const ImageOutSizeType &   outputRegionForThreadSize = outputRegionForThread.GetSize();
+    unsigned int nbCol = static_cast<unsigned int>(outputRegionForThreadSize[0]); 
+  
+    // Allocate Arrays for accumulations (size = nbCol)
+    double *complexMulConjTab_Re = new double[nbCol];
+    double *complexMulConjTab_Im = new double[nbCol];
+    double *complexMulConjTab_Mod = new double[nbCol];
+    
+    // Output Pixel (VectorImage Pixel)
+    ImageOutPixelType outPixel;
+    outPixel.Reserve(4);
+
+    //// Initialize the indLastL and indFirstL (for spending our averaging with a margin) ////
+    int indFirstL = 0;
+    int indLastL = m_MLAzi+ m_MarginAzi;
+    int marginAzi = (int) m_MarginAzi;
+    
+    if (inIt.GetIndex()[1] >= marginAzi && marginAzi != 0)
+      {
+	indFirstL = -marginAzi;
+      }
+	
+    if ((inIt.GetIndex()[1] + indLastL - indFirstL) > (m_nbLinesSAR-1))
+      {
+	indLastL = (m_nbLinesSAR-1) - inIt.GetIndex()[1] + indFirstL;
+      }
+
+    bool backTopreviousLines = false;
+    bool backTopreviousColunms = false;
+    
+
+    // For each Line
+    while ( !inIt.IsAtEnd() && !outIt.IsAtEnd())
+      {
+	// reinitialize
+	for (unsigned int j = 0; j < nbCol; j++) 
+	  {
+	    complexMulConjTab_Re[j] = 0;
+	    complexMulConjTab_Im[j] = 0;
+	    complexMulConjTab_Mod[j] = 0;
+	  }
+
+	if (backTopreviousLines)
+	  {
+	    //// Previous index in azimut (for averaging) ////
+	    if (inIt.GetIndex()[1] >= marginAzi && marginAzi != 0)
+	      {
+		indFirstL = -marginAzi;
+	      }
+
+	    ImageSARIndexType indSARL;
+	    indSARL[0] = inIt.GetIndex()[0];
+	    indSARL[1] = inIt.GetIndex()[1] - (2*marginAzi);
+	
+	    if ((indSARL[1] + indLastL - indFirstL) > (m_nbLinesSAR-1))
+	      {
+		indLastL = (m_nbLinesSAR - 1) - indSARL[1] + indFirstL;
+	      }
+
+	    // Check previous index
+	    if (indSARL[1] >= 0 && marginAzi != 0)
+	      {
+		// Previous index inputs
+		inIt.SetIndex(indSARL);
+		
+		indFirstL = -marginAzi;
+	      }
+	  }
+
+	backTopreviousLines = true;
+
+
+	// GenerateInputRequestedRegion function requires an input region with In_nbLine = m_Averaging * Out_nbLine
+	for (int i = indFirstL ; i < indLastL; i++) 
+	  {
+	    inIt.GoToBeginOfLine();
+	    outIt.GoToBeginOfLine();
+	  
+	    int colCounter = 0;
+	    backTopreviousColunms = false;
+
+	    //// Initialize the indLastC and indFirstC (for spending average with a margin) ////
+	    int indFirstC = 0;
+	    int indLastC = m_MLRan + m_MarginRan;
+	    int marginRan = (int) m_MarginRan;
+    
+	    if (inIt.GetIndex()[0] >= marginRan  && marginRan != 0)
+	      {
+		indFirstC = -marginRan;
+	      }
+	
+	    if ((inIt.GetIndex()[0] + indLastC - indFirstC) > (m_nbColSAR-1))
+	      {
+		indLastC = (m_nbColSAR-1) - inIt.GetIndex()[0] + indFirstC;
+	      }
+
+
+	    // For each column
+	    while (!inIt.IsAtEndOfLine() && !outIt.IsAtEndOfLine())
+	      {
+		if (backTopreviousColunms)
+		  {
+		    //// Previous index in range (for averaging) ////
+		    if (inIt.GetIndex()[0] >= marginRan  && marginRan != 0)
+		      {
+			indFirstC = -marginRan;
+		      }
+
+		    ImageSARIndexType indSARC;
+		    indSARC[0] = inIt.GetIndex()[0] - (2*marginRan);
+		    indSARC[1] = inIt.GetIndex()[1];
+		
+		    if ((indSARC[0] + indLastC - indFirstC) > (m_nbColSAR-1))
+		      {
+			indLastC = (m_nbColSAR-1) - indSARC[0] + indFirstC;
+		      }
+
+		    // Check previous index
+		    if (indSARC[0] >= 0 && marginRan != 0)
+		      {
+			// Previous index inputs
+			inIt.SetIndex(indSARC);
+		
+			indFirstC = -marginRan;
+		      }
+	
+		  }
+
+		backTopreviousColunms = true;
+
+		for (int k = indFirstC ; k < indLastC; k++) 
+		  {
+		    if (!inIt.IsAtEndOfLine())
+		      {
+			// Complex interferogram (master * conj(slave) : get from inout image)
+			double real_interfero = inIt.Get().real();
+
+			double imag_interfero = inIt.Get().imag();
+
+			///////////// Accumulations ///////////////
+			complexMulConjTab_Re[colCounter] += real_interfero;
+
+			
+			complexMulConjTab_Im[colCounter] += imag_interfero;
+			
+			
+			complexMulConjTab_Mod[colCounter] += sqrt(real_interfero*real_interfero +
+								  imag_interfero*imag_interfero);
+
+			// Next colunm inputs
+			++inIt;
+		      }
+		  }
+	      	      
+		// Estiamte and Assigne outIt
+		if (i == (indLastL-1))
+		  {
+		    ///////////// Estimations of amplitude, phase and coherency ///////////////
+		    double mod_Acc = sqrt(complexMulConjTab_Re[colCounter]*complexMulConjTab_Re[colCounter] + 
+					  complexMulConjTab_Im[colCounter]*complexMulConjTab_Im[colCounter]);
+
+		    int countPixel = (m_MLRan + 2*m_MarginRan) * (m_MLAzi + 2*m_MarginAzi);
+
+		    // Amplitude
+		    outPixel[0] = m_Gain * sqrt((complexMulConjTab_Mod[colCounter]/(countPixel)));
+		    
+		    // Phase
+		    outPixel[1] = std::atan2(complexMulConjTab_Im[colCounter], 
+					     complexMulConjTab_Re[colCounter]);
+
+		    // Mod 2*Pi
+		    outPixel[1] =  outPixel[1]-(2*M_PI)*floor(outPixel[1]/(2*M_PI));
+		    
+		    // Coherency
+		    outPixel[2] = mod_Acc / complexMulConjTab_Mod[colCounter] ;
+	  
+		    // IsData always set to 1
+		    outPixel[3] = 1;
+		    
+		    if (outPixel[0] == 0 && outPixel[1] == 0 && outPixel[2] == 0)
+		      {
+			outPixel[3] = 0;
+		      } 
+		    
+		    // Assigne Main output (ML geometry) 
+		    outIt.Set(outPixel);
+		    progress.CompletedPixel();		  
+		  }
+
+		// Next colunm output
+		++outIt;
+		colCounter++;
+	      }
+	  
+	    // Next line intputs
+	    inIt.NextLine();
+	  }
+      
+	// Next line output
+	outIt.NextLine();
+   
+      }
+    delete [] complexMulConjTab_Re;
+    delete[ ] complexMulConjTab_Im;
+    delete[ ] complexMulConjTab_Mod; 
+  }
+  
+
+} /*namespace otb*/
+
+#endif