Fast linear 2d interpolation

.clang-format

@@ -0,0 +1,10 @@
+BasedOnStyle: Google
+IndentWidth: 4
+AccessModifierOffset: -4
+ColumnLimit: 100
+AllowShortFunctionsOnASingleLine: None
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: false
+DerivePointerAlignment: true
+KeepEmptyLinesAtTheStartOfBlocks: true
+SortIncludes: false

.clang-tidy

@@ -0,0 +1,47 @@
+---
+Checks:          '-*,cppcoreguidelines-avoid-goto,cppcoreguidelines-pro-type-const-cast, google-readability-casting, google-runtime-int, modernize-replace-random-shuffle, modernize-use-nullptr, readability-braces-around-statements, readability-container-size-empty, readability-redundant-control-flow, readability-redundant-string-init, readability-identifier-naming, google-build-using-namespace'
+HeaderFilterRegex: '\.h$'
+WarningsAsErrors: '*'
+CheckOptions:
+  - key: readability-identifier-naming.NamespaceCase
+    value: lower_case
+  - key: readability-identifier-naming.ClassCase
+    value: CamelCase
+  - key: readability-identifier-naming.TypedefCase
+    value: CamelCase
+  - key: readability-identifier-naming.TypeAliasCase
+    value: CamelCase
+  - key: readability-identifier-naming.PrivateMemberSuffix
+    value: '_'
+  - key: readability-identifier-naming.StructCase
+    value: CamelCase
+  - key: readability-identifier-naming.FunctionCase
+    value: CamelCase
+  - key: readability-identifier-naming.VariableCase
+    value: lower_case
+  - key: readability-identifier-naming.PrivateMemberCase
+    value: lower_case
+  - key: readability-identifier-naming.ParameterCase
+    value: lower_case
+  - key: readability-identifier-naming.GlobalConstantPrefix
+    value: k
+  - key: readability-identifier-naming.GlobalConstantCase
+    value: CamelCase
+  - key: readability-identifier-naming.StaticConstantPrefix
+    value: k
+  - key: readability-identifier-naming.StaticConstantCase
+    value: CamelCase
+  - key: readability-identifier-naming.ConstexprVariableCase
+    value: CamelCase
+  - key: readability-identifier-naming.ConstexprVariablePrefix
+    value: k
+  - key: google-runtime-int.TypeSuffix
+    value: _t
+  - key: readability-identifier-naming.TypeTemplateParameterCase
+    value: CamelCase
+  - key: readability-identifier-naming.TypeTemplateParameterIgnoredRegexp
+    value: expr-type
+  - key: readability-identifier-naming.FunctionIgnoredRegexp
+    value: ^(begin|end)$
+  - key: readability-identifier-naming.TypeAliasIgnoredRegexp
+    value: ^(iterator_category|value_type|difference_type|pointer|reference|iterator_type|iterator_concept)$

.github/workflows/release.yml

@@ -50,8 +50,26 @@ jobs:
         run: |
           brew install llvm libomp
           echo $PATH
-          ln -sf /usr/local/bin/gcc-12 /usr/local/bin/gcc
+          ln -sf /usr/local/bin/gcc-11 /usr/local/bin/gcc
+          ln -sf /usr/local/bin/g++-11 /usr/local/bin/g++
           ls /usr/local/bin/gcc*
+          ls /usr/local/bin/g++*
+          gcc --version
+          g++ --version
+      - name: Install g++-11 for ubuntu
+        if: matrix.os == 'ubuntu-20.04'
+        id: install_cc
+        uses: rlalik/setup-cpp-compiler@master
+        with:
+          compiler: g++-11
+      - name: Check compilers for ubuntu
+        if: matrix.os == 'ubuntu-20.04'
+        run: |
+          ls /usr/bin/gcc*
+          ls /usr/bin/g++*
+          sudo ln -sf /usr/bin/gcc-11 /usr/bin/gcc
+          sudo ln -sf /usr/bin/g++-11 /usr/bin/g++
+          g++ --version
           gcc --version
       - name: Build wheels
         run: |

.github/workflows/test_build.yml

@@ -24,8 +24,26 @@ jobs:
         run: |
           brew install llvm libomp
           echo $PATH
