Skip to content
Snippets Groups Projects
Commit 5cdadc37 authored by Gaëlle USSEGLIO's avatar Gaëlle USSEGLIO
Browse files

UPDATE : Changing TopographicPhase Filter with the new OTB function : LineToSatPositionAndVelocity

parent 2f398c6a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
include/otbSARTopographicPhaseImageFilter.txx deleted 100644 → 0
+ 0
− 547
View file @ 2f398c6a
/*
* 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 otbSARTopographicPhaseImageFilter_txx
#define otbSARTopographicPhaseImageFilter_txx
#include "otbSARTopographicPhaseImageFilter.h"
#include "itkImageScanlineConstIterator.h"
#include "itkImageScanlineIterator.h"
#include "itkProgressReporter.h"
#include "itkNumericTraitsPointPixel.h"
#include "ossim/ossimSarSensorModel.h"
#include <cmath>
#include <algorithm>
#include <omp.h>
namespace otb
{
/**
* Constructor with default initialization
*/
template <class TImageIn, class TImageOut>
SARTopographicPhaseImageFilter< TImageIn, TImageOut >::SARTopographicPhaseImageFilter()
: m_MLRan(1), m_MLAzi(1), m_MaxShiftInRange(0), m_MaxShiftInAzimut(0), m_Factor(2),
m_MasterCopy(false), m_SarSensorModelAdapterForSlave(ITK_NULLPTR),
m_SarSensorModelAdapterForMaster(ITK_NULLPTR), m_OutputCounter(0), m_ApproxDiapason(false)
{
// Inputs required and/or needed
this->SetNumberOfRequiredInputs(2);
m_GridStep.Fill(1);
}
/**
* Destructor
*/
template <class TImageIn, class TImageOut>
SARTopographicPhaseImageFilter< TImageIn, TImageOut >::~SARTopographicPhaseImageFilter()
{
}
/**
* Print
*/
template<class TImageIn, class TImageOut>
void
SARTopographicPhaseImageFilter< TImageIn, 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 << "Grid Step : " << m_GridStep << std::endl;
}
/**
* Set Master Cartesian Mean Image
*/
template<class TImageIn, class TImageOut>
void
SARTopographicPhaseImageFilter< TImageIn, TImageOut >
::SetMasterCartesianMeanInput(const ImageInType* image )
{
// Process object is not const-correct so the const casting is required.
this->SetNthInput(0, const_cast<ImageInType *>(image));
}
/**
* Set Shift Grid
*/
template<class TImageIn, class TImageOut>
void
SARTopographicPhaseImageFilter< TImageIn, TImageOut >
::SetGridInput(const GridType* image )
{
// Process object is not const-correct so the const casting is required.
this->SetNthInput(1, const_cast<GridType *>(image));
}
/**
* Get Master Cartesian Mean Image
*/
template<class TImageIn, class TImageOut>
const typename SARTopographicPhaseImageFilter< TImageIn, TImageOut >::ImageInType *
SARTopographicPhaseImageFilter< TImageIn, TImageOut >
::GetMasterCartesianMeanInput() const
{
if (this->GetNumberOfInputs()<1)
{
return 0;
}
return static_cast<const ImageInType *>(this->itk::ProcessObject::GetInput(0));
}
/**
* Get Shift Grid
*/
template<class TImageIn, class TImageOut>
const typename SARTopographicPhaseImageFilter< TImageIn, TImageOut >::GridType *
SARTopographicPhaseImageFilter< TImageIn, TImageOut >
::GetGridInput() const
{
if (this->GetNumberOfInputs()<1)
{
return 0;
}
return static_cast<const GridType *>(this->itk::ProcessObject::GetInput(1));
}
template<class TImageIn, class TImageOut>
void
SARTopographicPhaseImageFilter< TImageIn, TImageOut >
::SetSlaveImageKeyWorList(ImageKeywordlist sarSlaveKWL)
{
m_SlaveKeyWordList = sarSlaveKWL;
}
/**
* Method GenerateOutputInformaton()
**/
template<class TImageIn, class TImageOut>
void
SARTopographicPhaseImageFilter< TImageIn, TImageOut >
::GenerateOutputInformation()
{
// Call superclass implementation
Superclass::GenerateOutputInformation();
// Get pointers to the input and output
ImageInConstPointer masterCartMeanPtr = this->GetMasterCartesianMeanInput();
ImageOutPointer outputPtr = this->GetOutput();
// KeyWordList
ImageKeywordlist masterKWL = masterCartMeanPtr->GetImageKeywordlist();
// Master SAR Dimensions
int nbColMasterCartMean = masterCartMeanPtr->GetLargestPossibleRegion().GetSize()[0];
int nbLinesMasterCartMean = masterCartMeanPtr->GetLargestPossibleRegion().GetSize()[1];
//////////////////////////// TopographicPhase into Master Cartesian Mean Geo ///////////////////////////
// Vector Image :
// At Least 5 Components :
// _ Topographic Phase
// _ IsData Mask
// _ Copy of XCart mean Master (if copy)
// _ Copy of YCart mean Master (if copy)
// _ Copy of ZCart mean Master (if copy)
if (m_MasterCopy)
{
outputPtr->SetNumberOfComponentsPerPixel(5);
}
else
{
outputPtr->SetNumberOfComponentsPerPixel(2);
}
// The output is defined with the Master Cartesian Mean
ImageOutSizeType outputSize;
outputSize[0] = nbColMasterCartMean;
outputSize[1] = nbLinesMasterCartMean;
ImageOutPointType outOrigin;
outOrigin = masterCartMeanPtr->GetOrigin();
ImageOutSpacingType outSP;
outSP = masterCartMeanPtr->GetSpacing();
// Define Output Largest Region
ImageOutRegionType outputLargestPossibleRegion = masterCartMeanPtr->GetLargestPossibleRegion();
outputLargestPossibleRegion.SetSize(outputSize);
outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
outputPtr->SetOrigin(outOrigin);
outputPtr->SetSpacing(outSP);
if (m_OutputCounter == 0)
{
////////// Checks (with input keywordlists/metadata) /////////////
// Check gridSteps
ImageKeywordlist gridKWL = this->GetGridInput()->GetImageKeywordlist();
if (gridKWL.HasKey("support_data.gridstep.range") && gridKWL.HasKey("support_data.gridstep.azimut"))
{
unsigned int gridStepRange = atoi(gridKWL.GetMetadataByKey("support_data.gridstep.range").c_str());
unsigned int gridStepAzimut = atoi(gridKWL.GetMetadataByKey("support_data.gridstep.azimut").c_str());
if (gridStepRange != m_GridStep[0] || gridStepAzimut != m_GridStep[1])
{
itkExceptionMacro(<<"Provided GridSteps are not consistent with grid keywordlist.");
}
}
// Adapt GridStep to output Geometry (ML)
if (m_MLRan != 1 || m_MLAzi != 1)
{
if (m_GridStep[0] % m_MLRan)
{
itkExceptionMacro(<<"GridSteps range mot a multiple of MLRan.");
}
else
{
m_GridStep[0] /= m_MLRan;
}
if (m_GridStep[1] % m_MLAzi)
{
itkExceptionMacro(<<"GridSteps azimut mot a multiple of MLAzi.");
}
else
{
m_GridStep[1] /= m_MLAzi;
}
}
// Check ML Factors
if (masterKWL.HasKey("support_data.ml_ran") && masterKWL.HasKey("support_data.ml_azi"))
{
// Get Master ML Factors
unsigned int master_MLRan = atoi(masterKWL.GetMetadataByKey("support_data.ml_ran").c_str());
unsigned int master_MLAzi = atoi(masterKWL.GetMetadataByKey("support_data.ml_azi").c_str());
if ((master_MLRan != m_MLRan) || (master_MLAzi != m_MLAzi))
{
itkExceptionMacro(<<"ML Factor betwwen master and inputs of this application are different.");
}
}
// Set new keyword list to output image with bands meaning and ML Factors
ImageKeywordlist outputKWL = masterKWL;
if (masterKWL.HasKey("support_data.ml_ran") && masterKWL.HasKey("support_data.ml_azi"))
{
outputKWL.AddKey("support_data.ml_ran", std::to_string(m_MLRan));
outputKWL.AddKey("support_data.ml_azi", std::to_string(m_MLAzi));
}
outputKWL.AddKey("support_data.band.Topographic", std::to_string(0));
outputKWL.AddKey("support_data.band.XCart", std::to_string(1));
outputKWL.AddKey("support_data.band.YCart", std::to_string(2));
outputKWL.AddKey("support_data.band.ZCart", std::to_string(3));
outputKWL.AddKey("support_data.band.isData", std::to_string(4));
outputPtr->SetImageKeywordList(outputKWL);
// Calculate lambda
double radarFreq = atof(masterKWL.GetMetadataByKey("support_data.radar_frequency").c_str());
const double C = 299792458.;
m_Lambda = C/radarFreq;
// Create and Initilaze SarSensorModelAdapters
m_SarSensorModelAdapterForMaster = SarSensorModelAdapter::New();
bool loadOk = m_SarSensorModelAdapterForMaster->LoadState(masterKWL);
if(!loadOk || !m_SarSensorModelAdapterForMaster->IsValidSensorModel())
{
itkExceptionMacro(<<"SAR image does not contain a valid SAR sensor model.");
}
m_SarSensorModelAdapterForSlave = SarSensorModelAdapter::New();
loadOk = m_SarSensorModelAdapterForSlave->LoadState(m_SlaveKeyWordList);
if(!loadOk || !m_SarSensorModelAdapterForSlave->IsValidSensorModel())
{
itkExceptionMacro(<<"SAR image does not contain a valid SAR sensor model.");
}
}
++m_OutputCounter;
}
/**
* Method OutputRegionToInputRegion for GenerateInputRequestedRegion
*/
template<class TImageIn, class TImageOut>
typename SARTopographicPhaseImageFilter< TImageIn, TImageOut >::GridRegionType
SARTopographicPhaseImageFilter< TImageIn, TImageOut >
::OutputRegionToInputGridRegion(const ImageOutRegionType& outputRegion) const
{
// Compute the input requested region (size and start index)
// Use the image transformations to insure an input grid requested region
// that will provide the proper range
const ImageOutSizeType & outputRequestedRegionSize = outputRegion.GetSize();
ImageOutIndexType outputRequestedRegionIndex = outputRegion.GetIndex();
// Define the index and size for grid requested region. The input grid is on output/master geometry with
// m_GridStep as factor
GridIndexType indexGrid;
indexGrid[0] = static_cast<GridIndexValueType>(outputRequestedRegionIndex[0]/m_GridStep[0]) - 1;
indexGrid[1] = static_cast<GridIndexValueType>(outputRequestedRegionIndex[1]/m_GridStep[1]) - 1;
GridSizeType sizeGrid;
sizeGrid[0] = static_cast<GridSizeValueType>(outputRequestedRegionSize[0]/m_GridStep[0]) + 3;
sizeGrid[1] = static_cast<GridSizeValueType>(outputRequestedRegionSize[1]/m_GridStep[1]) + 3;
// Check Index and Size
if (indexGrid[0] < this->GetGridInput()->GetLargestPossibleRegion().GetIndex()[0])
{
indexGrid[0] = this->GetGridInput()->GetLargestPossibleRegion().GetIndex()[0];
}
if (indexGrid[1] < this->GetGridInput()->GetLargestPossibleRegion().GetIndex()[1])
{
indexGrid[1] = this->GetGridInput()->GetLargestPossibleRegion().GetIndex()[1];
}
if ((sizeGrid[0] + indexGrid[0]) >
this->GetGridInput()->GetLargestPossibleRegion().GetSize()[0])
{
sizeGrid[0] = this->GetGridInput()->GetLargestPossibleRegion().GetSize()[0] -
indexGrid[0];
}
if ((sizeGrid[1] + indexGrid[1]) >
this->GetGridInput()->GetLargestPossibleRegion().GetSize()[1])
{
sizeGrid[1] = this->GetGridInput()->GetLargestPossibleRegion().GetSize()[1] -
indexGrid[1];
}
// Transform into a region1
GridRegionType gridRequestedRegion = outputRegion;
gridRequestedRegion.SetIndex(indexGrid);
gridRequestedRegion.SetSize(sizeGrid);
return gridRequestedRegion;
}
/**
* Method GenerateInputRequestedRegion
*/
template<class TImageIn, class TImageOut>
void
SARTopographicPhaseImageFilter< TImageIn, TImageOut >
::GenerateInputRequestedRegion()
{
// call the superclass' implementation of this method
Superclass::GenerateInputRequestedRegion();
// Get Output requested region
ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
///////////// For Master Cartesian Mean same region /////////////
ImageInPointer masterCartMeanPtr = const_cast< ImageInType * >( this->GetMasterCartesianMeanInput() );
masterCartMeanPtr->SetRequestedRegion(outputRequestedRegion);
///////////// Find the region into Shift Grid ////////////
GridRegionType gridRequestedRegion = OutputRegionToInputGridRegion(outputRequestedRegion);
GridPointer gridPtr = const_cast< GridType * >( this->GetGridInput() );
gridPtr->SetRequestedRegion(gridRequestedRegion);
}
/**
* Method ThreadedGenerateData
*/
template<class TImageIn, class TImageOut>
void
SARTopographicPhaseImageFilter< TImageIn, TImageOut >
::ThreadedGenerateData(const ImageOutRegionType & outputRegionForThread,
itk::ThreadIdType threadId)
{
// Compute corresponding input region for master and slave cartesian mean
ImageInRegionType inputMasterRegionForThread = outputRegionForThread;
// Compute corresponding input region for grid
GridRegionType inputGridRegionForThread = OutputRegionToInputGridRegion(outputRegionForThread);
// Iterator on output
OutputIterator OutIt(this->GetOutput(), outputRegionForThread);
OutIt.GoToBegin();
// Iterator on input master cartesian mean
InputIterator InMasterCartMeanIt(this->GetMasterCartesianMeanInput(), inputMasterRegionForThread);
InMasterCartMeanIt.GoToBegin();
// Allocate output pixel
ImageOutPixelType pixelOut;
if (m_MasterCopy)
{
pixelOut.Reserve(5);
}
else
{
pixelOut.Reserve(2);
}
double constMul = static_cast<double>(m_Factor*2*M_PI)/m_Lambda;
if (m_ApproxDiapason)
{
constMul = static_cast<double>(m_Factor*256)/m_Lambda;
}
Point3DType worldSlave;
Point3DType satPosSlave;
Point3DType satVelSlave;
Point3DType worldMaster;
Point3DType satPosMaster;
Point3DType satVelMaster;
// For each line
while (!OutIt.IsAtEnd() && !InMasterCartMeanIt.IsAtEnd())
{
OutIt.GoToBeginOfLine();
InMasterCartMeanIt.GoToBeginOfLine();
// Index of current line (into output Geometry)
int ind_Line = OutIt.GetIndex()[1] + int(this->GetOutput()->GetOrigin()[1]);
// Get Master Cartesian Mean Per line
ImageInIndexType indexMasterPerLine;
indexMasterPerLine[0] = 0; // Always 0 since Master/Slave Cartesain Per Line are vectors
indexMasterPerLine[1] = ind_Line;
// Get the index of current tile into grid to retrive the shifts (the closest (round) grid point
// at the center of current tile). Output Geo = Master Cart Mean Geo = (Grid geo / GridStep)
int Lgrid = std::round( ind_Line / m_GridStep[1]);
Lgrid = std::min (std::max (Lgrid, 0),
static_cast<int>(this->GetGridInput()->GetLargestPossibleRegion().GetSize()[1])-1);
GridIndexType gridIndex;
gridIndex[0] = 0;
gridIndex[1] = Lgrid;
double gridShift_Azi = this->GetGridInput()->GetPixel(gridIndex)[1];
// Apply on slave, the integer shifts
int Le = std::round(ind_Line + gridShift_Azi);
// Get Slave Cartesian Means
ImageInIndexType indexSlavePerLine;
indexSlavePerLine[0] = 0; // Always 0 since Master/Slave Cartesain Per Line are vectors
indexSlavePerLine[1] = Le;
////////// Estimate satellite positions for the current line //////////
bool checkSlave = m_SarSensorModelAdapterForSlave->LineToSatPositionAndVelocity(static_cast<double>(Le), satPosSlave,
satVelSlave);
bool checkMaster = m_SarSensorModelAdapterForMaster->LineToSatPositionAndVelocity(static_cast<double>(ind_Line),
satPosMaster, satVelMaster);
// For each colunm
while (!OutIt.IsAtEndOfLine() && !InMasterCartMeanIt.IsAtEndOfLine())
{
// Check slave Index
if (checkSlave)
{
// Check if Value into Master Cartesian Mean with IsData Mask
if (InMasterCartMeanIt.Get()[3] != 0)
{
float Xcart_Ground = InMasterCartMeanIt.Get()[0];
float Ycart_Ground = InMasterCartMeanIt.Get()[1];
float Zcart_Ground = InMasterCartMeanIt.Get()[2];
//////////// Estimate Topographic phase (P = (factor*256/lambda) * (De-Dm)) //////////
double De = sqrt(pow((Xcart_Ground - satPosSlave[0]), 2) +
pow((Ycart_Ground - satPosSlave[1]), 2) +
pow((Zcart_Ground - satPosSlave[2]), 2));
double Dm = sqrt(pow((Xcart_Ground - satPosMaster[0]), 2) +
pow((Ycart_Ground - satPosMaster[1]), 2) +
pow((Zcart_Ground - satPosMaster[2]), 2));
pixelOut[0] = constMul * (De-Dm);
// Mod 2*Pi
//pixelOut[0] = pixelOut[0]-(2*M_PI)*floor(pixelOut[0]/(2*M_PI));
// IsData set to 1
pixelOut[1] = 1;
if (m_MasterCopy)
{
////////////// Copy of Master Cartesian Mean //////////////
pixelOut[2] = Xcart_Ground;
pixelOut[3] = Ycart_Ground;
pixelOut[4] = Zcart_Ground;
}
}
else
{
// All components set to 0
pixelOut[0] = 0;
pixelOut[1] = 0;
if (m_MasterCopy)
{
pixelOut[2] = 0;
pixelOut[3] = 0;
pixelOut[4] = 0;
}
}
}
else
{
// All components set to 0
pixelOut[0] = 0;
pixelOut[1] = 0;
if (m_MasterCopy)
{
pixelOut[2] = 0;
pixelOut[3] = 0;
pixelOut[4] = 0;
}
}
//////////// Assign Output ////////////
OutIt.Set(pixelOut);
// Increment iterators
++OutIt;
++InMasterCartMeanIt;
} // End colunms (ouput)
// Next Line
OutIt.NextLine();
InMasterCartMeanIt.NextLine();
} // End lines (ouput)
}
} /*namespace otb*/
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment