otbSARDEMPolygonsAnalysisImageFilter.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 otbSARDEMPolygonsAnalysisImageFilter_txx
#define otbSARDEMPolygonsAnalysisImageFilter_txx

#include "otbSARDEMPolygonsAnalysisImageFilter.h"

#include "itkImageScanlineConstIterator.h"
#include "itkImageScanlineIterator.h"
#include "itkProgressReporter.h"
#include "itkNumericTraitsPointPixel.h"

#include "otbVectorImage.h"

#include <cmath>
#include <type_traits>


namespace otb
{
  /** 
   * Constructor with default initialization
   */
  template <class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction> 
  SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >::SARDEMPolygonsAnalysisImageFilter()
    : m_SarSensorModelAdapter(ITK_NULLPTR), m_Gain(100.), m_Margin(0), m_NoData(-32768), m_MLRan(1), m_MLAzi(1)
  {
  }
  
  /**
   * Set Sar Image keyWordList
   */ 
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void
  SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >
  ::SetSARImageKeyWorList(ImageKeywordlist sarImageKWL)
  {
    m_SarImageKwl = sarImageKWL;
  }
  
  /**
   * Set Sar Image Ptr
   */ 
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void
  SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >
  ::SetSARImagePtr(ImageSARPointer sarPtr)
  {
    // Check if sarImageKWL not NULL
    assert(sarPtr && "SAR Image doesn't exist.");
    m_SarImagePtr = sarPtr;
  }

  /**
   * Set Dem Image Ptr
   */ 
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void
  SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >
  ::SetDEMImagePtr(ImageDEMPointer demPtr)
  {
    // Check if sarImageKWL not NULL
    assert(demPtr && "DEM Image doesn't exist.");
    m_DemImagePtr = demPtr;
  }

  /**
   * Set Dem Image Ptr
   */ 
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void
  SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >
  ::SetFunctorPtr(FunctionPointer functor)
  {
    m_FunctionOnPolygon = functor;
  }


  /**
   * Set DEM Information
   */ 
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void
  SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >
  ::SetDEMInformation(double gain, int DEMDirC, int DEMDirL)
  {
    // Check if gain is greater than 0
    assert(gain >= 0 && "gain must be greater than 0.");
    // Check if directions are equal to 1 or -1
    assert((DEMDirC == 1 || DEMDirC == -1)  && "direction to scan the DEM in range must be 1 or -1.");
    assert((DEMDirL == 1 || DEMDirL == -1)  && "direction to scan the DEM in aimut must be 1 or -1.");
    
    m_Gain = gain;
    m_DEMdirL= DEMDirL;
    m_DEMdirC = DEMDirC;
  }

  /**
   * Initialize margin and GenericRSTransform
   */ 
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void
  SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >
  ::initializeMarginAndRSTransform()
  {
    // Initialize the GenericRSTransform
    m_RSTransform = RSTransformType2D::New();
    m_RSTransform->SetInputKeywordList( m_SarImagePtr->GetImageKeywordlist() );
    m_RSTransform->SetInputProjectionRef( m_SarImagePtr->GetProjectionRef() );
    m_RSTransform->SetOutputProjectionRef(m_DemImagePtr->GetProjectionRef());
    m_RSTransform->InstantiateTransform();

    // Set the margin
    int nbLinesSAR = m_SarImagePtr->GetLargestPossibleRegion().GetSize()[1];
    int nbLinesDEM = m_DemImagePtr->GetLargestPossibleRegion().GetSize()[1];
    if (nbLinesSAR/2 > nbLinesDEM)
      { 
	m_Margin = (nbLinesSAR/nbLinesDEM)*30;
      }
    else if (nbLinesSAR > nbLinesDEM)
      {
	m_Margin = (nbLinesSAR/nbLinesDEM)*60;
      }
    else
      {
	m_Margin = 100;
      }
  }

  /**
   * Print
   */
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void
  SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction >
  ::PrintSelf(std::ostream & os, itk::Indent indent) const
  {
    Superclass::PrintSelf(os, indent);

    os << indent << "Gain : " << m_Gain << std::endl;
    os << indent << "DEM scan direction in range : " << m_DEMdirC << std::endl;
    os << indent << "DEM scan direction in aeimut : " << m_DEMdirL << std::endl;
  }

