Resolve "Factorize Computations in SARDEMProjection"

Closes #7 (closed)

include/otbSARDEMProjectionImageFilter2.hxx

 /*
- * 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