Add common range and azimuth band filtering running on CPU to Crossmul

cxx/isce3/cuda/signal/gpuCrossMul.cu

@@ -119,8 +119,7 @@ __global__ void flatten_g(thrust::complex<T> *ifg,
 
 
 /** Set number of range looks */
-void
-gpuCrossmul::rangeLooks(int rngLks) {
+void gpuCrossmul::rangeLooks(int rngLks) {
     if (rngLks < 1) {
         std::string error_msg = "ERROR CUDA crossmul range multilook < 1";
         throw isce3::except::InvalidArgument(ISCE_SRCINFO(), error_msg);
@@ -134,8 +133,7 @@ gpuCrossmul::rangeLooks(int rngLks) {
 }
 
 /** Set number of azimuth looks */
-void
-gpuCrossmul::azimuthLooks(int azLks) {
+void gpuCrossmul::azimuthLooks(int azLks) {
     if (azLks < 1) {
         std::string error_msg = "ERROR CUDA crossmul azimuth multilook < 1";
         throw isce3::except::InvalidArgument(ISCE_SRCINFO(), error_msg);
@@ -148,9 +146,8 @@ gpuCrossmul::azimuthLooks(int azLks) {
         _multiLookEnabled = false;
 }
 
-void
-gpuCrossmul::doppler(isce3::core::LUT1d<double> refDoppler,
-                     isce3::core::LUT1d<double> secDoppler)
+void gpuCrossmul::doppler(isce3::core::LUT2d<double> refDoppler,
+                          isce3::core::LUT2d<double> secDoppler)
 {
     _refDoppler = refDoppler;
     _secDoppler = secDoppler;
@@ -214,7 +211,8 @@ gpuCrossmul::crossmul(isce3::io::Raster& refSlcRaster,
                       isce3::io::Raster& secSlcRaster,
                       isce3::io::Raster& ifgRaster,
                       isce3::io::Raster& coherenceRaster,
-                      isce3::io::Raster* rngOffsetRaster) const
+                      isce3::io::Raster* rngOffsetRaster,
+                      isce3::io::Raster* aziOffsetRaster) const
 {
     // set flatten flag based range offset raster ptr value
     bool flatten = rngOffsetRaster ? true : false;

cxx/isce3/cuda/signal/gpuCrossMul.h

@@ -2,9 +2,11 @@
 
 #include "forward.h"
 #include <isce3/io/forward.h> // Raster
+#include <isce3/core/forward.h>
 
 #include <isce3/core/Common.h>
-#include <isce3/core/LUT1d.h>
+#include <isce3/core/LUT2d.h>
+
 #include <thrust/complex.h>
 
 namespace isce3::cuda::signal {
@@ -24,28 +26,32 @@ class gpuCrossmul {
          * \param[out]  coherenceRaster output coherence raster
          * \param[in]  rngOffsetRaster  optional pointer to range offset raster
          *                              if provided, interferogram will be flattened
+         * \param[in]  aziOffsetRaster  optional pointer to azimuth offset raster
+         *                              it is a placeholder here, and will implement the
+         *                              azimuth common band filter based on it
          */
         void crossmul(isce3::io::Raster& refSlcRaster,
                 isce3::io::Raster& secSlcRaster,
                 isce3::io::Raster& ifgRaster,
                 isce3::io::Raster& coherenceRaster,
-                isce3::io::Raster* rngOffsetRaster = nullptr) const;
+                isce3::io::Raster* rngOffsetRaster = nullptr,
+                isce3::io::Raster* aziOffsetRaster = nullptr) const;
 
         /** Set doppler LUTs for reference and secondary SLCs*/
-        void doppler(isce3::core::LUT1d<double> refDoppler,
-                isce3::core::LUT1d<double> secDoppler);
+        void doppler(isce3::core::LUT2d<double> refDoppler,
+                isce3::core::LUT2d<double> secDoppler);
 
         /** Set reference doppler */
-        inline void refDoppler(isce3::core::LUT1d<double> refDopp) {_refDoppler = refDopp;};
+        inline void refDoppler(isce3::core::LUT2d<double> refDopp) {_refDoppler = refDopp;};
 
         /** Get reference doppler */
-        inline const isce3::core::LUT1d<double> & refDoppler() const {return _refDoppler;};
+        inline const isce3::core::LUT2d<double> & refDoppler() const {return _refDoppler;};
 
         /** Set secondary doppler */
-        inline void secDoppler(isce3::core::LUT1d<double> secDopp) {_secDoppler = secDopp;};
+        inline void secDoppler(isce3::core::LUT2d<double> secDopp) {_secDoppler = secDopp;};
 
         /** Get secondary doppler */
-        inline const isce3::core::LUT1d<double> & secDoppler() const {return _secDoppler;};
+        inline const isce3::core::LUT2d<double> & secDoppler() const {return _secDoppler;};
 
         /** Set reference and secondary starting range shift */
         inline void startingRangeShift(double rng_shift) { _offsetStartingRangeShift = rng_shift; }
@@ -93,11 +99,11 @@ class gpuCrossmul {
         inline bool multiLookEnabled() const { return _multiLookEnabled; }
 
     private:
-        //Doppler LUT for the reference SLC
-        isce3::core::LUT1d<double> _refDoppler;
+        //Doppler LUT for the refernce SLC
+        isce3::core::LUT2d<double> _refDoppler;
 
         //Doppler LUT for the secondary SLC
-        isce3::core::LUT1d<double> _secDoppler;
+        isce3::core::LUT2d<double> _secDoppler;
 
         // starting range shifts between the secondary and reference RSLC in meters
         double _offsetStartingRangeShift = 0.0;

cxx/isce3/geocode/GeocodeCov.cpp

@@ -491,7 +491,7 @@ void Geocode<T>::geocodeInterp(
                     std::make_unique<isce3::io::Raster>(
                         vsimem_ref, radar_grid.width(),
                         radar_grid.length(), 1, GDT_Float32, "ENVI");
-                rtc_sigma0_raster = 
+                rtc_sigma0_raster =
                     rtc_raster_sigma0_unique_ptr.get();
             } else {
                 rtc_sigma0_raster =  output_rtc_sigma;
@@ -542,7 +542,7 @@ void Geocode<T>::geocodeInterp(
         block_length = geogrid.length();
     } else {
 
-        if (geocode_memory_mode == 
+        if (geocode_memory_mode ==
                 isce3::core::GeocodeMemoryMode::BlocksGeogridAndRadarGrid) {
             warning << "WARNING the geocode memory mode"
                  << " BlocksGeogridAndRadarGrid is not available"
@@ -554,7 +554,7 @@ void Geocode<T>::geocodeInterp(
         isce3::core::getBlockProcessingParametersY(
             geogrid.length(), geogrid.width(), nbands, sizeof(T),
             &info, &block_length, &nBlocks, min_block_size, max_block_size);
-    } 
+    }
 
     info << "number of blocks: " << nBlocks << pyre::journal::newline;
     info << "block length: " << block_length << pyre::journal::newline;
@@ -1012,7 +1012,6 @@ inline void Geocode<T>::_interpolate(
         }
 
         if (sub_swaths != nullptr) {
-
             // read the sub-swath value at the center to save the output mask
             uint8_t sample_sub_swath_center = sub_swaths->getSampleSubSwath(
                 rdr_y_rslc, rdr_x_rslc);
@@ -1047,10 +1046,10 @@ inline void Geocode<T>::_interpolate(
            }
         }
 
-        /* 
+        /*
         check within the interpolation kernel (approximated by `interp_margin`)
         if any of the samples is marked as shadow or layover-and-shadow
-        in which case we skip to the next position, i.e., we "break" the 
+        in which case we skip to the next position, i.e., we "break" the
         2 inner for-loop bellow (vars: yy and xx) and "continue" from the parent
         for-loop (var: kk) above.
         */
@@ -2079,7 +2078,7 @@ void Geocode<T>::geocodeAreaProj(
                     output_rtc_sigma == nullptr) {
                 std::string vsimem_ref = (
                     "/vsimem/" + getTempString("geocode_cov_areaproj_rtc_sigma0"));
-                rtc_raster_sigma0_unique_ptr = 
+                rtc_raster_sigma0_unique_ptr =
                     std::make_unique<isce3::io::Raster>(
                         vsimem_ref, radar_grid_cropped.width(),
                         radar_grid_cropped.length(), 1, GDT_Float32, "ENVI");
@@ -2259,7 +2258,7 @@ void Geocode<T>::geocodeAreaProj(
     } else {
         isce3::core::getBlockProcessingParametersXY(
                 imax, jmax, nbands + nbands_off_diag_terms, sizeof(T_out),
-                &info, &block_size_with_upsampling_y, &nblocks_y, 
+                &info, &block_size_with_upsampling_y, &nblocks_y,
                 &block_size_with_upsampling_x, &nblocks_x,
                 min_block_size, max_block_size, geogrid_upsampling);
         block_size_x = block_size_with_upsampling_x / geogrid_upsampling;
@@ -2509,7 +2508,7 @@ void Geocode<T>::_runBlock(
         isce3::io::Raster& dem_raster, isce3::io::Raster* out_off_diag_terms,
         isce3::io::Raster* out_geo_rdr, isce3::io::Raster* out_geo_dem,
         isce3::io::Raster* out_geo_nlooks, isce3::io::Raster* out_geo_rtc,
-        isce3::io::Raster* out_geo_rtc_gamma0_to_sigma0, 
+        isce3::io::Raster* out_geo_rtc_gamma0_to_sigma0,
         isce3::core::ProjectionBase* proj, bool flag_apply_rtc,
         bool flag_rtc_raster_is_in_memory, bool flag_rtc_sigma0_raster_is_in_memory,
         isce3::io::Raster* rtc_raster, isce3::io::Raster* rtc_sigma0_raster,
@@ -2905,7 +2904,7 @@ void Geocode<T>::_runBlock(
                 this_block_size_y, this_block_size_x));
 
     nan_t_out *= std::numeric_limits<T_out_real>::quiet_NaN();
- 
+
     for (int band = 0; band < nbands; ++band)
         geoDataBlock[band]->fill(nan_t_out);
 
@@ -3286,7 +3285,7 @@ void Geocode<T>::_runBlock(
                     also need to add `offset_x` and `offset_y` that represent the offsets
                     in X- and Y- directions over the radar-grid coordinates.
 
-                    in which case we skip to the next position, i.e., we "break" the 
+                    in which case we skip to the next position, i.e., we "break" the
                     2 inner for-loop bellow (vars: yy and xx) and "continue" from the parent
                     for-loop (var: kk) above.
                     */
@@ -3647,7 +3646,7 @@ std::string _get_geocode_memory_mode_str(
 }
 
 template<class T>
-void Geocode<T>::_print_parameters(pyre::journal::info_t& channel, 
+void Geocode<T>::_print_parameters(pyre::journal::info_t& channel,
                                   isce3::core::GeocodeMemoryMode& geocode_memory_mode,
                                   const long long min_block_size,
                                   const long long max_block_size) {

cxx/isce3/math/Sinc.h

@@ -1,10 +1,15 @@
 #pragma once
 
 #include <isce3/core/Common.h>
+#include <complex>
 
 namespace isce3 { namespace math {
 
 /** sinc function defined as \f$ \frac{\sin(\pi x)}{\pi x} \f$ */
+template<typename T>
+CUDA_HOSTDEV
+std::complex<T> sinc(std::complex<T> t);
+
 template<typename T>
 CUDA_HOSTDEV
 T sinc(T t);

cxx/isce3/math/Sinc.icc

@@ -90,4 +90,29 @@ T sinc(T t)
     return std::sin(x) / x;
 }
 
+template<typename T>
+CUDA_HOSTDEV
+inline
+std::complex<T> sinc(std::complex<T> t)
+{
+    static_assert(std::is_floating_point<T>::value);
+
+    using U = typename std::remove_cv<T>::type;
+    static constexpr auto eps1 = detail::root_epsilon<U>();
+    static constexpr auto eps2 = detail::fourth_root_epsilon<U>();
+
+    auto u = T(M_PI) * t;
+    auto x = std::abs(u);
+
+    // approximate by Taylor series expansion near zero
+    if (x < eps2) {
+        std::complex<T> out = static_cast<T>(1);
+        if (x > eps1) {
+            out -= u * u / static_cast<T>(6);
+        }
+        return out;
+    }
+    return std::sin(u) / u;
+}
+
 }}

cxx/isce3/product/RadarGridProduct.h

@@ -27,7 +27,7 @@ namespace isce3 { namespace product {
 /**
  * Return the path to each child group of `group` that ends with the substring
  * `group_name`.
- * 
+ *
  * \param[in] group       Parent group
  * \param[in] group_name  Search string
  * \returns               List of child group paths
@@ -39,13 +39,13 @@ std::vector<std::string> findGroupPath(
  * Return grids or swaths group paths within the base_group.
  * Start by assigning an empty string to image_group_str in case
  * grids and swaths group are not found.
- * 
+ *
  * \param[in] file
  * \param[in] base_dir           Path to `base_group` object (e.g. '/science/')
  * \param[in] base_group         Base group
  * \param[in] key_vector         Vector containing possible image groups
  * (e.g., 'swaths', 'grids', or both) to look for
- * \param[out] image_group_str   Path to first image group found containing 
+ * \param[out] image_group_str   Path to first image group found containing
  * one of the `key_vector` keys (e.g., '/science/LSAR/RSLC/swaths')
  * \param[in] metadata_group_str Path to first metadata group found by
  * substituting `key` with `metadata` in `image_group_str`
@@ -60,10 +60,10 @@ void setImageMetadataGroupStr(
         std::string &metadata_group_str);
 
 /** RadarGridProduct class declaration
- * 
+ *
  * The Produt attribute Swaths map, i.e. _swaths, associates the
- * frequency (key) with the Swath object (value). The RadarGridProduct object 
- * is usually initiated with an empty map and the serialization of 
+ * frequency (key) with the Swath object (value). The RadarGridProduct object
+ * is usually initiated with an empty map and the serialization of
  * the SAR product is responsible for populating the Swath map
  * from the product's metadata.
  *
@@ -75,7 +75,7 @@ class RadarGridProduct {
         RadarGridProduct(isce3::io::IH5File &);
         /** Constructor with Metadata and Swath map. */
         inline RadarGridProduct(const Metadata &, const std::map<char, isce3::product::Swath> &);
-
+    
         /** Get a read-only reference to the metadata */
         inline const Metadata & metadata() const { return _metadata; }
         /** Get a reference to the metadata. */

cxx/isce3/product/Swath.h

@@ -118,7 +118,7 @@ class isce3::product::Swath {
         inline isce3::product::SubSwaths& subSwaths() {
             return _subSwaths;
         }
-        
+
         inline const isce3::product::SubSwaths& subSwaths() const {
             return _subSwaths;
         }