  /**
   * Method GenerateOutputInformaton()
   **/
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void
  SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction >
  ::GenerateOutputInformation()
  {
    // Call superclass implementation
    Superclass::GenerateOutputInformation();

    // Get pointers to the input and output
    ImageInConstPointer inputPtr = this->GetInput();
    ImageOutPointer     outputPtr = this->GetOutput();

    ImageKeywordlist inputKwl = inputPtr->GetImageKeywordlist();
    
    // SAR dimensions
    m_NbLinesSAR = m_SarImagePtr->GetLargestPossibleRegion().GetSize()[1];
    m_NbColSAR = m_SarImagePtr->GetLargestPossibleRegion().GetSize()[0];

    m_FunctionOnPolygon->SetSARDimensions(m_NbColSAR, m_NbLinesSAR);


    //////// Check Input/Output in function of Functor //////// 
    //////// The functor can be required a classic image or a vector Image as main output ////////
    if (strcmp(inputPtr->GetNameOfClass(), "VectorImage")  != 0)
      {
	itkExceptionMacro(<<"Input must be a Vector Image.");
	return;
      }

    unsigned int numberOfExpectedComponents = m_FunctionOnPolygon->GetNumberOfExpectedComponents();
    if (numberOfExpectedComponents > inputPtr->GetNumberOfComponentsPerPixel())
      {
	itkExceptionMacro(<<"Number of Components by pixels of the input is not enough for the used funtor.");
	return;
      }
    
    unsigned int numberOfEstimatedComponents = m_FunctionOnPolygon->GetNumberOfEstimatedComponents();
    if (numberOfEstimatedComponents > 1)
      {	
	if (strcmp(outputPtr->GetNameOfClass(), "VectorImage") != 0)
	  {
	    itkExceptionMacro(<<"Output must be a Vector Image.");
	    return;
	  }
	
	// Set the number of Components
	outputPtr->SetNumberOfComponentsPerPixel(numberOfEstimatedComponents);
      }
    else
      {
	if (strcmp(outputPtr->GetNameOfClass(), "Image") != 0)
	  {
	    itkExceptionMacro(<<"Output must be a Vector Image.");
	    return;
	  }
      }


    ////////// Check the nature of Components into projeted DEM //////////
    // The functor requires several components (with GetRequiredComponents). These components have to be into 
    // projeted DEM and correspond to a band into our input.
    std::string BandKey = "support_data.band.";
    std::vector<std::string> vecRequiredComponents = m_FunctionOnPolygon->GetRequiredComponents();
    std::map<std::string, int> mapRequiredComponents;
    
    for (unsigned int i = 0; i < vecRequiredComponents.size(); i++)
      {
	std::string component = BandKey + vecRequiredComponents.at(i);
	
	if (!inputKwl.HasKey(component))
	  {
	    itkExceptionMacro(<<"Missing Component into Projeted DEM.");
	    return;
	  }
	else
	  {
	    mapRequiredComponents[vecRequiredComponents.at(i)] = atoi(inputKwl.GetMetadataByKey(component).c_str());
	  }
      } 
    
    m_FunctionOnPolygon->SetRequiredComponentsToInd(mapRequiredComponents);
    

    ////////// Define geometry output with the used functor //////////
    // Two kinds of output geometries : SAR or ML
    m_FunctionOnPolygon->GetOutputGeometry(m_OutGeometry, m_MLRan, m_MLAzi);
    
    if (m_OutGeometry == "SAR")
      {	
	// The output is defined with the m_SarImageKwl 
	// Origin, Spacing and Size (SAR image geometry)
	ImageOutSizeType outputSize;
	outputSize[0] = m_NbColSAR;
	outputSize[1] = m_NbLinesSAR;
	ImageOutPointType outOrigin;
	outOrigin = m_SarImagePtr->GetOrigin();
	ImageOutSpacingType outSP;
	outSP = m_SarImagePtr->GetSignedSpacing();

	// Define Output Largest Region
	ImageOutRegionType outputLargestPossibleRegion = m_SarImagePtr->GetLargestPossibleRegion();
	outputLargestPossibleRegion.SetSize(outputSize);
	outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
	outputPtr->SetOrigin(outOrigin);
	outputPtr->SetSignedSpacing(outSP);
	outputPtr->SetDirection(m_SarImagePtr->GetDirection());

	// Set new keyword list to output image
	outputPtr->SetImageKeywordList(m_SarImageKwl);
      }
    else if (m_OutGeometry == "ML")
      {
	// The output is defined with the m_SarImageKwl and by ML factors 
	// Origin, Spacing and Size (SAR image 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 = m_SarImagePtr->GetOrigin();
	ImageOutSpacingType outSP;
	outSP = m_SarImagePtr->GetSignedSpacing();
	outSP[0] *= m_MLRan;
	outSP[1] *= m_MLAzi;

	// Define Output Largest Region
	ImageOutRegionType outputLargestPossibleRegion = m_SarImagePtr->GetLargestPossibleRegion();
	outputLargestPossibleRegion.SetSize(outputSize);
	outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
	outputPtr->SetOrigin(outOrigin);
	outputPtr->SetSignedSpacing(outSP);
	outputPtr->SetDirection(m_SarImagePtr->GetDirection());

	// Add ML factors into keyWordList
	m_SarImageKwl.AddKey("support_data.ml_ran", std::to_string(m_MLRan));
	m_SarImageKwl.AddKey("support_data.ml_azi", std::to_string(m_MLAzi));

	// Set new keyword list to output image
	outputPtr->SetImageKeywordList(m_SarImageKwl);
      }
    else
      {
	itkExceptionMacro(<<"Wrong Output Geometry.");
	return;
      }


    ////////// Add the estimated bands to keywordList //////////
    ImageKeywordlist outputKwl = outputPtr->GetImageKeywordlist();
    std::vector<std::string> vecEstimatedComponents = m_FunctionOnPolygon->GetEstimatedComponents();
	
    for (unsigned int i = 0; i < vecEstimatedComponents.size(); i++)
      {
	std::string key = BandKey + vecEstimatedComponents.at(i);
	outputKwl.AddKey(key, std::to_string(i)); 
      }
    outputPtr->SetImageKeywordList(outputKwl);

    // Change Projection to fit with SAR image
    outputPtr->SetProjectionRef(m_SarImagePtr->GetProjectionRef());
  }

