diff --git a/imops/cpp/src/interp2d/interpolator.h b/imops/cpp/src/interp2d/interpolator.h index b816134b..0b7b9372 100644 --- a/imops/cpp/src/interp2d/interpolator.h +++ b/imops/cpp/src/interp2d/interpolator.h @@ -5,14 +5,14 @@ class Interpolator { public: using pyarr_double = py::array_t; const Triangulator triangulation; - Interpolator(const Triangulator::pyarr_size_t& pypoints, int n_jobs, std::optional pytriangles) + Interpolator(const Triangulator::pyarr_size_t& pypoints, int n_jobs, + std::optional pytriangles) : triangulation(pypoints, n_jobs, pytriangles) { - // if (pytriangles.has_value()) - // std::cout << pytriangles -> shape(0) << ' ' << pytriangles -> shape(1) << '\n'; } - pyarr_double operator()(const Triangulator::pyarr_size_t& int_points, const pyarr_double& values, - const Triangulator::pyarr_size_t& neighbors, double fill_value = 0.0) { + pyarr_double operator()(const Triangulator::pyarr_size_t& int_points, + const pyarr_double& values, const Triangulator::pyarr_size_t& neighbors, + double fill_value = 0.0) { if (triangulation.points.size() / 2 != values.shape(0)) { throw std::invalid_argument("Length mismatch between known points and their values"); } @@ -21,7 +21,7 @@ class Interpolator { } size_t n = int_points.shape()[0]; - std::vector int_values(n, fill_value); + std::vector int_values(n); omp_set_dynamic(0); // Explicitly disable dynamic teams omp_set_num_threads(triangulation.n_jobs_); @@ -42,6 +42,8 @@ class Interpolator { double f2 = values.at(triangulation.triangles.at(t + 1)); double f3 = values.at(triangulation.triangles.at(t + 2)); int_values[i] = (f1 * lambda_1 + f2 * lambda_2 + f3 * lambda_3) * one_over_2area; + } else { + int_values[i] = fill_value; } } diff --git a/imops/cpp/src/interp2d/triangulator.h b/imops/cpp/src/interp2d/triangulator.h index 1644e553..f8a1b77d 100644 --- a/imops/cpp/src/interp2d/triangulator.h +++ b/imops/cpp/src/interp2d/triangulator.h @@ -22,7 +22,8 @@ class Triangulator { std::vector> point2tri; std::unordered_map> edge2tri; - Triangulator(const pyarr_size_t& pypoints, int n_jobs, std::optional pytriangles) { + Triangulator(const pyarr_size_t& pypoints, int n_jobs, + std::optional pytriangles) { if (n_jobs < 0 and n_jobs != -1) { throw std::invalid_argument( "Invalid number of workers, have to be -1 or positive integer"); @@ -42,16 +43,16 @@ class Triangulator { points.at(j) = pypoints.at(i, 0); points.at(j + 1) = pypoints.at(i, 1); } - - size_t m = pytriangles -> shape(0); + + size_t m = pytriangles->shape(0); triangles.resize(3 * m); - + #pragma parallel for for (size_t i = 0; i < m; ++i) { size_t j = 3 * i; - triangles.at(j) = pytriangles -> at(i, 0); - triangles.at(j + 1) = pytriangles -> at(i, 1); - triangles.at(j + 2) = pytriangles -> at(i, 2); + triangles.at(j) = pytriangles->at(i, 0); + triangles.at(j + 1) = pytriangles->at(i, 1); + triangles.at(j + 2) = pytriangles->at(i, 2); } } diff --git a/imops/interp2d.py b/imops/interp2d.py index 8dc3c31b..a38d39b9 100644 --- a/imops/interp2d.py +++ b/imops/interp2d.py @@ -7,20 +7,20 @@ class Linear2DInterpolator(Linear2DInterpolatorCpp): - """ - 2D delaunay triangulation and parallel linear interpolation - Example: - n, m = 1024, 2 - points = np.random.randint(low=0, high=1024, size=(n, m)) - points = np.unique(points, axis=0) - x_points = points[: n // 2] - values = np.random.uniform(low=0.0, high=1.0, size=(len(x_points),)) - interp_points = points[n // 2:] - n_jobs = -1 # will be equal to num of CPU cores - interpolator = Linear2DInterpolator(x_points, values, n_jobs) - # Also you can pass values to __call__ and rewrite the ones that were passed to __init__ - interp_values = interpolator(interp_points, values + 1.0, fill_value=0.0) - """ + # """ + # #### 2D delaunay triangulation and parallel linear interpolation + # Example: + # n, m = 1024, 2 + # points = np.random.randint(low=0, high=1024, size=(n, m)) + # points = np.unique(points, axis=0) + # x_points = points[: n // 2] + # values = np.random.uniform(low=0.0, high=1.0, size=(len(x_points),)) + # interp_points = points[n // 2:] + # n_jobs = -1 // will be equal to num of CPU cores + # interpolator = Linear2DInterpolator(x_points, values, n_jobs) + # // Also you can pass values to __call__ and rewrite the ones that were passed to __init__ + # interp_values = interpolator(interp_points, values + 1.0, fill_value=0.0) + # """ def __init__(self, points: np.ndarray, values: np.ndarray, n_jobs: int = 1, triangles: Optional[np.ndarray] = None, **kwargs): super().__init__(points, n_jobs, triangles) diff --git a/tests/test_data/arr_0.npy b/tests/test_data/arr_0.npy new file mode 100644 index 00000000..1226003d Binary files /dev/null and b/tests/test_data/arr_0.npy differ diff --git a/tests/test_data/arr_1.npy b/tests/test_data/arr_1.npy new file mode 100644 index 00000000..c360f9c3 Binary files /dev/null and b/tests/test_data/arr_1.npy differ diff --git a/tests/test_data/arr_2.npy b/tests/test_data/arr_2.npy new file mode 100644 index 00000000..08df177e Binary files /dev/null and b/tests/test_data/arr_2.npy differ diff --git a/tests/test_data/arr_3.npy b/tests/test_data/arr_3.npy new file mode 100644 index 00000000..a1c6a097 Binary files /dev/null and b/tests/test_data/arr_3.npy differ diff --git a/tests/test_data/arr_4.npy b/tests/test_data/arr_4.npy new file mode 100644 index 00000000..b81be593 Binary files /dev/null and b/tests/test_data/arr_4.npy differ diff --git a/tests/test_data/arr_5.npy b/tests/test_data/arr_5.npy new file mode 100644 index 00000000..b60c50a1 Binary files /dev/null and b/tests/test_data/arr_5.npy differ diff --git a/tests/test_data/arr_6.npy b/tests/test_data/arr_6.npy new file mode 100644 index 00000000..312b01b8 Binary files /dev/null and b/tests/test_data/arr_6.npy differ diff --git a/tests/test_data/val_2.npy b/tests/test_data/val_2.npy new file mode 100644 index 00000000..2ce3a380 Binary files /dev/null and b/tests/test_data/val_2.npy differ diff --git a/tests/test_data/val_3.npy b/tests/test_data/val_3.npy new file mode 100644 index 00000000..dd194b5a Binary files /dev/null and b/tests/test_data/val_3.npy differ diff --git a/tests/test_data/val_4.npy b/tests/test_data/val_4.npy new file mode 100644 index 00000000..2450b328 Binary files /dev/null and b/tests/test_data/val_4.npy differ diff --git a/tests/test_data/val_5.npy b/tests/test_data/val_5.npy new file mode 100644 index 00000000..8c9bca57 Binary files /dev/null and b/tests/test_data/val_5.npy differ diff --git a/tests/test_data/val_6.npy b/tests/test_data/val_6.npy new file mode 100644 index 00000000..00ad14b9 Binary files /dev/null and b/tests/test_data/val_6.npy differ diff --git a/tests/test_interp2d.py b/tests/test_interp2d.py new file mode 100644 index 00000000..93fa0021 --- /dev/null +++ b/tests/test_interp2d.py @@ -0,0 +1,42 @@ +import numpy as np +import scipy +from scipy.interpolate import griddata + +from imops.interp2d import Linear2DInterpolator +from deli import load_numpy + +np.random.seed(1337) + + +def test_test_data(): + for i in range(2): + x = load_numpy(f'test_data/arr_{i}.npy') + int_points = np.transpose((np.isnan(x)).nonzero()) + x_points: np.ndarray = np.transpose((~np.isnan(x)).nonzero()) + x_values = x[~np.isnan(x)] + imops_values = Linear2DInterpolator(x_points, x_values, n_jobs=1)(int_points, fill_value=0.0) + triangles = scipy.spatial.Delaunay(x_points).simplices + delaunay_values = Linear2DInterpolator(x_points, x_values, n_jobs=1, triangles=triangles)(int_points, fill_value=0.0) + scipy_values = griddata(x_points, x_values, int_points, method='linear', fill_value=0.0) + + delta_ds = np.abs(delaunay_values - scipy_values) + # delta_di = np.abs(delaunay_values - imops_values) + delta_si = np.abs(scipy_values - imops_values) + + assert delta_ds.max() <= 1e-8 and delta_si.max() <= 5, 'Failed with big cases, arr_0 or arr_1' + + for i in range(2, 7): + x = load_numpy(f'test_data/arr_{i}.npy') + x_values = load_numpy(f'test_data/val_{i}.npy') + x_points = np.transpose(x.nonzero()) + int_points = np.transpose((~x).nonzero()) + imops_values = Linear2DInterpolator(x_points, x_values, n_jobs=1)(int_points, fill_value=0.0) + triangles = scipy.spatial.Delaunay(x_points).simplices + delaunay_values = Linear2DInterpolator(x_points, x_values, n_jobs=1, triangles=triangles)(int_points, fill_value=0.0) + scipy_values = griddata(x_points, x_values, int_points, method='linear', fill_value=0.0) + + delta_ds = np.abs(delaunay_values - scipy_values) + delta_di = np.abs(delaunay_values - imops_values) + delta_si = np.abs(scipy_values - imops_values) + + assert delta_di.max() <= 1.5 and delta_ds.max() <= 1e-10 and delta_si.max() <= 1.5, 'Failed with small cases, arr_2 ...'