-          ln -sf /usr/local/bin/gcc-12 /usr/local/bin/gcc
+          ln -sf /usr/local/bin/gcc-11 /usr/local/bin/gcc
+          ln -sf /usr/local/bin/g++-11 /usr/local/bin/g++
           ls /usr/local/bin/gcc*
+          ls /usr/local/bin/g++*
+          gcc --version
+          g++ --version
+      - name: Install g++-11 for ubuntu
+        if: matrix.os == 'ubuntu-20.04'
+        id: install_cc
+        uses: rlalik/setup-cpp-compiler@master
+        with:
+          compiler: g++-11
+      - name: Check compilers for ubuntu
+        if: matrix.os == 'ubuntu-20.04'
+        run: |
+          ls /usr/bin/gcc*
+          ls /usr/bin/g++*
+          sudo ln -sf /usr/bin/gcc-11 /usr/bin/gcc
+          sudo ln -sf /usr/bin/g++-11 /usr/bin/g++
+          g++ --version
           gcc --version
       - name: Build wheels
         run: |

.gitignore

@@ -7,3 +7,5 @@ imops/src/_fast*.pyx
 *.egg-info/
 .asv/
 dist/
+*.so
+.vscode/

MANIFEST.in

@@ -3,3 +3,4 @@ include requirements.txt
 include pyproject.toml
 include _build_utils.py
 recursive-include imops *.py
+recursive-include imops/cpp *.h *.hpp *.cpp

README.md

@@ -60,6 +60,24 @@ Works faster only for `ndim<=4, dtype=float32 or float64 (and bool-int16-32-64 i
 from imops import interp1d  # same as `scipy.interpolate.interp1d`
 ```
 Works faster only for `ndim<=3, dtype=float32 or float64, order=1`
