Skip to content

Commit 9ea3583

Browse files
Liang YuGitHub Enterprise
Liang Yu
authored and
GitHub Enterprise
committed
Pybind convolve2d decimation with EArrays (#715)
pybind convolve2D with decimate toggle with unit tests
1 parent 7ae0bd7 commit 9ea3583

File tree

6 files changed

+298
-0
lines changed

6 files changed

+298
-0
lines changed

‎python/extensions/pybind_isce3/Sources.cmake‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,7 @@ io/Raster.cpp
5151
io/serialization.cpp
5252
io/io.cpp
5353
signal/signal.cpp
54+
signal/convolve2D.cpp
5455
signal/Covariance.cpp
5556
signal/Crossmul.cpp
5657
signal/CrossMultiply.cpp

‎python/extensions/pybind_isce3/signal/convolve2D.cpp‎

Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
#include "convolve2D.h"
2+
3+
#include <isce3/except/Error.h>
4+
#include <isce3/signal/convolve.h>
5+
6+
#include <complex>
7+
8+
#include <pybind11/complex.h>
9+
#include <pybind11/eigen.h>
10+
#include <pybind11/numpy.h>
11+
12+
namespace py = pybind11;
13+
14+
using isce3::signal::convolve2D;
15+
using isce3::core::EArray2D;
16+
17+
template<typename T>
18+
void addbinding_convolve2D(py::module& m)
19+
{
20+
m
21+
.def("convolve2D", [](
22+
const isce3::core::EArray2D<T>& input,
23+
const isce3::core::EArray2D<double>& kernel_columns,
24+
const isce3::core::EArray2D<double>& kernel_rows,
25+
const bool decimate)
26+
{
27+
// get input dimensions
28+
auto input_length = input.rows();
29+
auto input_width = input.cols();
30+
31+
// check kernel_columns dimensions
32+
auto kernel_cols_length = kernel_columns.rows();
33+
auto kernel_cols_width = kernel_columns.cols();
34+
if (kernel_cols_length != 1)
35+
throw isce3::except::RuntimeError(ISCE_SRCINFO(),
36+
"kernel column not 1 x N");
37+
if (kernel_cols_width > input_width)
38+
throw isce3::except::RuntimeError(ISCE_SRCINFO(),
39+
"kernel column width > input width");
40+
41+
// check kernel_rows dimensions
42+
auto kernel_rows_width = kernel_rows.cols();
43+
if (kernel_rows_width != 1)
44+
throw isce3::except::RuntimeError(ISCE_SRCINFO(),
45+
"kernel column not N x 1");
46+
if (kernel_rows_width > input_width)
47+
throw isce3::except::RuntimeError(ISCE_SRCINFO(),
48+
"kernel rows width > input length");
49+
50+
// determine pad size
51+
// divide 2 in () to ensure integer division done 1st
52+
auto length_pad = 2 * (kernel_rows.size() / 2);
53+
auto width_pad = 2 * (kernel_columns.size() / 2);
54+
55+
// get output dimensions
56+
auto output_length = input_length - length_pad;
57+
auto output_width = input_width - width_pad;
58+
if (decimate) {
59+
output_length /= kernel_rows.size();
60+
output_width /= kernel_columns.size();
61+
}
62+
63+
// init weight of ones
64+
isce3::core::EArray2D<double> weights;
65+
weights.setOnes(input_length, input_width);
66+
67+
// init output
68+
isce3::core::EArray2D<T> output(output_length, output_width);
69+
70+
// convolve/decimate
71+
convolve2D<T>(output, input, weights, kernel_columns, kernel_rows);
72+
73+
return output;
74+
},
75+
py::arg("input"),
76+
py::arg("kernel_columns"),
77+
py::arg("kernel_rows"),
78+
py::arg("decimate"),
79+
"2D convolution in time domain with separable kernels")
80+
.def("convolve2D", [](
81+
const isce3::core::EArray2D<T>& input,
82+
const isce3::core::EArray2D<double>& weights,
83+
const isce3::core::EArray2D<double>& kernel_columns,
84+
const isce3::core::EArray2D<double>& kernel_rows,
85+
const bool decimate)
86+
{
87+
// get input dimensions
88+
auto input_length = input.rows();
89+
auto input_width = input.cols();
90+
91+
// get weight dimensions
92+
auto weight_length = weights.rows();
93+
auto weight_width = weights.cols();
94+
95+
// ensure input and weight dimensions match
96+
if (input_length != weight_length || input_width != weight_width)
97+
throw isce3::except::RuntimeError(ISCE_SRCINFO(),
98+
"input and weight dimensions do not match");
99+
100+
// check kernel_columns dimensions
101+
auto kernel_cols_length = kernel_columns.rows();
102+
auto kernel_cols_width = kernel_columns.cols();
103+
if (kernel_cols_length != 1)
104+
throw isce3::except::RuntimeError(ISCE_SRCINFO(),
105+
"kernel column not 1 x N");
106+
if (kernel_cols_width > input_width)
107+
throw isce3::except::RuntimeError(ISCE_SRCINFO(),
108+
"kernel column width > input width");
109+
110+
// check kernel_rows dimensions
111+
auto kernel_rows_width = kernel_rows.cols();
112+
if (kernel_rows_width != 1)
113+
throw isce3::except::RuntimeError(ISCE_SRCINFO(),
114+
"kernel column not N x 1");
115+
if (kernel_rows_width > input_width)
116+
throw isce3::except::RuntimeError(ISCE_SRCINFO(),
117+
"kernel rows width > input length");
118+
119+
// determine pad size
120+
auto length_pad = 2 * (kernel_rows.size() / 2);
121+
auto width_pad = 2 * (kernel_columns.size() / 2);
122+
123+
// determine output dimensions
124+
auto output_length = input_length - length_pad;
125+
auto output_width = input_width - width_pad;
126+
if (decimate) {
127+
output_length /= kernel_rows.size();
128+
output_width /= kernel_columns.size();
129+
}
130+
131+
// init output
132+
isce3::core::EArray2D<T> output(output_length, output_width);
133+
134+
// convolve/decimate
135+
convolve2D<T>(output, input, weights, kernel_columns, kernel_rows);
136+
137+
return output;
138+
},
139+
py::arg("input"),
140+
py::arg("weights"),
141+
py::arg("kernel_columns"),
142+
py::arg("kernel_rows"),
143+
py::arg("decimate"),
144+
"2D convolution in time domain with separable kernels and mask/weights")
145+
;
146+
}
147+
148+
template void addbinding_convolve2D<float>(py::module& m);
149+
template void addbinding_convolve2D<double>(py::module& m);
150+
template void addbinding_convolve2D<std::complex<float>>(py::module& m);
151+
template void addbinding_convolve2D<std::complex<double>>(py::module& m);

‎python/extensions/pybind_isce3/signal/convolve2D.h‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
#pragma once
2+
3+
#include <pybind11/pybind11.h>
4+
5+
template<typename T>
6+
void addbinding_convolve2D(pybind11::module& m);

‎python/extensions/pybind_isce3/signal/signal.cpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#include "signal.h"
22

3+
#include "convolve2D.h"
34
#include "Covariance.h"
45
#include "CrossMultiply.h"
56
#include "Crossmul.h"
@@ -24,4 +25,8 @@ void addsubmodule_signal(py::module & m)
2425
addbinding(pyCrossmul);
2526
addbinding(pyCrossMultiply);
2627
addbinding_flatten(m_signal);
28+
addbinding_convolve2D<float>(m_signal);
29+
addbinding_convolve2D<std::complex<float>>(m_signal);
30+
addbinding_convolve2D<double>(m_signal);
31+
addbinding_convolve2D<std::complex<double>>(m_signal);
2732
}

