-
Gaëlle USSEGLIO authoredGaëlle USSEGLIO authored
otbSARDEMPolygonsAnalysisImageFilter.txx 39.80 KiB
/*
* 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>
#define sign(a,b) (((b)>=0)?fabs((double)(a)):-fabs((double)(a)))
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)
{
// Outputs required and/or needed (one required and one optionnal)
this->SetNumberOfRequiredOutputs(1);
this->SetNumberOfIndexedOutputs(2);
this->SetNthOutput(0, ImageOutType::New());
this->SetNthOutput(1, ImageOptionnalType::New());
m_OutputCounter = 0;
}
template <class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
TImageOut *
SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >
::GetOutput()
{
if (this->GetNumberOfOutputs() < 1)
{
return ITK_NULLPTR;
}
return static_cast<ImageOutType *>(this->itk::ProcessObject::GetOutput(0));
}
template <class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
const TImageOut * SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >
::GetOutput() const
{
if (this->GetNumberOfOutputs() < 1)
{
return 0;
}
return static_cast<ImageOutType *>(this->itk::ProcessObject::GetOutput(0));
}
template <class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
typename SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >::ImageOptionnalType *
SARDEMPolygonsAnalysisImageFilter< TImageIn ,TImageOut, TImageDEM, TImageSAR, TFunction >
::GetOptionnalOutput()
{
if (this->GetNumberOfOutputs() < 2)
{
return ITK_NULLPTR;
}
return static_cast<ImageOptionnalType *>(this->itk::ProcessObject::GetOutput(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)
{
// Check if sarImageKWL not NULL
assert(&sarImageKWL && "SAR Image Metadata don't exist.");
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 > nbLinesDEM)
{
m_Margin = (nbLinesSAR/nbLinesDEM)*5;
}
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 (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->GetSpacing();
// Define Output Largest Region
ImageOutRegionType outputLargestPossibleRegion = m_SarImagePtr->GetLargestPossibleRegion();
outputLargestPossibleRegion.SetSize(outputSize);
outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
outputPtr->SetOrigin(outOrigin);
outputPtr->SetSpacing(outSP);
// 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->GetSpacing();
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->SetSpacing(outSP);
// 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 (int i = 0; i < vecEstimatedComponents.size(); i++)
{
std::string key = BandKey + vecEstimatedComponents.at(i);
outputKwl.AddKey(key, std::to_string(i));
}
outputPtr->SetImageKeywordList(outputKwl);
////////// Check if the optionnal output is required //////////
if (m_FunctionOnPolygon->GetWithOptionnalOutput())
{
// This output is just a vector with 1xNbLine_MainOutput pixels.
// Pixels are the same nature than Main Output
ImageOptionnalPointer outputOptionnalPtr = this->GetOptionnalOutput();
if (numberOfEstimatedComponents > 1)
{
// Set the number of Components
outputOptionnalPtr->SetNumberOfComponentsPerPixel(numberOfEstimatedComponents);
}
ImageOptionnalRegionType region = outputPtr->GetLargestPossibleRegion();
ImageOptionnalSizeType sizeRegion;
sizeRegion[0] = 1;
sizeRegion[1] = outputPtr->GetLargestPossibleRegion().GetSize()[1];
region.SetSize(sizeRegion);
outputOptionnalPtr->SetLargestPossibleRegion(region);
outputOptionnalPtr->SetOrigin(outputPtr->GetOrigin());
outputOptionnalPtr->SetSpacing(outputPtr->GetSpacing());
// Complete Optionnal Output KeyWordList
ImageKeywordlist outputOptionnalKwl;
for (int i = 0; i < vecEstimatedComponents.size(); i++)
{
std::string key = BandKey + vecEstimatedComponents.at(i);
outputOptionnalKwl.AddKey(key, std::to_string(i));
}
outputOptionnalPtr->SetImageKeywordList(outputOptionnalKwl);
if (m_OutputCounter == 0)
{
// Allocate optionnal Results (Vector)
int nbLineOptionnal = this->GetOptionnalOutput()->GetLargestPossibleRegion().GetSize()[1];
m_OptionnalResults = new ImageOutPixelType[nbLineOptionnal];
}
}
// Increment Output Counter
++m_OutputCounter;
}
/**
* Method EnlargeOutputRequestedRegion
*/
template<class TImageIn, class TImageOut, class TImageDEM, class TImageSAR, class TFunction>
void
SARDEMPolygonsAnalysisImageFilter< TImageIn, TImageOut, TImageDEM, TImageSAR, TFunction >
::EnlargeOutputRequestedRegion( itk::DataObject *itkNotUsed(output) )
{
// This filter requires all of the optionnal output image (Nth Output = 1). The output is // vector with 1xNbLines_Main_Output as size
if (m_FunctionOnPolygon->GetWithOptionnalOutput())
{
if ( this->itk::ProcessObject::GetOutput(1) )
{
this->itk::ProcessObject::GetOutput(1)->SetRequestedRegionToLargestPossibleRegion();
}
}
}
/**
* 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
}
// 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();
if (m_FunctionOnPolygon->GetWithOptionnalOutput())
{
if (outputRequestedRegion.GetIndex()[0] != 0 &&
outputRequestedRegion.GetSize()[0] != this->GetOutput()->GetLargestPossibleRegion().GetSize()[0])
{
ImageOutIndexType outputIndex;
outputIndex[0] = 0;
outputIndex[1] = outputRequestedRegion.GetIndex()[1];
ImageOutSizeType outputSize;
outputSize[0] = this->GetOutput()->GetLargestPossibleRegion().GetSize()[0];
outputSize[1] = outputRequestedRegion.GetSize()[1];
// Region affectation (with Crop to respect image dimension)
outputRequestedRegion.SetSize(outputSize);
outputRequestedRegion.SetIndex(outputIndex);
//outputRequestedRegion.Crop(this->GetOutput()->GetLargestPossibleRegion());
this->GetOutput()->SetRequestedRegion(outputRequestedRegion);
}
}
///////////// 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>
bool
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 >
::BeforeThreadedGenerateData()
{
// Allocation for the optionnal image (if needed)
if (m_FunctionOnPolygon->GetWithOptionnalOutput())
{
// Allocate optionnal output
ImageOptionnalPointer outputOptionnalPtr = this->GetOptionnalOutput();
outputOptionnalPtr->SetBufferedRegion(outputOptionnalPtr->GetRequestedRegion() );
outputOptionnalPtr->Allocate();
}
}
/**
* 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[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 [nbDEMLines];
std::vector<ImageInPixelType> Polygon_SideInDown_VecLine [nbDEMLines];
std::vector<ImageInPixelType> Polygon_SideOutUp_VecLine [nbDEMLines];
std::vector<ImageInPixelType> Polygon_SideOutDown_VecLine [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 = ind_Line * m_MLAzi;
int ind_Line_end = (ind_Line+1) * m_MLAzi;
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());
Polygon_SideInDown_VecLine[countDEMLine].push_back(tabInIt[i_InSideDown]->Get());
Polygon_SideOutUp_VecLine[countDEMLine].push_back(tabInIt[i_OutSideUp]->Get());
Polygon_SideOutDown_VecLine[countDEMLine].push_back(tabInIt[i_OutSideDown]->Get());
}
}
}
// Next Column For input
++InIt_UL;
++InIt_UR;
++InIt_LR;
++InIt_LL;
} // End colunms (input)
// if vector of for the current Line of DEM is not empty => store this vector into our vectors
// of vectors
if ( !(Polygon_SideInUp_VecLine[countDEMLine].empty()))
{
Polygon_SideInUp_Vec.push_back(&Polygon_SideInUp_VecLine[countDEMLine]);
Polygon_SideInDown_Vec.push_back(&Polygon_SideInDown_VecLine[countDEMLine]);
Polygon_SideOutUp_Vec.push_back(&Polygon_SideOutUp_VecLine[countDEMLine]);
Polygon_SideOutDown_Vec.push_back(&Polygon_SideOutDown_VecLine[countDEMLine]);
}
++countDEMLine;
// Next Line For input
InIt_UL.NextLine();
InIt_UR.NextLine();
InIt_LR.NextLine();
InIt_LL.NextLine();
} // End of lines (input)
///////////////////////////// Scan of all selected polygons to extract wanted data //////////////
this->PolygonsScan(ind_Line_Into_SARGeo, &Polygon_SideInUp_Vec,
&Polygon_SideInDown_Vec, &Polygon_SideOutUp_Vec,
&Polygon_SideOutDown_Vec, firstCol, nbCol, outValueTab, threadId);
// Clear Vectors
for (int k = 0; k < nbDEMLines; k++)
{
Polygon_SideInUp_VecLine[k].clear();
Polygon_SideInDown_VecLine[k].clear();
Polygon_SideOutUp_VecLine[k].clear();
Polygon_SideOutDown_VecLine[k].clear();
}
Polygon_SideInUp_Vec.clear();
Polygon_SideInDown_Vec.clear();
Polygon_SideOutUp_Vec.clear();
Polygon_SideOutDown_Vec.clear();
} // End Loop on Line into SAR Geometry
// After the scan : Synthetize results
m_FunctionOnPolygon->synthetize(nbCol, outValueTab, threadId);
// Assignate the output with the contribution of the current polygons
int counter = 0;
while (!OutIt.IsAtEndOfLine())
{
// Affect OutIt
OutIt.Set(outValueTab[counter]);
++counter;
++OutIt;
} // End colunms (ouput)
// Next Line
OutIt.NextLine();
/// Optionnal Output ////
if (m_FunctionOnPolygon->GetWithOptionnalOutput())
{
// Assignate the optionnal image
ImageOutIndexType indexOptionnalImage;
indexOptionnalImage[0] = 0; // Always 0 since optionnal image is a vector
indexOptionnalImage[1] = ind_Line;
// Estimate the value for the current ind_Line and store it into m_OptionnalResults
ImageOutPixelType outOptionnalValue;
// If first output => Estimate ans store results if second => copy and write into optionnal image
// m_OptionnalResults is threadSafe because only one thread access to ind_Line index
if (m_OutputCounter == 1)
{
m_FunctionOnPolygon->estimateOptionnalImage(&m_OptionnalResults[ind_Line], threadId);
}
else
{
this->GetOptionnalOutput()->SetPixel(indexOptionnalImage, m_OptionnalResults[ind_Line]);
}
}
} // End lines (Main output)
}
} /*namespace otb*/
#endif