Merge pull request #31 from neuro-ml/dev

Gigafast morphology

README.md

@@ -5,7 +5,9 @@
 
 # Imops
 
-Efficient parallelizable algorithms for multidimensional arrays to speed up your data pipelines. Docs are [here](https://neuro-ml.github.io/imops/).
+Efficient parallelizable algorithms for multidimensional arrays to speed up your data pipelines.
+- [Documentation](https://neuro-ml.github.io/imops/)
+- [Benchmarks](https://neuro-ml.github.io/imops/benchmarks/)
 
 # Install
 
@@ -14,15 +16,33 @@ pip install imops  # default install with Cython backend
 pip install imops[numba]  # additionally install Numba backend
 ```
 
+# How fast is it?
+
+Time comparisons (ms) for Intel(R) Xeon(R) Silver 4114 CPU @ 2.20GHz using 8 threads. All inputs are C-contiguous NumPy arrays. For morphology functions `bool` dtype is used and `float64` for all others.
+| function / backend   |  Scipy()  |  Cython(fast=False)  |  Cython(fast=True)  |  Numba()  |
+|:----------------------:|:-----------:|:----------------------:|:---------------------:|:-----------:|
+| `zoom(..., order=0)` |   2072    |         1114         |         **867**         |   3590    |
+| `zoom(..., order=1)` |   6527    |         596          |         **575**         |   3757    |
+| `interp1d`           |    780    |         149          |         **146**         |    420    |
+| `radon`              |   59711   |         5982         |        **4837**         |      -     |
+| `inverse_radon`      |   52928   |         8254         |        **6535**         |         -  |
+| `binary_dilation`    |   2207    |         310          |         **298**         |        -   |
+| `binary_erosion`     |   2296    |         326          |         **304**         |        -   |
+| `binary_closing`     |   4158    |         544          |         **469**         |        -   |
+| `binary_opening`     |   4410    |         567          |         **522**         |        -   |
+| `center_of_mass`     |   2237    |          **64**          |         **64**          |        -   |
+
+We use [`airspeed velocity`](https://asv.readthedocs.io/en/stable/) to benchmark our code. For detailed results visit [benchmark page](https://neuro-ml.github.io/imops/benchmarks/).
+
 # Features
 
-## Fast Radon transform
+### Fast Radon transform
 
 ```python
 from imops import radon, inverse_radon
 ```
 
-## Fast linear/bilinear/trilinear zoom
+### Fast 0/1-order zoom
 
 ```python
 from imops import zoom, zoom_to_shape
@@ -33,20 +53,20 @@ y = zoom(x, 2, axis=[0, 1])
 # without the need to compute the scale factor
 z = zoom_to_shape(x, (4, 120, 67))
 ```
-Works faster only for `ndim<=3, dtype=float32 or float64, output=None, order=1, mode='constant', grid_mode=False`
-## Fast 1d linear interpolation
+Works faster only for `ndim<=4, dtype=float32 or float64 (and bool-int16-32-64 if order == 0), output=None, order=0 or 1, mode='constant', grid_mode=False`
+### Fast 1d linear interpolation
 
 ```python
 from imops import interp1d  # same as `scipy.interpolate.interp1d`
 ```
-Works faster only for `ndim<=3, dtype=float32 or float64, order=1 or 'linear'`
-## Fast binary morphology
+Works faster only for `ndim<=3, dtype=float32 or float64, order=1`
+### Fast binary morphology
 
 ```python
 from imops import binary_dilation, binary_erosion, binary_opening, binary_closing
 ```
 These functions mimic `scikit-image` counterparts
-## Padding
+### Padding
 
 ```python
 from imops import pad, pad_to_shape
@@ -59,7 +79,7 @@ y = pad(x, 10, axis=[0, 1])
 z = pad_to_shape(x, (4, 120, 67), ratio=0.25)
 ```
 
-## Cropping
+### Cropping
 
 ```python
 from imops import crop_to_shape
@@ -71,7 +91,7 @@ from imops import crop_to_shape
 z = crop_to_shape(x, (4, 120, 67), ratio=0.25)
 ```
 
-## Labeling
+### Labeling
 
 ```python
 from imops import label
@@ -81,7 +101,7 @@ labeled, num_components = label(x, background=1, return_num=True)
 ```
 
 # Backends
-For `zoom`, `zoom_to_shape`, `interp1d`, `radon`, `inverse_radon` you can specify which backend to use. Backend can be specified by a string or by an instance of `Backend` class. The latter allows you to customize some backend options:
+For all heavy image routines except `label` you can specify which backend to use. Backend can be specified by a string or by an instance of `Backend` class. The latter allows you to customize some backend options:
 ```python
 from imops import Cython, Numba, Scipy, zoom
 
@@ -102,17 +122,17 @@ with imops_backend('Cython'):  # sets Cython backend via context manager
 ```
 Note that for `Numba` backend setting `num_threads` argument has no effect for now and you should use `NUMBA_NUM_THREADS` environment variable.
 Available backends:
-|                   | Scipy   | Cython  | Numba   |
-|-------------------|---------|---------|---------|
+|         function / backend            | Scipy   | Cython  | Numba   |
+|:-------------------:|:---------:|:---------:|:---------:|
 | `zoom`            | &check; | &check; | &check; |
-| `zoom_to_shape`   | &check; | &check; | &check; |
 | `interp1d`        | &check; | &check; | &check; |
 | `radon`           | &cross; | &check; | &cross; |
 | `inverse_radon`   | &cross; | &check; | &cross; |
 | `binary_dilation` | &check; | &check; | &cross; |
 | `binary_erosion`  | &check; | &check; | &cross; |
 | `binary_closing`  | &check; | &check; | &cross; |
 | `binary_opening`  | &check; | &check; | &cross; |
+| `center_of_mass`  | &check; | &check; | &cross; |
 
 # Acknowledgements
 

_pyproject_build.py

@@ -5,7 +5,19 @@
 from setuptools.command.build_py import build_py
 
 
-class NumpyImport:
+class NumpyImport(dict):
+    """Hacky way to return Numpy'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):
+        return super().__init__(self, description=self.__doc__)
+
+    # Must be hashable due to
+    # https://github.com/cython/cython/blob/6ad6ca0e9e7d030354b7fe7d7b56c3f6e6a4bc23/Cython/Compiler/Main.py#L307
+    def __hash__(self):
+        return id(self)
+
     def __repr__(self):
         import numpy as np
 

asv.conf.json

@@ -11,7 +11,7 @@
     "matrix": {
         "req": {
             "numpy": [],
-            "Cython": [],
+            "Cython": ["0.29.36"],
             "scipy": [],
             "scikit-image": [],
             "numba": [],

benchmarks/benchmark_morphology.py

@@ -1,3 +1,5 @@
+from pathlib import Path
+
 import numpy as np
 
 
@@ -14,47 +16,91 @@
 
 from imops.morphology import binary_closing, binary_dilation, binary_erosion, binary_opening
 
-from .common import NUMS_THREADS_TO_BENCHMARK, discard_arg
+from .common import IMAGE_TYPE_BENCHMARK, IMAGE_TYPES_BENCHMARK, NUMS_THREADS_TO_BENCHMARK, discard_arg, load_npy_gz
 
 
 class MorphologySuite:
-    params = [morphology_configs, ('bool', 'int64'), NUMS_THREADS_TO_BENCHMARK]
-    param_names = ['backend', 'dtype', 'num_threads']
+    params = [morphology_configs, ('bool', 'int64'), NUMS_THREADS_TO_BENCHMARK, IMAGE_TYPES_BENCHMARK, (False, True)]
+    param_names = ['backend', 'dtype', 'num_threads', 'image_type', 'boxed']
 
     @discard_arg(1)
+    @discard_arg(2)
+    @discard_arg(-1)
     @discard_arg(-1)
     def setup(self, dtype):
-        self.image = np.random.randint(0, 5 if dtype == 'int64' else 2, (256, 256, 256)).astype(dtype)
+        real_images_path = Path(__file__).parent / 'data'
+
+        lungs_image = load_npy_gz(real_images_path / 'lungs.npy.gz').astype(dtype, copy=False)
+        bronchi_image = load_npy_gz(real_images_path / 'bronchi.npy.gz').astype(dtype, copy=False)
+        rand_image = np.random.randint(0, 5 if dtype == 'int64' else 2, (512, 512, 512)).astype(dtype, copy=False)
+
+        self.images = {
+            IMAGE_TYPE_BENCHMARK.RAND: rand_image,
+            IMAGE_TYPE_BENCHMARK.LUNGS: lungs_image,
+            IMAGE_TYPE_BENCHMARK.BRONCHI: bronchi_image,
+        }
 
     # FIXME: generalize this code somehow
     @discard_arg(2)
-    def time_closing(self, backend, num_threads):
-        binary_closing(self.image, num_threads=num_threads, backend=backend)
+    def time_erosion(self, backend, num_threads, image_type, boxed):
+        im = self.images[image_type]
+        try:
+            binary_erosion(im, num_threads=num_threads, backend=backend, boxed=boxed)
+        except TypeError:
+            binary_erosion(im, num_threads=num_threads, backend=backend)
 
     @discard_arg(2)
-    def time_dilation(self, backend, num_threads):
-        binary_dilation(self.image, num_threads=num_threads, backend=backend)
+    def time_dilation(self, backend, num_threads, image_type, boxed):
+        im = self.images[image_type]
+        try:
+            binary_dilation(im, num_threads=num_threads, backend=backend, boxed=boxed)
+        except TypeError:
+            binary_dilation(im, num_threads=num_threads, backend=backend)
 
     @discard_arg(2)
-    def time_erosion(self, backend, num_threads):
-        binary_erosion(self.image, num_threads=num_threads, backend=backend)
+    def time_opening(self, backend, num_threads, image_type, boxed):
+        im = self.images[image_type]
+        try:
+            binary_opening(im, num_threads=num_threads, backend=backend, boxed=boxed)
+        except TypeError:
+            binary_opening(im, num_threads=num_threads, backend=backend)
 
     @discard_arg(2)
-    def time_opening(self, backend, num_threads):
-        binary_opening(self.image, num_threads=num_threads, backend=backend)
+    def time_closing(self, backend, num_threads, image_type, boxed):
+        im = self.images[image_type]
+        try:
+            binary_closing(im, num_threads=num_threads, backend=backend, boxed=boxed)
+        except TypeError:
+            binary_closing(im, num_threads=num_threads, backend=backend)
 
     @discard_arg(2)
-    def peakmem_closing(self, backend, num_threads):
-        binary_closing(self.image, num_threads=num_threads, backend=backend)
+    def peakmem_erosion(self, backend, num_threads, image_type, boxed):
+        im = self.images[image_type]
+        try:
+            binary_erosion(im, num_threads=num_threads, backend=backend, boxed=boxed)
+        except TypeError:
+            binary_erosion(im, num_threads=num_threads, backend=backend)
 
     @discard_arg(2)
-    def peakmem_dilation(self, backend, num_threads):
-        binary_dilation(self.image, num_threads=num_threads, backend=backend)
+    def peakmem_dilation(self, backend, num_threads, image_type, boxed):
+        im = self.images[image_type]
+        try:
+            binary_dilation(im, num_threads=num_threads, backend=backend, boxed=boxed)
+        except TypeError:
+            binary_dilation(im, num_threads=num_threads, backend=backend)
 
     @discard_arg(2)
-    def peakmem_erosion(self, backend, num_threads):
-        binary_erosion(self.image, num_threads=num_threads, backend=backend)
+    def peakmem_opening(self, backend, num_threads, image_type, boxed):
+        im = self.images[image_type]
+        try:
+            binary_opening(im, num_threads=num_threads, backend=backend, boxed=boxed)
+        except TypeError:
+            binary_opening(im, num_threads=num_threads, backend=backend)
 
     @discard_arg(2)
-    def peakmem_opening(self, backend, num_threads):
-        binary_opening(self.image, num_threads=num_threads, backend=backend)
+    def peakmem_closing(self, backend, num_threads, image_type, boxed):
+        im = self.images[image_type]
+        try:
+            binary_closing(im, num_threads=num_threads, backend=backend, boxed=boxed)
+        except TypeError:
+            binary_closing(im, num_threads=num_threads, backend=backend)

benchmarks/common.py

@@ -1,4 +1,16 @@
-NUMS_THREADS_TO_BENCHMARK = list(range(1, 9))
+from enum import Enum
+from gzip import GzipFile
+
+from numpy import load
+
+
+class IMAGE_TYPE_BENCHMARK(Enum):
+    RAND = 1
+    LUNGS = 2
+    BRONCHI = 3
+
+    def __repr__(self):
+        return self.name
 
 
 def discard_arg(idx: int):
@@ -12,3 +24,12 @@ def wrapper(*args):
         return wrapper
 
     return inner
+
+
+def load_npy_gz(path):
+    with GzipFile(path, 'rb') as f:
+        return load(f)
+
+
+IMAGE_TYPES_BENCHMARK = list(IMAGE_TYPE_BENCHMARK)
+NUMS_THREADS_TO_BENCHMARK = list(range(1, 9))

imops/__version__.py

@@ -1 +1 @@
-__version__ = '0.8.1'
+__version__ = '0.8.2'

imops/box.py

@@ -1,7 +1,7 @@
 import itertools
 from copy import copy
 from functools import wraps
-from typing import Callable, Sequence, Tuple
+from typing import Callable, Tuple
 
 import numpy as np
 
@@ -10,34 +10,6 @@
 Box = np.ndarray
 
 
-def check_len(*args) -> None:
-    lengths = list(map(len, args))
-    if any(length != lengths[0] for length in lengths):
-        raise ValueError(f'Arguments of equal length are required: {", ".join(map(str, lengths))}')
-
-
-def build_slices(start: Sequence[int], stop: Sequence[int] = None, step: Sequence[int] = None) -> Tuple[slice, ...]:
-    """
-    Returns a tuple of slices built from `start` and `stop` with `step`.
-
-    Examples
-    --------
-    >>> build_slices([1, 2, 3], [4, 5, 6])
-    (slice(1, 4), slice(2, 5), slice(3, 6))
-    >>> build_slices([10, 11])
-    (slice(10), slice(11))
-    """
-
-    check_len(*filter(lambda x: x is not None, [start, stop, step]))
-    args = [
-        start,
-        stop if stop is not None else [None for _ in start],
-        step if step is not None else [None for _ in start],
-    ]
-
-    return tuple(map(slice, *args))
-
-
 def make_box(iterable) -> Box:
     """Returns a box, generated from copy of the `iterable`."""
     box = np.asarray(copy(iterable))
@@ -78,3 +50,22 @@ def mask_to_box(mask: np.ndarray) -> Box:
         stop.insert(0, right + 1)
 
     return start, stop
+
+
+@returns_box
+def shape_to_box(shape: Tuple) -> Box:
+    return make_box([(0,) * len(shape), shape])  # fmt: skip
+
+
+def box_to_shape(box: Box) -> Tuple:
+    return tuple(box[1] - box[0])
+
+
+@returns_box
+def add_margin(box: Box, margin) -> Box:
+    """
+    Returns a box with size increased by the ``margin`` (need to be broadcastable to the box)
+    compared to the input ``box``.
+    """
+    margin = np.broadcast_to(margin, box.shape)
+    return box[0] - margin[0], box[1] + margin[1]