  /** 
   * Method OutputRegionToInputRegion for GenerateInputRequestedRegion
   */
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  typename SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction >::ImageInRegionType 
  SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction >
  ::OutputRegionToInputRegion(const ImageOutRegionType& outputRegion) 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
    ImageOutSizeType  outputRequestedRegionSize = outputRegion.GetSize();
    ImageOutIndexType outputRequestedRegionIndex = outputRegion.GetIndex();
   
    // Here we are breaking traditional pipeline steps because we need to access the input field data
    // so as to compute the input image requested region
    // The inverse localisation (SAR to DEM) is performed in order to determine the DEM region (input) 
    // represented by the SAR region (output)

    int nbColDEM = m_DemImagePtr->GetLargestPossibleRegion().GetSize()[0];
    int nbLinesDEM = m_DemImagePtr->GetLargestPossibleRegion().GetSize()[1];   
 
    // Margin to apply on inverse localisation for the four sides.
    int margin = 100;
    if (std::max(nbColDEM, nbLinesDEM) > 2000)
      {
	margin = std::max(nbColDEM, nbLinesDEM)*0.01; // %1 of the largest dimension
      }
    if (margin < 100)
      {
	margin = 100;
      }

    // Indice for the SAR bloc (determined by the Pipeline)
    ImageOutIndexType id[4] ;
    if (m_OutGeometry == "SAR")
      {
	id[0][0] = outputRequestedRegionIndex[0];
	id[0][1] = outputRequestedRegionIndex[1];
	id[1][0] = outputRequestedRegionIndex[0];
	id[1][1] = outputRequestedRegionIndex[1] + outputRequestedRegionSize[1];
	id[2][0] = outputRequestedRegionIndex[0] + outputRequestedRegionSize[0];
	id[2][1] = outputRequestedRegionIndex[1];
	id[3][0] = outputRequestedRegionIndex[0] + outputRequestedRegionSize[0];
	id[3][1] = outputRequestedRegionIndex[1] + outputRequestedRegionSize[1];
      }
    else if (m_OutGeometry == "ML")
      {
	id[0][0] = outputRequestedRegionIndex[0]*m_MLRan;
	id[0][1] = outputRequestedRegionIndex[1]*m_MLAzi;
	id[1][0] = outputRequestedRegionIndex[0]*m_MLRan;
	id[1][1] = (outputRequestedRegionIndex[1] + outputRequestedRegionSize[1])*m_MLAzi;
	id[2][0] = (outputRequestedRegionIndex[0] + outputRequestedRegionSize[0])*m_MLRan;
	id[2][1] = outputRequestedRegionIndex[1]*m_MLAzi;
	id[3][0] = (outputRequestedRegionIndex[0] + outputRequestedRegionSize[0])*m_MLRan;
	id[3][1] = (outputRequestedRegionIndex[1] + outputRequestedRegionSize[1])*m_MLAzi;

	// Check Size
	if (id[1][1] > m_NbLinesSAR)
	  {
	    id[1][1] = m_NbLinesSAR;
	    id[3][1] = m_NbLinesSAR;
	  }
	if (id[2][0] > m_NbColSAR)
	  {
	    id[2][0] = m_NbColSAR;
	    id[3][0] = m_NbColSAR;
	  }
      }