+
+### Fast 2d linear interpolation
+```python
+import numpy as np
+from imops.interp2d import Linear2DInterpolator
+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:]
+num_threads = -1 # will be equal to num of CPU cores
+# You can optionally pass your own triangulation as an np.array of shape [num_triangles, 3], element at (i, j) position is an index of a point from x_points
+interpolator = Linear2DInterpolator(x_points, values, num_threads=num_threads, triangles=None)
+# 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)
+```
+
 ### Fast binary morphology
 
 ```python

_build_utils.py

@@ -6,12 +6,13 @@
 from setuptools.command.build_py import build_py
 
 
-class NumpyImport(dict):
-    """Hacky way to return Numpy's `include` path with lazy import."""
+class LazyImport(dict):
+    """Hacky way to return any module's `include` path with lazy import."""
 
     # Must be json-serializable due to
     # https://github.com/cython/cython/blob/6ad6ca0e9e7d030354b7fe7d7b56c3f6e6a4bc23/Cython/Compiler/ModuleNode.py#L773
-    def __init__(self):
+    def __init__(self, module_name):
+        self.module_name = module_name
         return super().__init__(self, description=self.__doc__)
 
     # Must be hashable due to
@@ -20,9 +21,10 @@ def __hash__(self):
         return id(self)
 
     def __repr__(self):
-        import numpy as np
+        scope = {}
+        exec(f'import {self.module_name} as module', scope)
 
-        return np.get_include()
+        return scope['module'].get_include()
 
     __fspath__ = __repr__
 
@@ -56,21 +58,25 @@ def get_ext_modules():
     name = 'imops'
     on_windows = platform.system() == 'Windows'
     args = ['/openmp' if on_windows else '-fopenmp']
+    cpp_args = [
+        '/std:c++20' if on_windows else '-std=c++2a',
+        '/O3' if on_windows else '-O3',
+    ]  # FIXME: account for higher gcc versions
 
     # Cython extension and .pyx source file names must be the same to compile
     # https://stackoverflow.com/questions/8024805/cython-compiled-c-extension-importerror-dynamic-module-does-not-define-init-fu
     modules = ['backprojection', 'measure', 'morphology', 'numeric', 'radon', 'zoom']
     modules_to_link_against_numpy_core_math_lib = ['numeric']
 
+    src_dir = Path(__file__).parent / name / 'src'
     for module in modules:
-        src_dir = Path(__file__).parent / name / 'src'
         shutil.copyfile(src_dir / f'_{module}.pyx', src_dir / f'_fast_{module}.pyx')
 
     return [
         Extension(
             f'{name}.src._{prefix}{module}',
             [f'{name}/src/_{prefix}{module}.pyx'],
-            include_dirs=[NumpyImport()],
+            include_dirs=[LazyImport('numpy')],
             library_dirs=[NumpyLibImport()] if module in modules_to_link_against_numpy_core_math_lib else [],
             libraries=['npymath'] + ['m'] * (not on_windows)
             if module in modules_to_link_against_numpy_core_math_lib
@@ -84,4 +90,13 @@ def get_ext_modules():
         # fallback to standard `-O2` compiled versions until https://github.com/neuro-ml/imops/issues/37 is resolved
         # for prefix, additional_args in zip(['', 'fast_'], [[], ['-ffast-math']])
         for prefix, additional_args in zip(['', 'fast_'], [[], []])
+    ] + [
+        Extension(
+            f'{name}.cpp.cpp_modules',
+            [f'{name}/cpp/main.cpp'],
+            include_dirs=[LazyImport('pybind11')],
+            extra_compile_args=args + cpp_args,
+            extra_link_args=args + cpp_args,
+            language='c++',
+        )
     ]

asv.conf.json

@@ -17,7 +17,8 @@
             "scikit-image": [],
             "numba": [],
             "pip+connected-components-3d": [],
-            "pip+fastremap": []
+            "pip+fastremap": [],
+            "pybind11": []
         },
         "env_nobuild": {
             "OMP_NUM_THREADS": "8",

benchmarks/benchmark_interp2d.py

@@ -0,0 +1,41 @@
+from pathlib import Path
+
+import numpy as np
+from scipy.interpolate import LinearNDInterpolator
+
+from imops.interp2d import Linear2DInterpolator
+
+from .common import NUMS_THREADS_TO_BENCHMARK, discard_arg, load_npy_gz
+
+
+class Interp2dSuite:
+    params = [('C++', 'Scipy'), ('float32', 'float64'), NUMS_THREADS_TO_BENCHMARK]
+    param_names = ['backend', 'dtype', 'num_threads']
+
+    @discard_arg(1)
+    @discard_arg(-1)
+    def setup(self, dtype):
+        x = load_npy_gz(Path(__file__).parent / 'data' / 'ribs.npy.gz').astype(dtype)
+        add_cols = x.shape[1] // 4
+        distances = np.concatenate((x[..., -add_cols:], x, x[..., :add_cols]), axis=1)
+        self.int_points = np.transpose((np.isnan(distances)).nonzero())
+        self.x_points = np.transpose((~np.isnan(distances)).nonzero())
+        self.x_values = distances[~np.isnan(distances)]
+
+    @discard_arg(2)
+    def time_interp2d(self, backend, num_threads):
+        if backend == 'C++':
+            Linear2DInterpolator(self.x_points, self.x_values, num_threads=num_threads)(self.int_points, fill_value=0.0)
+        elif backend == 'Scipy':
+            LinearNDInterpolator(self.x_points, self.x_values)(self.int_points)
+        else:
+            raise NotImplementedError
+
+    @discard_arg(2)
+    def peakmem_interp2d(self, backend, num_threads):
+        if backend == 'C++':
+            Linear2DInterpolator(self.x_points, self.x_values, num_threads=num_threads)(self.int_points, fill_value=0.0)
+        elif backend == 'Scipy':
+            LinearNDInterpolator(self.x_points, self.x_values)(self.int_points)
+        else:
+            raise NotImplementedError

docs/index.md

@@ -33,6 +33,11 @@ pip install imops[numba]  # additionally install Numba backend
       members:
         - __call__
 
+::: imops.interp2d.Linear2DInterpolator
+    options:
+      members:
+        - __call__
+
 ::: imops.morphology.binary_dilation
 
 ::: imops.morphology.binary_erosion

imops/__version__.py

@@ -1 +1 @@
-__version__ = '0.8.4'
+__version__ = '0.8.5'

imops/cpp/interp2d/delaunator/delaunator-header-only.hpp

@@ -0,0 +1,33 @@
+// MIT License
+
+// Copyright (c) 2018 Volodymyr Bilonenko
+
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+
+
+#pragma once
+
+#define DELAUNATOR_HEADER_ONLY
+
+#include "delaunator.hpp"
+
+#include "delaunator.cpp"
+
+#undef DELAUNATOR_HEADER_ONLY