Skip to content
Snippets Groups Projects

Resolve "Factorize Computations in SARDEMProjection"

Merged Luc Hermitte requested to merge 7-factorize-computations-in-sardemprojection into master
1 file
+ 3
3
Compare changes
  • Side-by-side
  • Inline
/*
* Copyright (C) 2005-2023 Centre National d'Etudes Spatiales (CNES)
* Copyright (C) 2005-2024 Centre National d'Etudes Spatiales (CNES)
*
* This file is part of Orfeo Toolbox
*
@@ -55,6 +55,17 @@
#include <thread>
#include <cassert>
#if defined(NDEBUG)
# define assert_op(a, op, b) ((void)0)
#else
# define assert_op(a, op, b) \
if (!( a op b)) { \
std::cerr << __FILE__ ":" << __LINE__ <<": assertion failed: " << a << " is not " #op " " << b << std::endl; \
std::abort(); \
}
#endif
namespace otb
{
/**
@@ -132,7 +143,7 @@ SARDEMProjectionImageFilter2< TImageIn, TImageOut >
// _ L : Line of DEM into SAR Image
// _ Z : Coordonates Z
// _ Y : Coordonates Y
// (Z,Y) coordinates, defined into ossimSarSensorMode (if otb version < 8) or SarSensorMode as follows:
// (Z,Y) coordinates, defined into ossimSarSensorModel (if otb version < 8) or SarSensorModel as follows:
// Let n = |sensorPos|,
// ps2 = scalar_product(sensorPos,worldPoint)
// d = distance(sensorPos,worldPoint)
@@ -213,7 +224,7 @@ SARDEMProjectionImageFilter2< TImageIn, TImageOut >
#if OTB_VERSION_MAJOR >= 8
// Create and Initilaze the SarSensorModel
m_SarSensorModel = std::make_unique<SarSensorModel>(m_SarImageMetadata);
m_SarSensorModel = std::make_unique<FastSarSensorModel>(m_SarImageMetadata);
#else
// Create and Initilaze the SarSensorModelAdapter
m_SarSensorModel = SarSensorModelAdapter::New();
@@ -226,11 +237,11 @@ SARDEMProjectionImageFilter2< TImageIn, TImageOut >
#endif
// Compute the output image size, spacing and origin (same as the input)
const ImageInSizeType & inputSize = inputPtr->GetLargestPossibleRegion().GetSize();
ImageOutRegionType outputLargestPossibleRegion = inputPtr->GetLargestPossibleRegion();
const ImageInSizeType & inputSize = outputLargestPossibleRegion.GetSize();
ImageInPointType origin = inputPtr->GetOrigin();
const ImageInSpacingType & inSP = inputPtr->GetSpacing();
ImageOutRegionType outputLargestPossibleRegion = inputPtr->GetLargestPossibleRegion();
outputLargestPossibleRegion.SetSize(static_cast<ImageOutSizeType>(inputSize));
outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
outputPtr->SetOrigin(static_cast<ImageOutPointType>(origin));
@@ -292,11 +303,6 @@ SARDEMProjectionImageFilter2< TImageIn, TImageOut >
Point3DType demGeoPoint;
Point2DType demLatLonPoint;
Point2DType col_row(0);
Point2DType y_z(0);
Point3DType xyzCart;
Point3DType satPos;
Point3DType satVel;
assert(m_NbComponents ==
4 + (m_withXYZ ? 3 : 0) + (m_withH ? 1 : 0) + (m_withSatPos ? 3 : 0));
@@ -366,7 +372,24 @@ SARDEMProjectionImageFilter2< TImageIn, TImageOut >
// Correct the height with earth geoid
demGeoPoint[2] = InIt.Get() + h;
// Call the method of sarSensorModel
// Always compute the ground point is ECEF
#if OTB_VERSION_MAJOR >= 8
Point3DType xyzCart = Projection::WorldToEcef(demGeoPoint); // Fastest version!
auto const zero_doppler_info = m_SarSensorModel->ZeroDopplerLookup(xyzCart);
// Only required in 2 cases!
double const range_time = m_SarSensorModel->CalculateRangeTime(xyzCart, zero_doppler_info.sensorPos);
auto const line_sample_yz = m_SarSensorModel->Doppler0ToLineSampleYZ(xyzCart, zero_doppler_info, range_time);
//
// Fill the four channels for output
pixelOutput[0] = line_sample_yz.col_row[0];
pixelOutput[1] = line_sample_yz.col_row[1];
pixelOutput[2] = line_sample_yz.yz[1];
pixelOutput[3] = line_sample_yz.yz[0];
#else
Point2DType col_row(0);
Point2DType y_z(0);
m_SarSensorModel->WorldToLineSampleYZ(demGeoPoint, col_row, y_z);
// Fill the four channels for output
@@ -374,14 +397,14 @@ SARDEMProjectionImageFilter2< TImageIn, TImageOut >
pixelOutput[1] = col_row[1];
pixelOutput[2] = y_z[1];
pixelOutput[3] = y_z[0];
#endif
// Add channels if required
if (m_withXYZ)
{
assert(pixelOutput.Size() >= 6);
#if OTB_VERSION_MAJOR >= 8
xyzCart = Projection::WorldToEcef(demGeoPoint);
#else
#if OTB_VERSION_MAJOR < 8
Point3DType xyzCart;
SarSensorModelAdapter::WorldToCartesian(demGeoPoint, xyzCart);
#endif
pixelOutput[4] = xyzCart[0];
@@ -396,10 +419,18 @@ SARDEMProjectionImageFilter2< TImageIn, TImageOut >
if (m_withSatPos)
{
assert(pixelOutput.Size() >= m_indSatPos + 3);
#if OTB_VERSION_MAJOR >= 8
pixelOutput[m_indSatPos+0] = zero_doppler_info.sensorPos[0];
pixelOutput[m_indSatPos+1] = zero_doppler_info.sensorPos[1];
pixelOutput[m_indSatPos+2] = zero_doppler_info.sensorPos[2];
#else
Point3DType satPos;
Point3DType satVel;
m_SarSensorModel->WorldToSatPositionAndVelocity(demGeoPoint, satPos, satVel);
pixelOutput[m_indSatPos] = satPos[0];
pixelOutput[m_indSatPos+1] = satPos[1];
pixelOutput[m_indSatPos+2] = satPos[2];
#endif
}
}
else
Loading