    Point2DType pixelSAR_Into_MNT_LatLon(0);
    Point2DType pixelSAR(0);
    //  ImageType::IndexType pixelSAR_Into_MNT;
    itk::ContinuousIndex<double,2> pixelSAR_Into_MNT;

    // Initialie the first and last DEM indice
    int firstL, firstC, lastL, lastC; 
  
    firstL = nbLinesDEM;
    lastL = 0;
    firstC = nbColDEM;
    lastC = 0;


    // For each side of the output region
    for (int i = 0; i < 4; i++) 
      {
	// Index to physical 
	m_SarImagePtr->TransformIndexToPhysicalPoint(id[i], pixelSAR);
	// Call TransformPoint with physical point
	pixelSAR_Into_MNT_LatLon = m_RSTransform->TransformPoint(pixelSAR);
	// Transform Lat/Lon to index 
	m_DemImagePtr->TransformPhysicalPointToContinuousIndex(pixelSAR_Into_MNT_LatLon, pixelSAR_Into_MNT);
      
	// Assigne the limits
	if (firstL > pixelSAR_Into_MNT[1])
	  {
	    firstL = pixelSAR_Into_MNT[1];
	  }
      
	if (lastL < pixelSAR_Into_MNT[1])
	  {
	    lastL = pixelSAR_Into_MNT[1];
	  }

	if (firstC > pixelSAR_Into_MNT[0])
	  {
	    firstC = pixelSAR_Into_MNT[0];
	  }
      
	if (lastC < pixelSAR_Into_MNT[0])
	  {
	    lastC = pixelSAR_Into_MNT[0];
	  }

      }

    // Apply the marge
    firstL -= margin;
    firstC -= margin;
    lastL +=  margin;
    lastC +=  margin;

    // Check the limits
    if (firstC < 0)
      {
	firstC = 0;
      } 
    if (firstL < 0)
      {
	firstL = 0;
      } 
    if (lastC > nbColDEM-1)
      {
	lastC = nbColDEM-1;
      } 
    if (lastL > nbLinesDEM-1)
      {
	lastL = nbLinesDEM-1;
      } 


    // Transform sides to region
    ImageInIndexType inputRequestedRegionIndex; 
    inputRequestedRegionIndex[0] = static_cast<ImageInIndexValueType>(firstC);
    inputRequestedRegionIndex[1] = static_cast<ImageInIndexValueType>(firstL);
    ImageInSizeType inputRequestedRegionSize;
    inputRequestedRegionSize[0] = static_cast<ImageInIndexValueType>(lastC - firstC)+1;
    inputRequestedRegionSize[1] = static_cast<ImageInIndexValueType>(lastL - firstL)+1;  

      
    ImageInRegionType inputRequestedRegion = outputRegion;
    inputRequestedRegion.SetIndex(inputRequestedRegionIndex);
    inputRequestedRegion.SetSize(inputRequestedRegionSize);

    return inputRequestedRegion;  
  }

  /** 
   * Method GenerateInputRequestedRegion
   */
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void
  SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction >
  ::GenerateInputRequestedRegion()
  {
    // call the superclass' implementation of this method
    Superclass::GenerateInputRequestedRegion();
  
    ///////////// Set Output Requestion Region (Strips) ///////////////
    ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
    
    ///////////// With the new output requested region, find the region into Projeted DEM /////////////
    ImageInRegionType inputRequestedRegion = OutputRegionToInputRegion(outputRequestedRegion);
    ImageInPointer  inputPtr = const_cast< ImageInType * >( this->GetInput() );
    inputPtr->SetRequestedRegion(inputRequestedRegion);
  }


