Compute statistics (#811)

* add Stats struct and computeRasterStats() function at C++ level

* update cmake files with new Stats files

* add pybind11 bindings to Stats struct and computeRasterStats() function

* add stats computation to the GCOV workflow

* add unitest

* update complex operations

* update RTC block mode enums to follow recently added Stats block mode enums

* fix CUDA complex_operations division operator

* update enum corresponding to block processing in the Y-direction

* remove prefix "k" from enum class memoryModeBlockY()

* rename pybind_nisar to nisar

* rename pybind_nisar to nisar (2)

* merge branch develop into add_stats (2)

* merge branch develop into add_stats (3)

* rename pybind_nisar to nisar (3)

* add computation of complex-valued statistics considering real and imaginary parts separately rather than using polar representation

* add stats to GSLC products

* add stats to GSLC products (2)

* add stats to GUNW products

* add stats to RSLC products

* update stats for GCOV products

* fix standard deviation formula for sample standard deviation

* update unitest thresholds

* add Python unitest to compare results with NumPy

* do not update Raster's metadata statistics

* make BlockStats struct private (remove it from the header file)

* remove unnecessary include

* uncomment complex Stats pybind declaration

* fix python docstrings identation

* substitute GDAL GetStatistics() with ComputeStatistics() and make error thresholds more strict

* remove complexMultiply and make complexOperations lowercase

* remove complexMultiply (2)

* fix identation in BlockStats

* improve Stats docstrings

* enable parallelization

* move declaration of _compute_stats() before its use

* remove file merged by mistake

* run stats unitest also for GeocodeCov interp results

* move _compute_stats() back to the end of the file

* include pyre/journal.h

* use templates in CUDA version of ComplexOperations

* update docstrings for Stats and StatsRealImag structs

* remove remaining ComplexMultiply.h

* add docstrings to SLC writer compute_stats()

* simplify assignment of block length

* separate declaration of real- and complex-valued Stats struct

* remove unnecessary if statement from Stats unitest in geocode.cpp

* fix CUDA version of ComplexOperations

* fix the namespace of CUDA complexOperations, remove operation* from GeocodeHelpers

* remove unused variable output_raster_dataset

* remove complexOperations from Interp1d.icc

* add omp critical to getBlock() to avoid concurrence issues, improve info messages

* make compute_raster_stats() template instantiation explicit by adding a suffix to the function name

* update python bindings with new function names

* lowercase ComplexOperation.h files

* remove layoverShadowMask stats from GUNW products

* reduce the number of function calls to signedRealOrComplexModulus()

* Update cxx/isce3/math/Stats.h

Co-authored-by: Brian P Hawkins <[email protected]>

* remove schedule(dynamic) from omp directive

* Update cxx/isce3/math/Stats.cpp

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Update cxx/isce3/math/Stats.cpp

Co-authored-by: Geoffrey M Gunter <[email protected]>

* Update cxx/isce3/math/Stats.cpp

Co-authored-by: Geoffrey M Gunter <[email protected]>

* limit scope of "using operator*" to functions require it

Co-authored-by: Gustavo H Shiroma <[email protected]>
Co-authored-by: Brian P Hawkins <[email protected]>
Co-authored-by: Geoffrey M Gunter <[email protected]>

Author: 4 people

cxx/isce3/Headers.cmake

@@ -127,7 +127,8 @@ io/Raster.h
 io/Raster.icc
 io/Serialization.h
 math/Bessel.h
-math/ComplexMultiply.h
+math/complexOperations.h
+math/Stats.h
 math/detail/RootFind1dBase.h
 math/polyfunc.h
 math/RootFind1dNewton.h

cxx/isce3/Sources.cmake

@@ -63,6 +63,7 @@ io/IH5.cpp
 io/IH5Dataset.cpp
 io/Raster.cpp
 math/Bessel.cpp
+math/Stats.cpp
 math/polyfunc.cpp
 math/RootFind1dNewton.cpp
 math/RootFind1dSecant.cpp

cxx/isce3/core/Interp1d.icc

@@ -1,12 +1,12 @@
-#include <isce3/math/ComplexMultiply.h>
+#include <isce3/math/complexOperations.h>
 
 namespace isce3 { namespace core {
 
 template<typename TK, typename TD>
 TD interp1d(const Kernel<TK>& kernel, const TD* x, size_t length, size_t stride,
             double t, bool periodic)
 {
-    using namespace isce3::math::complex_multiply;
+    using namespace isce3::math::complex_operations;
     int _width = int(ceil(kernel.width()));
     long i0 = 0;
     if (_width % 2 == 0) {

cxx/isce3/cuda/Headers.cmake

@@ -44,7 +44,7 @@ image/gpuResampSlc.h
 image/ResampSlc.h
 io/DataStream.h
 io/DataStream.icc
-math/ComplexMultiply.h
+math/complexOperations.h
 signal/forward.h
 signal/gpuCrossMul.h
 signal/gpuFilter.h

cxx/isce3/cuda/core/Interp1d.icc

@@ -1,15 +1,12 @@
 #include <cmath>
 #include <type_traits>
 
-#include <isce3/cuda/math/ComplexMultiply.h>
-
 namespace isce3 { namespace cuda { namespace core {
 
 template<class Kernel, typename T>
 CUDA_HOSTDEV inline T interp1d(const Kernel& kernel, const T* x, size_t length,
                                size_t stride, double t, bool periodic)
 {
-    using isce3::math::complex_multiply::operator*;
 
     auto width = static_cast<int>(std::ceil(kernel.width()));
     auto i0 = (width % 2 == 0) ? static_cast<long>(std::ceil(t))

cxx/isce3/cuda/math/ComplexMultiply.h

@@ -0,0 +1,39 @@
+#pragma once
+
+#include <thrust/complex.h>
+
+#include <isce3/core/Common.h>
+#include <isce3/math/complexOperations.h>
+
+/** Promotion rules so T * complex<U> works for {T,U} in {float,double}. */
+namespace isce3 { namespace cuda { namespace math { namespace complex_operations {
+
+template<typename T1, typename T2>
+CUDA_HOSTDEV inline auto operator*(T1 t, thrust::complex<T2> u)
+{
+    using U = typename std::common_type_t<T1, T2>;
+    return U(t) * thrust::complex<U>(u);
+}
+
+template<typename T1, typename T2>
+CUDA_HOSTDEV inline auto operator*(thrust::complex<T1> t, T2 u)
+{
+    using U = typename std::common_type_t<T1, T2>;
+    return thrust::complex<U>(t) * U(u);
+}
+
+template<typename T1, typename T2>
+CUDA_HOSTDEV inline auto operator/(T1 t, thrust::complex<T2> u)
+{
+    using U = typename std::common_type_t<T1, T2>;
+    return U(t) / thrust::complex<U>(u);
+}
+
+template<typename T1, typename T2>
+CUDA_HOSTDEV inline auto operator/(thrust::complex<T1> t, T2 u)
+{
+    using U = typename std::common_type_t<T1, T2>;
+    return thrust::complex<U>(t) / U(u);
+}
+
+}}}} // namespace isce3::cuda::math::complex_operations

cxx/isce3/cuda/math/complexOperations.h

@@ -374,8 +374,8 @@ void Geocode<T>::geocodeInterp(
             if (std::isnan(rtc_geogrid_upsampling))
                 rtc_geogrid_upsampling = 1;
 
-            isce3::geometry::rtcMemoryMode rtc_memory_mode =
-                    isce3::geometry::RTC_AUTO;
+            isce3::core::MemoryModeBlockY rtc_memory_mode = 
+                isce3::core::MemoryModeBlockY::AutoBlocksY;
             int radar_grid_nlooks = 1;
             computeRtc(demRaster, *rtc_raster, radar_grid, _orbit, _doppler,
                     _geoGridStartY, _geoGridSpacingY, _geoGridStartX,
@@ -729,6 +729,8 @@ inline void Geocode<T>::_interpolate(
         isce3::core::Matrix<float>& out_geo_rtc_array)
 {
 
+    using isce3::math::complex_operations::operator*;
+
     size_t length = geoDataBlock.length();
     size_t width = geoDataBlock.width();
     // Add extra margin for interpolation
@@ -1677,13 +1679,13 @@ void Geocode<T>::geocodeAreaProj(
             else if (std::isnan(rtc_geogrid_upsampling))
                 rtc_geogrid_upsampling = 2 * geogrid_upsampling;
 
-            isce3::geometry::rtcMemoryMode rtc_memory_mode;
+            isce3::core::MemoryModeBlockY rtc_memory_mode;
             if (geocode_memory_mode == geocodeMemoryMode::AUTO)
-                rtc_memory_mode = isce3::geometry::RTC_AUTO;
+                rtc_memory_mode = isce3::core::MemoryModeBlockY::AutoBlocksY;
             else if (geocode_memory_mode == geocodeMemoryMode::SINGLE_BLOCK)
-                rtc_memory_mode = isce3::geometry::RTC_SINGLE_BLOCK;
+                rtc_memory_mode = isce3::core::MemoryModeBlockY::SingleBlockY;
             else
-                rtc_memory_mode = isce3::geometry::RTC_BLOCKS_GEOGRID;
+                rtc_memory_mode = isce3::core::MemoryModeBlockY::MultipleBlocksY;
 
             computeRtc(dem_raster, *rtc_raster, radar_grid_cropped, _orbit,
                     _doppler, _geoGridStartY, _geoGridSpacingY, _geoGridStartX,
@@ -2015,6 +2017,8 @@ void Geocode<T>::_runBlock(
         pyre::journal::info_t& info)
 {
 
+    using isce3::math::complex_operations::operator*;
+
     // start (az) and r0 at the outer edge of the first pixel
     const double pixazm = radar_grid.azimuthTimeInterval();
     const double start = radar_grid.sensingStart() - 0.5 * pixazm;

cxx/isce3/geocode/GeocodeCov.cpp

@@ -1,21 +1,8 @@
 #pragma once
 
-namespace isce3 { namespace geocode {
-
-
-template<typename T1, typename T2>
-auto operator*(const std::complex<T1>& lhs, const T2& rhs)
-{
-    using U = typename std::common_type_t<T1, T2>;
-    return std::complex<U>(lhs) * U(rhs);
-}
+#include <isce3/math/complexOperations.h>
 
-template<typename T1, typename T2>
-auto operator*(const T1& lhs, const std::complex<T2>& rhs)
-{
-    using U = typename std::common_type_t<T1, T2>;
-    return U(lhs) * std::complex<U>(rhs);
-}
+namespace isce3 { namespace geocode {
 
 template<typename T, typename T_out>
 void _convertToOutputType(T a, T_out& b)
@@ -40,6 +27,9 @@ void _accumulate(T_out& band_value, T a, double b)
 {
     if (b == 0)
         return;
+
+    using isce3::math::complex_operations::operator*;
+
     T_out a2;
     _convertToOutputType(a, a2);
     band_value += a2 * b;

cxx/isce3/geocode/GeocodeHelpers.h

@@ -212,8 +212,8 @@ void GeocodePolygon<T>::getPolygonMean(
         isce3::geometry::rtcAlgorithm rtc_algorithm =
                 isce3::geometry::rtcAlgorithm::RTC_AREA_PROJECTION;
 
-        isce3::geometry::rtcMemoryMode rtc_memory_mode =
-                isce3::geometry::rtcMemoryMode::RTC_SINGLE_BLOCK;
+        isce3::core::MemoryModeBlockY rtc_memory_mode =
+                isce3::core::MemoryModeBlockY::SingleBlockY;
 
         computeRtc(radar_grid_cropped, _orbit, input_dop, dem_raster,
                    *rtc_raster, input_terrain_radiometry,
@@ -288,6 +288,8 @@ void GeocodePolygon<T>::_getPolygonMean(
 
     pyre::journal::info_t _info("isce.geometry._getPolygonMean");
 
+    using isce3::math::complex_operations::operator*;
+
     // number of bands in the input raster
     const int nbands = input_raster.numBands();
     _info << "nbands: " << nbands << pyre::journal::endl;