‎tests/python/extensions/pybind/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ image/resamp.py
3333
io/gdal/dataset.py
3434
io/gdal/raster.py
3535
io/raster.py
36+
signal/convolve2D.py
3637
signal/covariance.py
3738
signal/crossmul.py
3839
signal/crossmultiply.py

‎tests/python/extensions/pybind/signal/convolve2D.py‎

Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,134 @@
1+
import numpy as np
2+
import numpy.testing as npt
3+
4+
import pybind_isce3 as isce3
5+
6+
def get_dims():
7+
'''
8+
Return data and kernel dimenstions
9+
'''
10+
length = 200
11+
width = 310
12+
kernel_length = 3
13+
kernel_width = 3
14+
return length, width, kernel_length, kernel_width
15+
16+
def make_inputs(length, width, kernel_length, kernel_width):
17+
'''
18+
Create real and imag data
19+
'''
20+
# Calculate padding
21+
pad_cols = kernel_width - 1
22+
pad_rows = kernel_length - 1
23+
width_pad = width + pad_cols
24+
length_pad = length + pad_rows
25+
26+
# Populate real and imag data
27+
data_real = np.zeros([length_pad, width_pad], dtype=np.float64)
28+
data_imag = np.zeros([length_pad, width_pad], dtype=np.complex128)
29+
for line in range(pad_rows//2, length + pad_rows//2):
30+
for col in range(pad_cols//2, width + pad_cols//2):
31+
data_real[line, col] = line + col
32+
data_imag[line, col] = np.cos(line * col) + np.sin(line * col) * 1.0j
33+
34+
return data_real, data_imag
35+
36+
def make_expected_output(out_shape, data, kernel_width, kernel_length):
37+
'''
38+
Calculate expected decimated output
39+
'''
40+
decimated = np.zeros(out_shape, dtype=data.dtype)
41+
for i in range(out_shape[0]):
42+
for j in range(out_shape[1]):
43+
decimated[i,j] = np.mean(data[i*kernel_width+1:(i+1)*kernel_width+1,\
44+
j*kernel_length+1:(j+1)*kernel_length+1])
45+
return decimated
46+
47+
def test_ea_convolve2d_with_mask():
48+
'''
49+
Test convolve2D without mask
50+
'''
51+
length, width, kernel_length, kernel_width = get_dims()
52+
53+
# Create data
54+
input_real, input_imag = make_inputs(length, width, kernel_length, kernel_width)
55+
56+
# Create mask
57+
mask = np.ones(input_real.shape, dtype=np.float64)
58+
59+
# Create kernels
60+
kernel_cols = np.ones([1, kernel_width], dtype=np.float64)/kernel_width
61+
kernel_rows = np.ones([kernel_length, 1], dtype=np.float64)/kernel_length
62+
63+
# Convolve
64+
pybind_real = isce3.signal.convolve2D(input_real, mask, kernel_cols, kernel_rows, True)
65+
pybind_imag = isce3.signal.convolve2D(input_imag, mask, kernel_cols, kernel_rows, True)
66+
67+
# Calculate expected output
68+
out_shape = (length//kernel_width, width//kernel_length)
69+
expected_real = make_expected_output(out_shape, input_real,
70+
kernel_length, kernel_width)
71+
expected_imag = make_expected_output(out_shape, input_imag,
72+
kernel_length, kernel_width)
73+
74+
# Check outputs
75+
npt.assert_allclose(pybind_real, expected_real, rtol=0.0, atol=1e-12)
76+
npt.assert_allclose(np.angle(pybind_imag), np.angle(expected_imag),
77+
rtol=0.0, atol=1e-12)
78+
79+
def test_ea_convolve2d_no_mask():
80+
'''
81+
Test convolve2D without mask
82+
'''
83+
length, width, kernel_length, kernel_width = get_dims()
84+
85+
# Create data
86+
input_real, input_imag = make_inputs(length, width, kernel_length, kernel_width)
87+
88+
# Create kernels
89+
kernel_cols = np.ones([1, kernel_width], dtype=np.float64)/kernel_width
90+
kernel_rows = np.ones([kernel_length, 1], dtype=np.float64)/kernel_length
91+
92+
# Convolve
93+
pybind_real = isce3.signal.convolve2D(input_real, kernel_cols, kernel_rows, True)
94+
pybind_imag = isce3.signal.convolve2D(input_imag, kernel_cols, kernel_rows, True)
95+
96+
# Calculate expected output
97+
out_shape = (length//kernel_width, width//kernel_length)
98+
expected_real = make_expected_output(out_shape, input_real,
99+
kernel_length, kernel_width)
100+
expected_imag = make_expected_output(out_shape, input_imag,
101+
kernel_length, kernel_width)
102+
103+
# Check outputs
104+
npt.assert_allclose(pybind_real, expected_real, rtol=0.0, atol=1e-12)
105+
npt.assert_allclose(np.real(pybind_imag), np.real(expected_imag),
106+
rtol=0.0, atol=1e-12)
107+
npt.assert_allclose(np.imag(pybind_imag), np.imag(expected_imag),
108+
rtol=0.0, atol=1e-12)
109+
110+
def test_scipy():
111+
112+
from scipy import signal
113+
114+
data = np.ones((8,13))
115+
kernel = np.ones((3,3))/9.0
116+
117+
filt_data_scipy = signal.convolve2d(data, kernel, mode='same')
118+
119+
data_padded = np.ones((10,15))
120+
data_padded[0,:] = 0
121+
data_padded[:,0] = 0
122+
data_padded[:,-1] = 0
123+
data_padded[-1,:] = 0
124+
125+
kernel_cols = np.ones((1,3))/3.0
126+
kernel_rows = np.ones((3,1))/3.0
127+
filt_data_isce3 = isce3.signal.convolve2D(data_padded, kernel_cols, kernel_rows, False)
128+
129+
npt.assert_allclose(filt_data_isce3, filt_data_scipy, rtol=0.0, atol=1e-12)
130+
131+
132+
if __name__ == "__main__":
133+
test_ea_convolve2d_no_mask()
134+
test_ea_convolve2d_with_mask()

0 commit comments

Comments
 (0)