/**
   * Method PolygonAndSarLineIntersection
   */
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  bool 
  SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction >
  ::PolygonAndSarLineIntersection(const int ind_LineSAR, const double ligUL, const double colUL, 
				  const double ligUR, const double colUR, 
				  const double ligLR, const double colLR, const double ligLL, const double colLL,
				  int & i_InSideUp, int & i_InSideDown,
				  int & i_OutSideUp, int & i_OutSideDown)
  {
    double Lmax_ulur = std::max(ligUL, ligUR);
    double Lmin_ulur = std::min(ligUL, ligUR);
 
    double Lmax_urlr = std::max(ligLR, ligUR);
    double Lmin_urlr = std::min(ligLR, ligUR);
  
    double Lmax_lllr = std::max(ligLL, ligLR);
    double Lmin_lllr = std::min(ligLL, ligLR);
  
    double Lmax_ulll = std::max(ligUL, ligLL);
    double Lmin_ulll = std::min(ligUL, ligLL);
  
    // Selection for ind_LineSAR intersection
    if ( (Lmax_ulur >= ind_LineSAR && Lmin_ulur <= ind_LineSAR) ||
	 (Lmax_urlr >= ind_LineSAR && Lmin_urlr <= ind_LineSAR) ||
	 (Lmax_lllr >= ind_LineSAR && Lmin_lllr <= ind_LineSAR) ||
	 (Lmax_ulll >= ind_LineSAR && Lmin_ulll <= ind_LineSAR))
      {
	int countSides = 0;
	double CLZY_Up [4][2];
	double CLZY_Down [4][2];
	int i_Up [4];
	int i_Down [4];

	// Select if two sides is intersected by the current lineSar (ind_LineSAR)
	if (Lmax_ulur >= ind_LineSAR && Lmin_ulur <= ind_LineSAR)
	  {
	    if (ligUR > ligUL)
	      {
		CLZY_Up[countSides][0] = colUR;
		CLZY_Up[countSides][1] = ligUR;
		i_Up[countSides] = 1;
		CLZY_Down[countSides][0] = colUL;
		CLZY_Down[countSides][1] = ligUL;
		i_Down[countSides] = 0;
	      }
	    else
	      {
		CLZY_Up[countSides][0] = colUL;
		CLZY_Up[countSides][1] = ligUL;
		i_Up[countSides] = 0;
		CLZY_Down[countSides][0] = colUR;
		CLZY_Down[countSides][1] = ligUR;
		i_Down[countSides] = 1;
	      }
	
	    ++countSides;
	  }
		   
	if (Lmax_urlr >= ind_LineSAR && Lmin_urlr <= ind_LineSAR)
	  {
		       
	    if (ligLR > ligUR)
	      {
		CLZY_Up[countSides][0] = colLR;
		CLZY_Up[countSides][1] = ligLR;
		i_Up[countSides] = 2;
		CLZY_Down[countSides][0] = colUR;
		CLZY_Down[countSides][1] = ligUR;
		i_Down[countSides] = 1;
	      }
	    else
	      {
		CLZY_Up[countSides][0] = colUR;
		CLZY_Up[countSides][1] = ligUR;
		i_Up[countSides] = 1;
		CLZY_Down[countSides][0] = colLR;
		CLZY_Down[countSides][1] = ligLR;
		i_Down[countSides] = 2;
	      }
		       
	    ++countSides;
	  }
		   

	if (Lmax_lllr >= ind_LineSAR && Lmin_lllr <= ind_LineSAR)
	  {
	    if (ligLR > ligLL)
	      {
		CLZY_Up[countSides][0] = colLR;
		CLZY_Up[countSides][1] = ligLR;
		i_Up[countSides] = 2;
		CLZY_Down[countSides][0] = colLL;
		CLZY_Down[countSides][1] = ligLL;
		i_Down[countSides] = 3;
	      }
	    else
	      {
		CLZY_Up[countSides][0] = colLL;
		CLZY_Up[countSides][1] = ligLL;
		i_Up[countSides] = 3;
		CLZY_Down[countSides][0] = colLR;
		CLZY_Down[countSides][1] = ligLR;
		i_Down[countSides] = 2;
	      }

	    ++countSides;
	  }
		   
	if (Lmax_ulll >= ind_LineSAR && Lmin_ulll <= ind_LineSAR)
	  {
		       
	    if (ligUL > ligLL)
	      {
		CLZY_Up[countSides][0] = colUL;
		CLZY_Up[countSides][1] = ligUL;
		i_Up[countSides] = 0;
		CLZY_Down[countSides][0] = colLL;
		CLZY_Down[countSides][1] = ligLL;
		i_Down[countSides] = 3;
	      }
	    else
	      {
		CLZY_Up[countSides][0] = colLL;
		CLZY_Up[countSides][1] = ligLL;
		i_Up[countSides] = 3;
		CLZY_Down[countSides][0] = colUL;
		CLZY_Down[countSides][1] = ligUL;
		i_Down[countSides] = 0;
	      }
	    ++countSides;
	  }

	// If good Intesection (select input and output side)
	if (countSides ==  2 || countSides == 3)
	  {
	    i_InSideUp = i_Up[0];
	    i_InSideDown  = i_Down[0];
	    i_OutSideUp = i_Up[1];
	    i_OutSideDown = i_Down[1];
	
	    if (countSides == 3)
	      {
		// Select extermes colunms among the three CLZY points to cover the biggest part of SAR Image
		double Cmax_0 = std::max(CLZY_Up[0][0], CLZY_Down[0][0]);
		double Cmin_0 = std::min(CLZY_Up[0][0], CLZY_Down[0][0]);
		double Cmax_1 = std::max(CLZY_Up[1][0], CLZY_Down[1][0]);
		double Cmin_1 = std::min(CLZY_Up[1][0], CLZY_Down[1][0]);
		double Cmax_2 = std::max(CLZY_Up[2][0], CLZY_Down[2][0]);
		double Cmin_2 = std::min(CLZY_Up[2][0], CLZY_Down[2][0]);
	      
		// Max colunm
		if (Cmax_0 > Cmax_1 && Cmax_0 > Cmax_2)
		  {
		    i_InSideUp = i_Up[0];
		    i_InSideDown  = i_Down[0];

		    // Min colunm
		    if (Cmin_0 < Cmin_1 && Cmin_0 < Cmin_2)
		      {
			if (Cmin_1 < Cmin_2)
			  {
			    i_OutSideUp = i_Up[1];
			    i_OutSideDown  = i_Down[1];
			  }
			else
			  {
			    i_OutSideUp = i_Up[2];
			    i_OutSideDown  = i_Down[2];
			  }
		      }
		    else if (Cmin_1 < Cmin_0 && Cmin_1 < Cmin_2)
		      {
			i_OutSideUp = i_Up[1];
			i_OutSideDown  = i_Down[1];
		      }
		    else
		      {
			i_OutSideUp = i_Up[2];
			i_OutSideDown  = i_Down[2];
		      }
		  } 
		else if (Cmax_1 > Cmax_0 && Cmax_1 > Cmax_2)
		  {
		    i_InSideUp = i_Up[1];
		    i_InSideDown  = i_Down[1];
		    
		    // Min colunm
		     if (Cmin_0 < Cmin_1 && Cmin_0 < Cmin_2)
		      {
			i_OutSideUp = i_Up[0];
			i_OutSideDown  = i_Down[0];	
		      }
		     else if (Cmin_1 < Cmin_0 && Cmin_1 < Cmin_2)
		       {
			 if (Cmin_0 < Cmin_2)
			  {
			    i_OutSideUp = i_Up[0];
			    i_OutSideDown  = i_Down[0];
			  }
			else
			  {
			    i_OutSideUp = i_Up[2];
			    i_OutSideDown  = i_Down[2];
			  }
		       }
		     else
		       {
			 i_OutSideUp = i_Up[2];
			 i_OutSideDown  = i_Down[2];
		       }
		  }
		else
		  {
		    i_InSideUp = i_Up[2];
		    i_InSideDown  = i_Down[2];

		    // Min colunm
		    if (Cmin_0 < Cmin_1 && Cmin_0 < Cmin_2)
		      {
			i_OutSideUp = i_Up[0];
			i_OutSideDown  = i_Down[0];	
		      }
		     else if (Cmin_1 < Cmin_0 && Cmin_1 < Cmin_2)
		       {
			 i_OutSideUp = i_Up[1];
			 i_OutSideDown  = i_Down[1];
		       }
		     else
		       {
		    	if (Cmin_0 < Cmin_1)
			  {
			    i_OutSideUp = i_Up[0];
			    i_OutSideDown  = i_Down[0];
			  }
			else
			  {
			    i_OutSideUp = i_Up[1];
			    i_OutSideDown  = i_Down[1];
			  }
		       }
		   }
	      }
	    // Intersection with two sides or three sides
	    return true;
	    }
	else
	  {
	    // Intersection with one or four sides
	    return false;
	    }
      }
    // No intersection
    else
      {
	return false;
      }
   
  }

  /**
   * Method DEMPolygonsAnalysis
   */
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void 
  SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction >
  ::PolygonsScan(const int ind_LineSAR, 
		 std::vector<std::vector<ImageInPixelType> *> * CLZY_InSideUp_Vec, 
		 std::vector<std::vector<ImageInPixelType> *> * CLZY_InSideDown_Vec, 
		 std::vector<std::vector<ImageInPixelType> *> * CLZY_OutSideUp_Vec, 
		 std::vector<std::vector<ImageInPixelType> *> * CLZY_OutSideDown_Vec, 
		 int firstCol_into_outputRegion, int nbCol_into_outputRegion, 
		 ImageOutPixelType * outValueTab, int threadId)
  {
    bool isFirstPolygon = true;
    
    // Elevation angle max (for shadows) 
    double tgElevationMaxForCurrentLine = 0;

    int nbLinesIntoDEM_withSelectedPolygon = CLZY_InSideUp_Vec->size();

    int ind_LDEM, ind_CDEM;

    // Loop on polygons (for lines into DEM with al least one selected polygon) 
    // the loop direction depends on m_DEMdirL
    for (int i = 0 ; i < nbLinesIntoDEM_withSelectedPolygon; i++) 
      {
	// If 1, the lines increase
	if (m_DEMdirL == 1)
	  {
	    ind_LDEM = i;
	  }
	// If -1, the lines decrease
	else
	  {
	    ind_LDEM = nbLinesIntoDEM_withSelectedPolygon - (i+1);
	  } 

	// for colunms into DEM with selected polygons, depends on m_DEMdirC
	int nbColIntoDEM_withSelectedPolygon =  CLZY_InSideUp_Vec->at(ind_LDEM)->size();
	for (int j = 0; j <  nbColIntoDEM_withSelectedPolygon; j++)
	  {
		// If 1, the columns increase
		if (m_DEMdirC == 1)
		  {
		    ind_CDEM = j;
		  }
		// If -1, the colunms decrease
		else
		  {
		    ind_CDEM = nbColIntoDEM_withSelectedPolygon - (j+1);
		  }
		
		// Select the CLY_InSideUp, CLY_InSideDown, CLY_OutSideUp and CLY_OutSideDown and
		// estimate the wanted value (for example amplitude)
		m_FunctionOnPolygon->estimate(ind_LineSAR,
					      &CLZY_InSideUp_Vec->at(ind_LDEM)->at(ind_CDEM), 
					      &CLZY_InSideDown_Vec->at(ind_LDEM)->at(ind_CDEM), 
					      &CLZY_OutSideUp_Vec->at(ind_LDEM)->at(ind_CDEM), 
					      &CLZY_OutSideDown_Vec->at(ind_LDEM)->at(ind_CDEM), 
					      firstCol_into_outputRegion, 
					      nbCol_into_outputRegion, outValueTab, 
					      isFirstPolygon, tgElevationMaxForCurrentLine, threadId);

		if (isFirstPolygon)
		  {
		    isFirstPolygon = false;
		  }
	  }
      }
}
 
  /**
   * Method ThreadedGenerateData
   */
  template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
  void
  SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction >
  ::ThreadedGenerateData(const ImageOutRegionType & outputRegionForThread,
			 itk::ThreadIdType threadId)
  {
    // Compute corresponding input region
    ImageInRegionType inputRegionForThread = OutputRegionToInputRegion(outputRegionForThread);
    
    // Iterator on output (SAR geometry)
    OutputIterator OutIt(this->GetOutput(), outputRegionForThread);
  
    OutIt.GoToBegin();
   
    const ImageOutSizeType &  outputRegionForThreadSize = outputRegionForThread.GetSize();
    unsigned int nbCol = static_cast<unsigned int>(outputRegionForThreadSize[0]); 
    unsigned int firstCol = static_cast<unsigned int>(outputRegionForThread.GetIndex()[0]+this->GetOutput()->GetOrigin()[0]); 

    const ImageInSizeType & inputRegionForThreadSize = inputRegionForThread.GetSize();
    unsigned int nbDEMLines =  static_cast<unsigned int>(inputRegionForThreadSize[1]);

    
    // Allocate outValueTab (size = nbCol)
    ImageOutPixelType *outValueTab = new ImageOutPixelType[nbCol];

    // 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;
    std::vector<std::vector<ImageInPixelType> * > Polygon_SideInDown_Vec;
    std::vector<std::vector<ImageInPixelType> * > Polygon_SideOutUp_Vec;
    std::vector<std::vector<ImageInPixelType> * > Polygon_SideOutDown_Vec;

    int countDEMLine = 0;
  
    // A polygon is composed of four lines 
    // => four sides (UL UR LR LL) : UL = current_point
    //                               UR = current_point + 1 colunm
    //                               LR = current_point + 1 colunm and 1 line
    //                               LL = current_point + 1 line                                              
    // To construct all polygons of input region, four iterators are necessary 
    InputIterator  InIt_UL(this->GetInput(), inputRegionForThread);
    InputIterator  InIt_UR(this->GetInput(), inputRegionForThread);
    InputIterator  InIt_LR(this->GetInput(), inputRegionForThread);
    InputIterator  InIt_LL(this->GetInput(), inputRegionForThread);

    // Allocate Vectors for polygon storage for this line of DEM
    std::vector<ImageInPixelType> * Polygon_SideInUp_VecLine = new std::vector<ImageInPixelType>[nbDEMLines];
    std::vector<ImageInPixelType> * Polygon_SideInDown_VecLine  = new std::vector<ImageInPixelType>[nbDEMLines];
    std::vector<ImageInPixelType> * Polygon_SideOutUp_VecLine  = new std::vector<ImageInPixelType>[nbDEMLines];
    std::vector<ImageInPixelType> * Polygon_SideOutDown_VecLine  = new std::vector<ImageInPixelType>[nbDEMLines];


    // Pixel of otb::VectorImage (input) : 4 Components : C, L, Z and Y
    int i_InSideUp, i_InSideDown, i_OutSideUp, i_OutSideDown;


    InputIterator * tabInIt [4];
    tabInIt[0] = &InIt_UL;
    tabInIt[1] = &InIt_UR;
    tabInIt[2] = &InIt_LR;
    tabInIt[3] = &InIt_LL;

    double UL_line, UR_line, LR_line, LL_line; 
    double UL_col, UR_col, LR_col, LL_col;

    bool hasIntersection = false;

    // For each line of output : construct polygon after polygon, determine if polygon intersects the current 
    // line and estimate the amplitude with the contribution of all selected polygon
    while ( !OutIt.IsAtEnd())
      {
	OutIt.GoToBeginOfLine();
	// indice of current line (into output Geometry)
	//~ int ind_Line = OutIt.GetIndex()[1] + int(this->GetOutput()->GetOrigin()[1]);
     
	// The Projeted DEM Pixel/Bands is always into SAR geometry => must be adapted for mlRan and mlAzi != 1
	// If mlRan and mlAzi == 1 same as SAR geometry (loop on ind_Line_start=ind_Line only)
	int ind_Line_start = OutIt.GetIndex()[1] * m_MLAzi + int(this->GetOutput()->GetOrigin()[1]);
	int ind_Line_end =  (OutIt.GetIndex()[1]+1) * m_MLAzi + int(this->GetOutput()->GetOrigin()[1]);
	if (ind_Line_end > (m_NbLinesSAR + int(this->GetOutput()->GetOrigin()[1])))
	  {
	    ind_Line_end = m_NbLinesSAR + int(this->GetOutput()->GetOrigin()[1]);
	  }

	// Initialize outValueTab
	m_FunctionOnPolygon->initalize(nbCol, outValueTab, threadId);

	// Loop on Line into SAR Geometry
	for (int ind_Line_Into_SARGeo = ind_Line_start; ind_Line_Into_SARGeo < ind_Line_end; 
	     ind_Line_Into_SARGeo++)
	  {
	    ///////////////////////////// Polygon Construction and Check intersection /////////////////////////
	    InIt_UL.GoToBegin();
	    InIt_UR.GoToBegin();
	    InIt_LR.GoToBegin();
	    InIt_LL.GoToBegin();
     
	    // Put each Iterator to this start line
	    // +1 line for LL and for LR
	    InIt_LL.NextLine();
	    InIt_LR.NextLine();
     
	    countDEMLine = 0;
     
	    // For each Line of input 
	    while ( !InIt_LR.IsAtEnd() && !InIt_LL.IsAtEnd())
	      {
		InIt_UR.GoToBeginOfLine();
		InIt_UL.GoToBeginOfLine();
		InIt_LR.GoToBeginOfLine();
		InIt_LL.GoToBeginOfLine();
	 
		// Put each Iterator to this start colunm
		// +1 colunm for UR and LR
		++InIt_UR;
		++InIt_LR;
	 
		// For each column
		while (!InIt_UR.IsAtEndOfLine() && !InIt_LR.IsAtEndOfLine())
		  {
		    UL_line = InIt_UL.Get()[1];

		    // Protect Memory access and input access
		    if (UL_line > ind_Line_Into_SARGeo - m_Margin && UL_line < ind_Line_Into_SARGeo + m_Margin)
		      { 	 
			UL_col = InIt_UL.Get()[0];
			UR_col = InIt_UR.Get()[0];
			LR_col = InIt_LR.Get()[0];
			LL_col = InIt_LL.Get()[0];
		
			// Check if No Data		    
			if (!(UL_col == m_NoData || UR_col == m_NoData || LR_col == m_NoData 
			      || LL_col == m_NoData))
			  {
			    // Determine if intersection
			    UR_line = InIt_UR.Get()[1];
			    LR_line = InIt_LR.Get()[1];
			    LL_line = InIt_LL.Get()[1];
		
			    // With lig and col of UL, UR, LR and LL select iterators for InSide and Outside of 
			    // each polygon (if intersection)
			    hasIntersection = this->PolygonAndSarLineIntersection(ind_Line_Into_SARGeo, 
										  UL_line, UL_col, 
										  UR_line, UR_col, 
										  LR_line, LR_col, 
										  LL_line, LL_col, 
										  i_InSideUp, 
										  i_InSideDown,
										  i_OutSideUp, 
										  i_OutSideDown);
			    
			    // If intersection == true : Retrive the intersected sides (two intersected among 
			    // the four sides of the polygon)
			    if (hasIntersection)
			      {
				// Store the polygons
				Polygon_SideInUp_VecLine[countDEMLine].push_back(tabInIt[i_InSideUp]->Get());