Ant geom (#726)

* Add a new module "antenna::geometryfunc" and update  "geometry::geometry" module

- Add pybind11 module for "geometryfunc"

- Add a python test suite for "geometryfunc"

- Modify "srPosFromLookVecDem" in "geometry" with
  new initial value

- Fix typo in doc of pybind "ltpcoordinates.cpp" of
  "geometry"

- Replace all instances of scalar input argument DEM height
  "dem_hgt" with "isce3::geometry::DEMInterpolator" object

- Update all respective C++ and pybind11 scripts as well as their
  test suites.

- Update doxygen and numpy docstrings accordingly

- Run clang-format and autopep8

Author: rad-eng-59

cxx/isce3/Headers.cmake

@@ -1,5 +1,6 @@
 set(HEADERS
 antenna/forward.h
+antenna/geometryfunc.h
 antenna/Frame.h
 antenna/Frame.icc
 antenna/SphGridType.h

cxx/isce3/Sources.cmake

@@ -1,4 +1,5 @@
 set(SRCS
+antenna/geometryfunc.cpp
 core/Attitude.cpp
 core/Baseline.cpp
 core/Basis.cpp

cxx/isce3/antenna/geometryfunc.cpp

@@ -0,0 +1,158 @@
+// Definition of the antenna-related geometry functions
+#include "geometryfunc.h"
+
+#include <isce3/core/Quaternion.h>
+#include <isce3/core/Vector.h>
+#include <isce3/except/Error.h>
+#include <isce3/geometry/geometry.h>
+
+// Aliases
+using namespace isce3::core;
+namespace geom = isce3::geometry;
+namespace ant = isce3::antenna;
+using VecXd = Eigen::VectorXd;
+
+// Helper local functions
+
+/**
+ * @internal
+ * Helper function to get slant range and LLH and convergence flag
+ * @param[in] el_theta : either elevation or theta angle in radians
+ * depending on the "frame" object.
+ * @param[in] az_phi : either azimuth or phi angle in radians depending
+ * on the "frame" object.
+ * @param[in] pos_ecef : antenna/spacecraft position in ECEF (m,m,m)
+ * @param[in] vel_ecef : spacecraft velocity scaled by -2/wavalength in
+ * ECEF (m/s,m/s,m/s)
+ * @param[in] quat : isce3 quaternion object for transformation from antenna
+ * body-fixed to ECEF
+ * @param[in] dem_interp (optional): isce3 DEMInterpolator object
+ * w.r.t ellipsoid. Default is zero height.
+ * @param[in] abs_tol (optional): Abs error/tolerance in height estimation (m)
+ * between desired input height and final output height. Default is 0.5.
+ * @param[in] max_iter (optional): Max number of iterations in height
+ * estimation. Default is 10.
+ * @param[in] frame (optional): isce3 Frame object to define antenna spherical
+ *  coordinate system. Default is based on "EL_AND_AZ" spherical grid.
+ * @param[in] ellips (optional): isce3 Ellipsoid object defining the
+ * ellipsoidal planet. Default is WGS84 ellipsoid.
+ * @return a tuple of three values : slantrange (m), Doppler (Hz), a bool
+ * which is true if height tolerance is met, false otherwise.
+ * @exception InvalidArgument, RuntimeError
+ */
+static std::tuple<double, double, bool> _get_sr_dop_conv(double el_theta,
+        double az_phi, const Vec3& pos_ecef, const Vec3& vel_ecef,
+        const Quaternion& quat, const geom::DEMInterpolator& dem_interp,
+        double abs_tol, int max_iter, const ant::Frame& frame,
+        const Ellipsoid& ellips)
+{
+    // pointing sphercial to cartesian in Antenna frame
+    auto pnt_xyz = frame.sphToCart(el_theta, az_phi);
+    // pointing from Ant XYZ to global ECEF
+    auto pnt_ecef = quat.rotate(pnt_xyz);
+    // get slant range and target position on/above the ellipsoid
+    Vec3 tg_ecef, tg_llh;
+    double sr;
+    auto iter_info = geom::srPosFromLookVecDem(sr, tg_ecef, tg_llh, pos_ecef,
+            pnt_ecef, dem_interp, abs_tol, max_iter, ellips);
+    bool convergence {true};
+    if (iter_info.second > abs_tol)
+        convergence = false;
+
+    double doppler = vel_ecef.dot(pnt_ecef);
+
+    return {sr, doppler, convergence};
+}
+
+// Antenna to Radar functions
+std::tuple<double, double, bool> ant::ant2rgdop(double el_theta, double az_phi,
+        const Vec3& pos_ecef, const Vec3& vel_ecef, const Quaternion& quat,
+        double wavelength, const geom::DEMInterpolator& dem_interp,
+        double abs_tol, int max_iter, const ant::Frame& frame,
+        const Ellipsoid& ellips)
+{
+    if (!(wavelength > 0.0))
+        throw isce3::except::InvalidArgument(
+                ISCE_SRCINFO(), "Bad value for wavelength!");
+
+    const auto vel_ecef_cst = (-2. / wavelength) * vel_ecef;
+
+    return _get_sr_dop_conv(el_theta, az_phi, pos_ecef, vel_ecef_cst, quat,
+            dem_interp, abs_tol, max_iter, frame, ellips);
+}
+
+std::tuple<VecXd, VecXd, bool> ant::ant2rgdop(
+        const Eigen::Ref<const VecXd>& el_theta, double az_phi,
+        const Vec3& pos_ecef, const Vec3& vel_ecef, const Quaternion& quat,
+        double wavelength, const geom::DEMInterpolator& dem_interp,
+        double abs_tol, int max_iter, const ant::Frame& frame,
+        const Ellipsoid& ellips)
+{
+    if (wavelength <= 0.0)
+        throw isce3::except::InvalidArgument(
+                ISCE_SRCINFO(), "Bad value for wavelength!");
+
+    // initialization and vector allocations
+    const auto vel_ecef_cst = (-2. / wavelength) * vel_ecef;
+    auto ang_size = el_theta.size();
+    VecXd slantrange(ang_size);
+    VecXd doppler(ang_size);
+    bool converge {true};
+
+#pragma omp parallel for
+    for (decltype(ang_size) idx = 0; idx < ang_size; ++idx) {
+        auto [sr, dop, flag] =
+                _get_sr_dop_conv(el_theta[idx], az_phi, pos_ecef, vel_ecef_cst,
+                        quat, dem_interp, abs_tol, max_iter, frame, ellips);
+        slantrange(idx) = sr;
+        doppler(idx) = dop;
+        if (!flag)
+            converge = false;
+    }
+    return {slantrange, doppler, converge};
+}
+
+// Antenna to Geometry
+std::tuple<Vec3, bool> ant::ant2geo(double el_theta, double az_phi,
+        const Vec3& pos_ecef, const Quaternion& quat,
+        const geom::DEMInterpolator& dem_interp, double abs_tol, int max_iter,
+        const ant::Frame& frame, const Ellipsoid& ellips)
+{
+    // pointing sphercial to cartesian in Antenna frame
+    auto pnt_xyz = frame.sphToCart(el_theta, az_phi);
+    // pointing from Ant XYZ to global ECEF
+    auto pnt_ecef = quat.rotate(pnt_xyz);
+    // get slant range and target position on/above the ellipsoid
+    Vec3 tg_ecef, tg_llh;
+    double sr;
+    auto iter_info = geom::srPosFromLookVecDem(sr, tg_ecef, tg_llh, pos_ecef,
+            pnt_ecef, dem_interp, abs_tol, max_iter, ellips);
+    bool convergence {true};
+    if (iter_info.second > abs_tol)
+        convergence = false;
+
+    return {tg_llh, convergence};
+}
+
+std::tuple<std::vector<Vec3>, bool> ant::ant2geo(
+        const Eigen::Ref<const VecXd>& el_theta, double az_phi,
+        const Vec3& pos_ecef, const Quaternion& quat,
+        const geom::DEMInterpolator& dem_interp, double abs_tol, int max_iter,
+        const ant::Frame& frame, const Ellipsoid& ellips)
+{
+    // initialize and allocate vectors
+    auto ang_size {el_theta.size()};
+    std::vector<Vec3> tg_llh_vec(ang_size);
+    bool converge {true};
+
+#pragma omp parallel for
+    for (decltype(ang_size) idx = 0; idx < ang_size; ++idx) {
+        auto [tg_llh, flag] = ant::ant2geo(el_theta(idx), az_phi, pos_ecef,
+                quat, dem_interp, abs_tol, max_iter, frame, ellips);
+
+        if (!flag)
+            converge = false;
+        tg_llh_vec[idx] = tg_llh;
+    }
+    return {tg_llh_vec, converge};
+}

cxx/isce3/antenna/geometryfunc.h

@@ -0,0 +1,171 @@
+/** @file geometryfunc.h
+ * A collection of antenna-related geometry functions.
+ */
+#pragma once
+
+#include <isce3/core/forward.h>
+
+#include <tuple>
+#include <vector>
+
+#include <Eigen/Dense>
+
+#include <isce3/antenna/Frame.h>
+#include <isce3/core/Ellipsoid.h>
+#include <isce3/geometry/DEMInterpolator.h>
+
+/** @namespace isce3::antenna */
+namespace isce3 { namespace antenna {
+
+// Antenna to Radar
+
+/**
+ * Estimate Radar products, Slant range and Doppler centroid, from
+ * spherical angles in antenna body-fixed domain for a certain spacecraft
+ * position, velocity,and attitude at a certain height w.r.t. an ellipsoid.
+ * @param[in] el_theta : either elevation or theta angle in radians
+ * depending on the "frame" object.
+ * @param[in] az_phi : either azimuth or phi angle in radians depending
+ * on the "frame" object.
+ * @param[in] pos_ecef : antenna/spacecraft position in ECEF (m,m,m)
+ * @param[in] vel_ecef : spacecraft velocity in ECEF (m/s,m/s,m/s)
+ * @param[in] quat : isce3 quaternion object for transformation from antenna
+ * body-fixed to ECEF
+ * @param[in] wavelength : Radar wavelength in (m).
+ * @param[in] dem_interp (optional): isce3 DEMInterpolator object
+ * w.r.t ellipsoid. Default is zero height.
+ * @param[in] abs_tol (optional): Abs error/tolerance in height estimation (m)
+ * between desired input height and final output height. Default is 0.5.
+ * @param[in] max_iter (optional): Max number of iterations in height
+ * estimation. Default is 10.
+ * @param[in] frame (optional): isce3 Frame object to define antenna spherical
+ *  coordinate system. Default is based on "EL_AND_AZ" spherical grid.
+ * @param[in] ellips (optional): isce3 Ellipsoid object defining the
+ * ellipsoidal planet. Default is WGS84 ellipsoid.
+ * @return slantrange (m)
+ * @return Doppler (Hz)
+ * @return a bool which is true if height tolerance is met, false otherwise.
+ * @exception InvalidArgument, RuntimeError
+ * @cite ReeTechDesDoc
+ */
+
+std::tuple<double, double, bool> ant2rgdop(double el_theta, double az_phi,
+        const isce3::core::Vec3& pos_ecef, const isce3::core::Vec3& vel_ecef,
+        const isce3::core::Quaternion& quat, double wavelength,
+        const isce3::geometry::DEMInterpolator& dem_interp = {},
+        double abs_tol = 0.5, int max_iter = 10,
+        const isce3::antenna::Frame& frame = {},
+        const isce3::core::Ellipsoid& ellips = {});
+
+/**
+ * Overloaded function to estimate Radar products, Slant ranges and Doppler
+ * centroids, from spherical angles in antenna body-fixed domain for a
+ * certain spacecraft position, velocity,and attitude at a certain height
+ * w.r.t. an ellipsoid.
+ * @param[in] el_theta : a vector of either elevation or theta angle in
+ * radians depending on the "frame" object.
+ * @param[in] az_phi : either azimuth or phi angle in radians depending
+ * on the "frame" object.
+ * @param[in] pos_ecef : antenna/spacecraft position in ECEF (m,m,m)
+ * @param[in] vel_ecef : spacecraft velocity in ECEF (m/s,m/s,m/s)
+ * @param[in] quat : isce3 quaternion object for transformation from antenna
+ * body-fixed to ECEF
+ * @param[in] wavelength : Radar wavelength in (m).
+ * @param[in] dem_interp (optional): isce3 DEMInterpolator object
+ * w.r.t ellipsoid. Default is zero height.
+ * @param[in] abs_tol (optional): Abs error/tolerance in height estimation (m)
+ * between desired input height and final output height. Default is 0.5.
+ * @param[in] max_iter (optional): Max number of iterations in height
+ * estimation. Default is 10.
+ * @param[in] frame (optional): isce3 Frame object to define antenna spherical
+ *  coordinate system. Default is based on "EL_AND_AZ" spherical grid.
+ * @param[in] ellips (optional): isce3 Ellipsoid object defining the
+ * ellipsoidal planet. Default is WGS84 ellipsoid.
+ * @return an Eigen::VectorXd of slant ranges (m)
+ * @return an Eigen::VectorXd of Doppler values (Hz)
+ * @return a bool which is true if height tolerance is met, false otherwise.
+ * @exception InvalidArgument, RuntimeError
+ * @cite ReeTechDesDoc
+ */
+std::tuple<Eigen::VectorXd, Eigen::VectorXd, bool> ant2rgdop(
+        const Eigen::Ref<const Eigen::VectorXd>& el_theta, double az_phi,
+        const isce3::core::Vec3& pos_ecef, const isce3::core::Vec3& vel_ecef,
+        const isce3::core::Quaternion& quat, double wavelength,
+        const isce3::geometry::DEMInterpolator& dem_interp = {},
+        double abs_tol = 0.5, int max_iter = 10,
+        const isce3::antenna::Frame& frame = {},
+        const isce3::core::Ellipsoid& ellips = {});
+
+// Antenna to Geometry
+
+/**
+ * Estimate geodetic geolocation (longitude, latitude, height) from
+ * spherical angles in antenna body-fixed domain for a certain spacecraft
+ * position and attitude at a certain height w.r.t. an ellipsoid.
+ * @param[in] el_theta : either elevation or theta angle in radians
+ * depending on the "frame" object.
+ * @param[in] az_phi : either azimuth or phi angle in radians depending
+ * on the "frame" object.
+ * @param[in] pos_ecef : antenna/spacecraft position in ECEF (m,m,m)
+ * @param[in] quat : isce3 quaternion object for transformation from antenna
+ * body-fixed to ECEF
+ * @param[in] dem_interp (optional): isce3 DEMInterpolator object
+ * w.r.t ellipsoid. Default is zero height.
+ * @param[in] abs_tol (optional): Abs error/tolerance in height estimation (m)
+ * between desired input height and final output height. Default is 0.5.
+ * @param[in] max_iter (optional): Max number of iterations in height
+ * estimation. Default is 10.
+ * @param[in] frame (optional): isce3 Frame object to define antenna spherical
+ *  coordinate system. Default is based on "EL_AND_AZ" spherical grid.
+ * @param[in] ellips (optional): isce3 Ellipsoid object defining the
+ * ellipsoidal planet. Default is WGS84 ellipsoid.
+ * @return an isce3::core::Vec3 of geodetic LLH(longitude,latitude,height)
+ * in (rad,rad,m)
+ * @return a bool which is true if height tolerance is met, false otherwise.
+ * @exception InvalidArgument, RuntimeError
+ * @cite ReeTechDesDoc
+ */
+std::tuple<isce3::core::Vec3, bool> ant2geo(double el_theta, double az_phi,
+        const isce3::core::Vec3& pos_ecef, const isce3::core::Quaternion& quat,
+        const isce3::geometry::DEMInterpolator& dem_interp = {},
+        double abs_tol = 0.5, int max_iter = 10,
+        const isce3::antenna::Frame& frame = {},
+        const isce3::core::Ellipsoid& ellips = {});
+
+/**
+ * Overloaded function to estimate geodetic geolocation
+ * (longitude, latitude, height) from spherical angles in antenna
+ * body-fixed domain for a certain spacecraft position and attitude
+ * at a certain height w.r.t. an ellipsoid.
+ * @param[in] el_theta : a vector of either elevation or theta angle
+ * in radians depending on the "frame" object.
+ * @param[in] az_phi : either azimuth or phi angle in radians depending
+ * on the "frame" object.
+ * @param[in] pos_ecef : antenna/spacecraft position in ECEF (m,m,m)
+ * @param[in] quat : isce3 quaternion object for transformation from antenna
+ * body-fixed to ECEF
+ * @param[in] dem_interp (optional): isce3 DEMInterpolator object
+ * w.r.t ellipsoid. Default is zero height.
+ * @param[in] abs_tol (optional): Abs error/tolerance in height estimation (m)
+ * between desired input height and final output height. Default is 0.5.
+ * @param[in] max_iter (optional): Max number of iterations in height
+ * estimation. Default is 10.
+ * @param[in] frame (optional): isce3 Frame object to define antenna spherical
+ *  coordinate system. Default is based on "EL_AND_AZ" spherical grid.
+ * @param[in] ellips (optional): isce3 Ellipsoid object defining the
+ * ellipsoidal planet. Default is WGS84 ellipsoid.
+ * @return a vector of isce3::core::Vec3 of geodetic
+ * LLH(longitude,latitude,height) in (rad,rad,m)
+ * @return a bool which is true if height tolerance is met, false otherwise.
+ * @exception InvalidArgument, RuntimeError
+ * @cite ReeTechDesDoc
+ */
+std::tuple<std::vector<isce3::core::Vec3>, bool> ant2geo(
+        const Eigen::Ref<const Eigen::VectorXd>& el_theta, double az_phi,
+        const isce3::core::Vec3& pos_ecef, const isce3::core::Quaternion& quat,
+        const isce3::geometry::DEMInterpolator& dem_interp = {},
+        double abs_tol = 0.5, int max_iter = 10,
+        const isce3::antenna::Frame& frame = {},
+        const isce3::core::Ellipsoid& ellips = {});
+
+}} // namespace isce3::antenna

cxx/isce3/geometry/geometry.cpp

@@ -452,19 +452,28 @@ double isce3::geometry::slantRangeFromLookVec(
 
 std::pair<int, double> isce3::geometry::srPosFromLookVecDem(double& sr,
         Vec3& tg_pos, Vec3& llh, const Vec3& sc_pos, const Vec3& lkvec,
-        double dem_hgt, double hgt_err, int num_iter, const Ellipsoid& ellips)
+        const DEMInterpolator& dem_interp, double hgt_err, int num_iter,
+        const Ellipsoid& ellips)
 {
     if (hgt_err <= 0.0 || num_iter <= 0)
         throw isce3::except::InvalidArgument(ISCE_SRCINFO(),
                 "Height Error and number of iteration "
                 "must be non-zero positive values!");
-    sr = slantRangeFromLookVec(sc_pos, lkvec, ellips);
+
+    // form a new ellipsoid whose radii adjusted by mean height
+    const auto mean_hgt = dem_interp.meanHeight();
+    const auto a_new = ellips.a() + mean_hgt;
+    const auto b_new = ellips.b() + mean_hgt;
+    const auto e2_new = 1.0 - (b_new * b_new) / (a_new * a_new);
+    // initial guess of slant range per new ellipsoid
+    sr = slantRangeFromLookVec(sc_pos, lkvec, Ellipsoid(a_new, e2_new));
     int cnt {0};
     double abs_hgt_dif;
     do {
         tg_pos = sr * lkvec + sc_pos;
         llh = ellips.xyzToLonLat(tg_pos);
-        auto hgt_dif {dem_hgt - llh(2)};
+        auto dem_hgt = dem_interp.interpolateLonLat(llh(0), llh(1));
+        auto hgt_dif = dem_hgt - llh(2);
         sr += hgt_dif / _downVal(llh(0), llh(1), tg_pos);
         abs_hgt_dif = std::abs(hgt_dif);
         ++cnt;