Pybind interface for filter2D (#714)

* Add bindings for filter2D

* add filter2D bindings to Sources.make

* Add filter2D bindings in signal and correct typos

* Introduce default number of block_rows

* Add filter2D unit test

* Correcting typo in pybind test

* Change pybind implementation to lambda function

* Modify unit test for complex data

* Reformat pybind code to abide pep8

* Format unit test to abide pep8 notation

* Add pybind args and descriptive string

* Introduce kernel and input/output dimension checks

* Add IndexError as exception in filter2D unit test

* Specifying band number

* Insert RunTimeError exception as default option in switch block

* Switching the include filter2D statement in filter2D.cpp

* Simplify unit test using scipy.signal.convolve2D, assertion statements, add magnitude difference check for complex data

* Adding const to kernel_rows kernel_columns at C++ level

* Adding const to kernel variable in pybind11 interface

* Removing unused variables and inputs

* Correcting comment on data type

* Address input/output different size in pybind11 interface

Co-authored-by: vbrancat <[email protected]>

Author: 2 people

cxx/isce3/signal/filter2D.cpp

@@ -9,8 +9,8 @@
 #include <isce3/signal/convolve.h>
 #include <isce3/signal/decimate.h>
 
-void check_kernels(std::valarray<double>& kernel_columns,
-                   std::valarray<double>& kernel_rows)
+void check_kernels(const std::valarray<double>& kernel_columns,
+                   const std::valarray<double>& kernel_rows)
 {
 
     size_t ncols_kernel = kernel_columns.size();
@@ -128,8 +128,8 @@ void setup_block_parameters(const int nrows, const int blockRows,
 template<typename T>
 void isce3::signal::filter2D(isce3::io::Raster& output_raster,
                              isce3::io::Raster& input_raster,
-                             std::valarray<double>& kernel_columns,
-                             std::valarray<double>& kernel_rows, int block_rows)
+                             const std::valarray<double>& kernel_columns,
+                             const std::valarray<double>& kernel_rows, int block_rows)
 {
 
     bool do_decimate = false;
@@ -164,8 +164,8 @@ template<typename T>
 void isce3::signal::filter2D(isce3::io::Raster& output_raster,
                              isce3::io::Raster& input_raster,
                              isce3::io::Raster& mask_raster,
-                             std::valarray<double>& kernel_columns,
-                             std::valarray<double>& kernel_rows, int block_rows)
+                             const std::valarray<double>& kernel_columns,
+                             const std::valarray<double>& kernel_rows, int block_rows)
 {
 
     std::cout << "A mask is provided. The input will be masked before filtering"
@@ -188,8 +188,8 @@ template<typename T>
 void isce3::signal::filter2D(isce3::io::Raster& output_raster,
                              isce3::io::Raster& input_raster,
                              isce3::io::Raster& mask_raster,
-                             std::valarray<double>& kernel_columns,
-                             std::valarray<double>& kernel_rows,
+                             const std::valarray<double>& kernel_columns,
+                             const std::valarray<double>& kernel_rows,
                              const bool do_decimate, const bool mask_data,
                              int block_rows)
 {
@@ -334,18 +334,18 @@ void isce3::signal::filter2D(isce3::io::Raster& output_raster,
     template void isce3::signal::filter2D<T>(                                  \
             isce3::io::Raster & output_raster,                                 \
             isce3::io::Raster & input_raster,                                  \
-            std::valarray<double> & kernel_columns,                            \
-            std::valarray<double> & kernel_rows, int block_rows);              \
+            const std::valarray<double> & kernel_columns,                      \
+            const std::valarray<double> & kernel_rows, int block_rows);        \
     template void isce3::signal::filter2D<T>(                                  \
             isce3::io::Raster & output_raster,                                 \
             isce3::io::Raster & input_raster, isce3::io::Raster & mask_raster, \
-            std::valarray<double> & kernel_columns,                            \
-            std::valarray<double> & kernel_rows, int block_rows);              \
+            const std::valarray<double> & kernel_columns,                      \
+            const std::valarray<double> & kernel_rows, int block_rows);        \
     template void isce3::signal::filter2D<T>(                                  \
             isce3::io::Raster & output_raster,                                 \
             isce3::io::Raster & input_raster, isce3::io::Raster & mask_raster, \
-            std::valarray<double> & kernel_columns,                            \
-            std::valarray<double> & kernel_rows, const bool do_decimate,       \
+            const std::valarray<double> & kernel_columns,                      \
+            const std::valarray<double> & kernel_rows, const bool do_decimate, \
             const bool mask, int block_rows)
 
 SPECIALIZE_FILTER(float);

cxx/isce3/signal/filter2D.h

@@ -16,8 +16,8 @@ namespace isce3 { namespace signal {
  */
 template<typename T>
 void filter2D(isce3::io::Raster& output_raster, isce3::io::Raster& input_raster,
-              std::valarray<double>& kernel_columns,
-              std::valarray<double>& kernel_rows, int block_rows = 1000);
+              const std::valarray<double>& kernel_columns,
+              const std::valarray<double>& kernel_rows, int block_rows = 1000);
 
 /**
  * filters real or complex data by convolving two 1D separable kernels in
@@ -32,8 +32,8 @@ void filter2D(isce3::io::Raster& output_raster, isce3::io::Raster& input_raster,
 template<typename T>
 void filter2D(isce3::io::Raster& output_raster, isce3::io::Raster& input_raster,
               isce3::io::Raster& mask_raster,
-              std::valarray<double>& kernel_columns,
-              std::valarray<double>& kernel_rows, int block_rows = 1000);
+              const std::valarray<double>& kernel_columns,
+              const std::valarray<double>& kernel_rows, int block_rows = 1000);
 
 /**
  * filters real or complex data by convolving two 1D separable kernels in
@@ -50,8 +50,8 @@ void filter2D(isce3::io::Raster& output_raster, isce3::io::Raster& input_raster,
 template<typename T>
 void filter2D(isce3::io::Raster& output_raster, isce3::io::Raster& input_raster,
               isce3::io::Raster& mask_raster,
-              std::valarray<double>& kernel_columns,
-              std::valarray<double>& kernel_rows, const bool do_decimate,
+              const std::valarray<double>& kernel_columns,
+              const std::valarray<double>& kernel_rows, const bool do_decimate,
               const bool mask = true, int block_rows = 1000);
 
 }} // namespace isce3::signal

python/extensions/pybind_isce3/Sources.cmake

@@ -56,6 +56,7 @@ signal/Covariance.cpp
 signal/Crossmul.cpp
 signal/CrossMultiply.cpp
 signal/flatten.cpp
+signal/filter2D.cpp
 product/GeoGridParameters.cpp
 product/product.cpp
 product/RadarGridParameters.cpp

python/extensions/pybind_isce3/signal/filter2D.cpp

@@ -0,0 +1,77 @@
+#include "filter2D.h"
+
+#include <gdal_priv.h>
+#include <valarray>
+#include <pybind11/stl.h>
+#include <pybind11/complex.h>
+
+#include <isce3/signal/filter2D.h>
+#include <isce3/except/Error.h>
+#include <isce3/io/Raster.h>
+
+namespace py = pybind11;
+
+using isce3::signal::filter2D;
+using isce3::io::Raster;
+
+void addbinding_filter2D(py::module& m)
+{
+    m
+    .def("filter2D", [](
+        Raster& output,
+        Raster& input,
+        const std::valarray<double> & kernel_columns,
+        const std::valarray<double> & kernel_rows,
+        int block_rows) {
+            int band=1;
+            auto in_type = input.dtype(band);
+            auto out_type = output.dtype(band);
+            // ensure input and output data type match
+            if (in_type != out_type)
+                throw isce3::except::RuntimeError(ISCE_SRCINFO(),
+                      "Input and output data type do not match");
+            // check kernel_columns dimensions
+            if (kernel_columns.size() > input.width())
+                throw isce3::except::RuntimeError(ISCE_SRCINFO(),
+                       "kernel column width > input width");
+            // check kernel_rows dimensions
+            if (kernel_rows.size() > input.length())
+               throw isce3::except::RuntimeError(ISCE_SRCINFO(),
+                       "kernel rows width > input length");
+            if (output.length() != input.length())
+                if (output.length() != input.length() / kernel_rows.size())
+                    throw isce3::except::RuntimeError(ISCE_SRCINFO(),
+                          "output length not equal to input length or not"
+                          "equal to input length divided by kernel row size");
+            if (output.width() != input.width())
+                if (output.width() != output.width() / kernel_columns.size())
+                    throw isce3::except::RuntimeError(ISCE_SRCINFO(),
+                          "output width not equal to input width or not"
+                          "equal to input width divided by kernel columns size");
+            // ensure types match
+            switch (in_type) {
+            case GDT_Float32: filter2D<float>(output, input, kernel_columns,
+                                      kernel_rows, block_rows);
+                                return;
+            case GDT_Float64: filter2D<double>(output, input, kernel_columns,
+                                      kernel_rows, block_rows);
+                                return;
+            case GDT_CFloat32: filter2D<std::complex<float>>(output, input,
+                                      kernel_columns, kernel_rows, block_rows);
+                                return;
+            case GDT_CFloat64: filter2D<std::complex<double>>(output, input,
+                                      kernel_columns, kernel_rows, block_rows);
+                                return;
+            default: throw isce3::except::RuntimeError(ISCE_SRCINFO(),
+                    "unsupported GDAL datatype");
+            }
+
+        },
+        py::arg("output"),
+        py::arg("input"),
+        py::arg("kernel_columns"),
+        py::arg("kernel_rows"),
+        py::arg("block_rows"),
+        "Filter real or complex data by convolving two separable 1D kernels")
+;
+}

python/extensions/pybind_isce3/signal/filter2D.h

@@ -0,0 +1,5 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+
+void addbinding_filter2D(pybind11::module& m);

python/extensions/pybind_isce3/signal/signal.cpp

@@ -5,6 +5,7 @@
 #include "CrossMultiply.h"
 #include "Crossmul.h"
 #include "flatten.h"
+#include "filter2D.h"
 
 namespace py = pybind11;
 
@@ -25,6 +26,7 @@ void addsubmodule_signal(py::module & m)
     addbinding(pyCrossmul);
     addbinding(pyCrossMultiply);
     addbinding_flatten(m_signal);
+    addbinding_filter2D(m_signal);
     addbinding_convolve2D<float>(m_signal);
     addbinding_convolve2D<std::complex<float>>(m_signal);
     addbinding_convolve2D<double>(m_signal);

tests/python/extensions/pybind/CMakeLists.txt

@@ -37,6 +37,7 @@ signal/convolve2D.py
 signal/covariance.py
 signal/crossmul.py
 signal/crossmultiply.py
+signal/filter2D.py
 product/radargridparameters.py
 product/swath.py
 unwrap/icu.py

tests/python/extensions/pybind/signal/filter2D.py

@@ -0,0 +1,104 @@
+import pybind_isce3 as isce3
+import numpy as np
+from scipy.signal import convolve2d
+
+from osgeo import gdal
+from osgeo import gdal_array
+
+
+def to_gdal_dataset(outpath, array):
+    driver = gdal.GetDriverByName("GTiff")
+    dtype = gdal_array.NumericTypeCodeToGDALTypeCode(array.dtype)
+    length, width = array.shape
+    dset = driver.Create(outpath, xsize=width, ysize=length, bands=1,
+                         eType=dtype)
+    dset.GetRasterBand(1).WriteArray(array)
+
+
+def open_raster(filepath):
+    ds = gdal.Open(filepath, gdal.GA_ReadOnly)
+    array = ds.GetRasterBand(1).ReadAsArray()
+
+    return array
+
+
+def test_run_filter2D():
+    # Data dimension
+    length = 200
+    width = 311
+
+    # Kernel dimensions
+    kernel_length = 3
+
+    block = 20
+
+    # Create filter kernels
+    kernel1d = np.ones([kernel_length, 1], dtype=np.float64) / kernel_length
+
+    # Create real data to filter
+    data_real = np.zeros([length, width], dtype=np.float64)
+    data_cpx = np.zeros([length, width], dtype=np.complex128)
+
+    for line in range(0, length):
+        for col in range(0, width):
+            data_real[line, col] = line + col
+            data_cpx[line, col] = np.cos(line * col) + 1.0j * np.sin(line * col)
+
+    # Save data
+    to_gdal_dataset('data.real', data_real)
+    to_gdal_dataset('data.cpx', data_cpx)
+
+    # Filter data
+    filt_data_real = isce3.io.Raster('data_real.filt', width, length, 1,
+                                     gdal.GDT_Float64, "ENVI")
+    filt_data_cpx = isce3.io.Raster('data_cpx.filt', width, length, 1,
+                                    gdal.GDT_CFloat64, "ENVI")
+
+    data_raster_real = isce3.io.Raster('data.real')
+    data_raster_cpx = isce3.io.Raster('data.cpx')
+
+    isce3.signal.filter2D(filt_data_real, data_raster_real, kernel1d,
+                          kernel1d, block)
+    isce3.signal.filter2D(filt_data_cpx, data_raster_cpx, kernel1d,
+                          kernel1d, block)
+
+
+def test_validate_filter_real():
+    # Create 2D kernel
+    kernel_length = 3
+    kernel_width = 3
+    kernel = np.ones((kernel_length, kernel_width), dtype=np.float64) / (
+            kernel_width * kernel_length)
+
+    # Validate filter2D for real data
+    data = open_raster('data.real')
+
+    filt_data = open_raster('data_real.filt')
+
+    # Regenerate filtered data
+    out = convolve2d(data, kernel, mode='same')
+
+    diff = np.abs(out - filt_data)
+    assert diff.max() < 1e-12
+
+
+def test_validate_filter_complex():
+    # Create kernel 2 D
+    kernel_width = 3
+    kernel_length = 3
+    kernel = np.ones((kernel_length, kernel_width), dtype=np.float64) / (
+            kernel_width * kernel_length)
+
+    # Open data
+    data = open_raster('data.cpx')
+
+    filt_data = open_raster('data_cpx.filt')
+
+    # Regenerate filter data
+    out = convolve2d(data, kernel, mode='same')
+
+    diff_pha = np.angle(filt_data * np.conj(out))
+    diff_amp = np.abs(filt_data) - np.abs(out)
+
+    assert diff_pha.max() < 1e-12
+    assert diff_amp.max() < 1e-12