diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml new file mode 100644 index 00000000..7b50bf31 --- /dev/null +++ b/.github/workflows/lint.yml @@ -0,0 +1,27 @@ +name: Lint + +on: [ pull_request ] + +env: + MODULE_NAME: imops + +jobs: + lint: + runs-on: ubuntu-22.04 + steps: + - uses: actions/checkout@v4 + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.12' + + - name: Check python code style + run: | + pip install -r requirements-linters.txt + flake8 . + isort --check . + black --check . + - name: Check Cython code style + run: | + pip install cython-lint + cython-lint imops/src diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index acd64184..cd1d9baa 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -35,7 +35,7 @@ jobs: needs: [ check_version ] strategy: matrix: - os: [ ubuntu-22.04, windows-2019, macOS-12 ] + os: [ ubuntu-22.04, windows-2019, macOS-13 ] name: Build wheels on ${{ matrix.os }} runs-on: ${{ matrix.os }} @@ -47,17 +47,10 @@ jobs: python-version: '3.9' - name: Install cibuildwheel run: python -m pip install cibuildwheel==2.17.0 - - name: Install gcc for mac - if: matrix.os == 'macOS-12' + - name: Install llvm for mac + if: matrix.os == 'macOS-13' run: | - brew install llvm libomp - echo $PATH - 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 + brew install llvm - name: Install g++-11 for ubuntu if: matrix.os == 'ubuntu-22.04' id: install_cc @@ -79,7 +72,7 @@ jobs: python -m cibuildwheel --output-dir wheelhouse env: CIBW_ENVIRONMENT_MACOS: > - PATH="/usr/local/opt/llvm/bin:$PATH" LDFLAGS="-L/usr/local/opt/llvm/lib" CPPFLAGS="-I/usr/local/opt/llvm/include" + CC="$(brew --prefix llvm)/bin/clang" CXX="$(brew --prefix llvm)/bin/clang++" CIBW_BUILD: cp37-* cp38-* cp39-* cp310-* cp311-* cp312-* CIBW_BEFORE_BUILD_LINUX: 'if [ $(python -c "import sys; print(sys.version_info[1])") -ge 9 ]; then python -m pip install "numpy<3.0.0" --config-settings=setup-args="-Dallow-noblas=true"; fi' - uses: actions/upload-artifact@v3 diff --git a/.github/workflows/test_build.yml b/.github/workflows/test_build.yml index 8f0d1718..d5022263 100644 --- a/.github/workflows/test_build.yml +++ b/.github/workflows/test_build.yml @@ -9,7 +9,7 @@ jobs: build_wheels: strategy: matrix: - os: [ubuntu-22.04, windows-2019, macOS-12 ] + os: [ubuntu-22.04, windows-2019, macOS-13 ] name: Build wheels on ${{ matrix.os }} runs-on: ${{ matrix.os }} @@ -21,17 +21,10 @@ jobs: python-version: '3.9' - name: Install cibuildwheel run: python -m pip install cibuildwheel==2.17.0 - - name: Install gcc for mac - if: matrix.os == 'macOS-12' + - name: Install llvm for mac + if: matrix.os == 'macOS-13' run: | - brew install llvm libomp - echo $PATH - 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 + brew install llvm - name: Install g++-11 for ubuntu if: matrix.os == 'ubuntu-22.04' id: install_cc @@ -53,7 +46,7 @@ jobs: python -m cibuildwheel --output-dir wheelhouse env: CIBW_ENVIRONMENT_MACOS: > - PATH="/usr/local/opt/llvm/bin:$PATH" LDFLAGS="-L/usr/local/opt/llvm/lib" CPPFLAGS="-I/usr/local/opt/llvm/include" + CC="$(brew --prefix llvm)/bin/clang" CXX="$(brew --prefix llvm)/bin/clang++" CIBW_BUILD: cp37-* cp39-* cp312-* CIBW_SKIP: "*manylinux_x86_64" CIBW_BEFORE_BUILD_LINUX: 'if [ $(python -c "import sys; print(sys.version_info[1])") -ge 9 ]; then python -m pip install "numpy<3.0.0" --config-settings=setup-args="-Dallow-noblas=true"; fi' diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 627ada08..634a26d5 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -43,17 +43,6 @@ jobs: echo $MODULE_PARENT echo "MODULE_PARENT=$(echo $MODULE_PARENT)" >> $GITHUB_ENV - - name: Check python code style - run: | - pip install -r requirements-linters.txt - flake8 . - isort --check . - black --check . - - name: Check Cython code style - run: | - pip install cython-lint - cython-lint imops/src - - name: Test with pytest run: | pytest tests -m "not nonumba" --junitxml=reports/junit-${{ matrix.python-version }}.xml --cov="$MODULE_PARENT/$MODULE_NAME" --cov-report=xml --cov-branch diff --git a/.gitignore b/.gitignore index 87e1733b..8d93aa0a 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,4 @@ imops/src/_fast*.pyx dist/ *.so .vscode/ +.idea/ diff --git a/MANIFEST.in b/MANIFEST.in index d25ac5dc..c77969d0 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,6 +1,6 @@ include *.md include requirements.txt include pyproject.toml -include _build_utils.py recursive-include imops *.py recursive-include imops/cpp *.h *.hpp *.cpp +exclude tests/* diff --git a/docs/index.md b/docs/index.md index 90ca24f5..3ffc0fae 100644 --- a/docs/index.md +++ b/docs/index.md @@ -48,6 +48,8 @@ pip install imops[numba] # additionally install Numba backend ::: imops.morphology.distance_transform_edt +::: imops.morphology.convex_hull_image + ::: imops.measure.label ::: imops.measure.center_of_mass diff --git a/imops/__version__.py b/imops/__version__.py index e4e49b3b..8969d496 100644 --- a/imops/__version__.py +++ b/imops/__version__.py @@ -1 +1 @@ -__version__ = '0.9.0' +__version__ = '0.9.1' diff --git a/_build_utils.py b/imops/_build_utils.py similarity index 98% rename from _build_utils.py rename to imops/_build_utils.py index c5f032ca..12c61668 100644 --- a/_build_utils.py +++ b/imops/_build_utils.py @@ -63,10 +63,10 @@ def get_ext_modules(): '/O3' if on_windows else '-O3', ] # FIXME: account for higher gcc versions - modules = ['backprojection', 'measure', 'morphology', 'numeric', 'radon', 'zoom'] + modules = ['backprojection', 'measure', 'morphology', 'numeric', 'radon', 'zoom', 'convex_hull'] modules_to_link_against_numpy_core_math_lib = ['numeric'] - src_dir = Path(__file__).parent / name / 'src' + src_dir = Path(__file__).parent / 'src' for module in modules: # 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 diff --git a/imops/compat.py b/imops/compat.py index 07a8a1df..5cfa7b8a 100644 --- a/imops/compat.py +++ b/imops/compat.py @@ -12,3 +12,13 @@ from scipy.ndimage._morphology import _ni_support except ImportError: from scipy.ndimage.morphology import _ni_support + +try: + from scipy.spatial import QhullError +except ImportError: + from scipy.spatial.qhull import QhullError + +from scipy.ndimage._nd_image import euclidean_feature_transform # noqa + + +normalize_sequence = _ni_support._normalize_sequence # noqa diff --git a/imops/crop.py b/imops/crop.py index d1e1a772..dbb2c331 100644 --- a/imops/crop.py +++ b/imops/crop.py @@ -1,3 +1,5 @@ +from typing import Optional + import numpy as np from .backend import BackendLike @@ -6,7 +8,9 @@ from .utils import AxesLike, AxesParams, assert_subdtype, broadcast_axis, fill_by_indices -def crop_to_shape(x: np.ndarray, shape: AxesLike, axis: AxesLike = None, ratio: AxesParams = 0.5) -> np.ndarray: +def crop_to_shape( + x: np.ndarray, shape: AxesLike, axis: Optional[AxesLike] = None, ratio: AxesParams = 0.5 +) -> np.ndarray: """ Crop `x` to match `shape` along `axis`. @@ -57,8 +61,8 @@ def crop_to_shape(x: np.ndarray, shape: AxesLike, axis: AxesLike = None, ratio: def crop_to_box( x: np.ndarray, box: np.ndarray, - axis: AxesLike = None, - padding_values: AxesParams = None, + axis: Optional[AxesLike] = None, + padding_values: Optional[AxesParams] = None, num_threads: int = _NUMERIC_DEFAULT_NUM_THREADS, backend: BackendLike = None, ) -> np.ndarray: diff --git a/imops/interp1d.py b/imops/interp1d.py index 269cf217..861a0b69 100644 --- a/imops/interp1d.py +++ b/imops/interp1d.py @@ -1,4 +1,4 @@ -from typing import Union +from typing import Optional, Union from warnings import warn import numpy as np @@ -71,7 +71,7 @@ def __init__( kind: Union[int, str] = 'linear', axis: int = -1, copy: bool = True, - bounds_error: bool = None, + bounds_error: Optional[bool] = None, fill_value: Union[float, str] = np.nan, assume_sorted: bool = False, num_threads: int = -1, diff --git a/imops/interp2d.py b/imops/interp2d.py index 5a3df754..306801c4 100644 --- a/imops/interp2d.py +++ b/imops/interp2d.py @@ -1,4 +1,5 @@ from platform import python_version +from typing import Optional import numpy as np from scipy.spatial import KDTree @@ -47,9 +48,9 @@ class Linear2DInterpolator(Linear2DInterpolatorCpp): def __init__( self, points: np.ndarray, - values: np.ndarray = None, + values: Optional[np.ndarray] = None, num_threads: int = 1, - triangles: np.ndarray = None, + triangles: Optional[np.ndarray] = None, **kwargs, ): if triangles is not None: @@ -77,7 +78,7 @@ def __init__( # FIXME: add backend dispatch self.num_threads = normalize_num_threads(num_threads, Cython(), warn_stacklevel=3) - def __call__(self, points: np.ndarray, values: np.ndarray = None, fill_value: float = 0.0) -> np.ndarray: + def __call__(self, points: np.ndarray, values: Optional[np.ndarray] = None, fill_value: float = 0.0) -> np.ndarray: """ Evaluate the interpolant diff --git a/imops/measure.py b/imops/measure.py index 056f2286..1a716826 100644 --- a/imops/measure.py +++ b/imops/measure.py @@ -1,6 +1,6 @@ from collections import namedtuple from platform import python_version -from typing import List, NamedTuple, Sequence, Tuple, Union +from typing import List, NamedTuple, Optional, Sequence, Tuple, Union from warnings import warn import numpy as np @@ -32,12 +32,12 @@ # TODO: Make it work and test on immutable arrays as soon as `cc3d` package is fixed def label( label_image: np.ndarray, - background: int = None, - connectivity: int = None, + background: Optional[int] = None, + connectivity: Optional[int] = None, return_num: bool = False, return_labels: bool = False, return_sizes: bool = False, - dtype: type = None, + dtype: Optional[type] = None, ) -> Union[np.ndarray, NamedTuple]: """ Fast version of `skimage.measure.label` which optionally returns number of connected components, labels and sizes. @@ -139,8 +139,8 @@ def label( def center_of_mass( array: np.ndarray, - labels: np.ndarray = None, - index: Union[int, Sequence[int]] = None, + labels: Union[np.ndarray, None] = None, + index: Union[int, Sequence[int], None] = None, num_threads: int = -1, backend: BackendLike = None, ) -> Union[Tuple[float, ...], List[Tuple[float, ...]]]: diff --git a/imops/morphology.py b/imops/morphology.py index 357d1243..e9e2bcc1 100644 --- a/imops/morphology.py +++ b/imops/morphology.py @@ -1,22 +1,24 @@ -from typing import Callable, Tuple, Union +from typing import Callable, Optional, Tuple, Union from warnings import warn import numpy as np from edt import edt from scipy.ndimage import distance_transform_edt as scipy_distance_transform_edt, generate_binary_structure -from scipy.ndimage._nd_image import euclidean_feature_transform +from scipy.spatial import ConvexHull from skimage.morphology import ( binary_closing as scipy_binary_closing, binary_dilation as scipy_binary_dilation, binary_erosion as scipy_binary_erosion, binary_opening as scipy_binary_opening, ) +from skimage.util import unique_rows from .backend import BackendLike, Cython, Scipy, resolve_backend from .box import add_margin, box_to_shape, mask_to_box, shape_to_box -from .compat import _ni_support +from .compat import QhullError, euclidean_feature_transform, normalize_sequence from .crop import crop_to_box from .pad import restore_crop +from .src._convex_hull import _grid_points_in_poly, _left_right_bounds, _offset_unique from .src._fast_morphology import ( _binary_dilation as cython_fast_binary_dilation, _binary_erosion as cython_fast_binary_erosion, @@ -30,8 +32,8 @@ def morphology_op_wrapper( ) -> Callable: def wrapped( image: np.ndarray, - footprint: np.ndarray = None, - output: np.ndarray = None, + footprint: Optional[np.ndarray] = None, + output: Optional[np.ndarray] = None, boxed: bool = False, num_threads: int = -1, backend: BackendLike = None, @@ -161,8 +163,8 @@ def wrapped( def binary_dilation( image: np.ndarray, - footprint: np.ndarray = None, - output: np.ndarray = None, + footprint: Optional[np.ndarray] = None, + output: Optional[np.ndarray] = None, boxed: bool = False, num_threads: int = -1, backend: BackendLike = None, @@ -215,8 +217,8 @@ def binary_dilation( def binary_erosion( image: np.ndarray, - footprint: np.ndarray = None, - output: np.ndarray = None, + footprint: Optional[np.ndarray] = None, + output: Optional[np.ndarray] = None, boxed: bool = False, num_threads: int = -1, backend: BackendLike = None, @@ -269,8 +271,8 @@ def binary_erosion( def binary_closing( image: np.ndarray, - footprint: np.ndarray = None, - output: np.ndarray = None, + footprint: Optional[np.ndarray] = None, + output: Optional[np.ndarray] = None, boxed: bool = False, num_threads: int = -1, backend: BackendLike = None, @@ -324,8 +326,8 @@ def binary_closing( def binary_opening( image: np.ndarray, - footprint: np.ndarray = None, - output: np.ndarray = None, + footprint: Optional[np.ndarray] = None, + output: Optional[np.ndarray] = None, boxed: bool = False, num_threads: int = -1, backend: BackendLike = None, @@ -369,7 +371,7 @@ def binary_opening( def distance_transform_edt( image: np.ndarray, - sampling: Tuple[float] = None, + sampling: Optional[Tuple[float]] = None, return_distances: bool = True, return_indices: bool = False, num_threads: int = -1, @@ -489,7 +491,7 @@ def distance_transform_edt( if image.dtype != bool: image = np.atleast_1d(np.where(image, 1, 0)) if sampling is not None: - sampling = _ni_support._normalize_sequence(sampling, image.ndim) + sampling = normalize_sequence(sampling, image.ndim) sampling = np.asarray(sampling, dtype=np.float64) if not sampling.flags.contiguous: sampling = sampling.copy() @@ -517,3 +519,73 @@ def distance_transform_edt( return result[0] return None + + +def convex_hull_image(image: np.ndarray, offset_coordinates: bool = True) -> np.ndarray: + """ + Fast convex hull of an image. Similar to skimage.morphology.convex_hull_image with include_borders=True + + Parameters + ---------- + image: np.ndarray + input image, must be 2D + offset_coordinates: bool + If True, a pixel at coordinate, e.g., (4, 7) will be represented by coordinates + (3.5, 7), (4.5, 7), (4, 6.5), and (4, 7.5). + This adds some “extent” to a pixel when computing the hull. + + Returns + ------- + output: np.ndarray + resulting convex hull of the input image + + Examples + -------- + ```python + chull = convex_hull_image(x) + ``` + """ + + ndim = image.ndim + if ndim != 2: + raise ValueError(f'convex_hull_image is currently implemented only for 2D arrays, got {ndim}D array') + + if np.count_nonzero(image) == 0: + return np.zeros(image.shape, dtype=bool) + + # In 2D, we do an optimisation by choosing only pixels that are + # the starting or ending pixel of a row or column. This vastly + # limits the number of coordinates to examine for the virtual hull. + image = np.ascontiguousarray(image, dtype=np.uint8) + + coords = _left_right_bounds(image) + + if offset_coordinates: + coords = _offset_unique(coords) + else: + coords = unique_rows(coords) + + # Find the convex hull + try: + hull = ConvexHull(coords) + except QhullError as err: + warn(f'Failed to get convex hull image. ' f'Returning empty image, see error message below: \n' f'{err}') + return np.zeros(image.shape, dtype=bool) + + vertices = hull.points[hull.vertices] + + # return vertices + + # If 2D, use fast Cython function to locate convex hull pixels + labels = _grid_points_in_poly( + np.ascontiguousarray(vertices[:, 0], dtype=np.float32), + np.ascontiguousarray(vertices[:, 1], dtype=np.float32), + image.shape[0], + image.shape[1], + len(vertices), + ) + + # edge points and vertices are included + mask = labels.view(bool) + + return mask diff --git a/imops/numeric.py b/imops/numeric.py index 8cca57f4..e9259589 100644 --- a/imops/numeric.py +++ b/imops/numeric.py @@ -1,4 +1,4 @@ -from typing import Callable, Sequence, Union +from typing import Callable, Optional, Sequence, Union import numpy as np @@ -99,7 +99,7 @@ def _choose_cython_copy(ndim: int, is_fp16: bool, fast: bool) -> Callable: def pointwise_add( nums: np.ndarray, summand: Union[np.array, int, float], - output: np.ndarray = None, + output: Optional[np.ndarray] = None, num_threads: int = _NUMERIC_DEFAULT_NUM_THREADS, backend: BackendLike = None, ) -> np.ndarray: @@ -256,7 +256,7 @@ def fill_( def full( shape: Union[int, Sequence[int]], fill_value: Union[np.number, int, float], - dtype: Union[type, str] = None, + dtype: Union[type, str, None] = None, order: str = 'C', num_threads: int = _NUMERIC_DEFAULT_NUM_THREADS, backend: BackendLike = None, @@ -302,7 +302,7 @@ def full( def copy( nums: np.ndarray, - output: np.ndarray = None, + output: Optional[np.ndarray] = None, order: str = 'K', num_threads: int = _NUMERIC_DEFAULT_NUM_THREADS, backend: BackendLike = None, diff --git a/imops/pad.py b/imops/pad.py index b277dcea..5aa0e64d 100644 --- a/imops/pad.py +++ b/imops/pad.py @@ -1,4 +1,4 @@ -from typing import Callable, Sequence, Union +from typing import Callable, Optional, Sequence, Union import numpy as np @@ -10,7 +10,7 @@ def pad( x: np.ndarray, padding: Union[AxesLike, Sequence[Sequence[int]]], - axis: AxesLike = None, + axis: Optional[AxesLike] = None, padding_values: Union[AxesParams, Callable] = 0, num_threads: int = _NUMERIC_DEFAULT_NUM_THREADS, backend: BackendLike = None, @@ -76,7 +76,7 @@ def pad( def pad_to_shape( x: np.ndarray, shape: AxesLike, - axis: AxesLike = None, + axis: Optional[AxesLike] = None, padding_values: Union[AxesParams, Callable] = 0, ratio: AxesParams = 0.5, num_threads: int = _NUMERIC_DEFAULT_NUM_THREADS, @@ -135,7 +135,7 @@ def pad_to_shape( def pad_to_divisible( x: np.ndarray, divisor: AxesLike, - axis: AxesLike = None, + axis: Optional[AxesLike] = None, padding_values: Union[AxesParams, Callable] = 0, ratio: AxesParams = 0.5, remainder: AxesLike = 0, diff --git a/imops/radon.py b/imops/radon.py index a3d122b5..5c068f58 100644 --- a/imops/radon.py +++ b/imops/radon.py @@ -1,4 +1,4 @@ -from typing import Sequence, Tuple, Union +from typing import Optional, Sequence, Tuple, Union import numpy as np from scipy.fftpack import fft, ifft @@ -15,7 +15,7 @@ def radon( image: np.ndarray, - axes: Tuple[int, int] = None, + axes: Optional[Tuple[int, int]] = None, theta: Union[int, Sequence[float]] = 180, return_fill: bool = False, num_threads: int = -1, @@ -104,8 +104,8 @@ def radon( def inverse_radon( sinogram: np.ndarray, - axes: Tuple[int, int] = None, - theta: Union[int, Sequence[float]] = None, + axes: Optional[Tuple[int, int]] = None, + theta: Union[int, Sequence[float], None] = None, fill_value: float = 0, a: float = 0, b: float = 1, diff --git a/imops/src/_convex_hull.pyx b/imops/src/_convex_hull.pyx new file mode 100644 index 00000000..351a70b1 --- /dev/null +++ b/imops/src/_convex_hull.pyx @@ -0,0 +1,262 @@ +# cython: cdivision=True +# cython: boundscheck=False +# cython: nonecheck=False +# cython: wraparound=False +import numpy as np + +cimport numpy as cnp +from libc.math cimport ceilf, floorf + + +cnp.import_array() + +FP_BOUND_DTYPE = np.dtype([ + ('lb', np.float32), + ('rb', np.float32), + ('assigned', np.uint8) +]) + +INT_BOUND_DTYPE = np.dtype([ + ('lb', np.int32), + ('rb', np.int32), + ('assigned', np.uint8) +]) + + +def _grid_points_in_poly(float[:] vx, float[:] vy, Py_ssize_t M, Py_ssize_t N, Py_ssize_t nr_verts): + cdef intBound tmp_int_bound + cdef Py_ssize_t m, n, i + cdef float prev_x, prev_y, curr_x, curr_y + cdef float tmp_from_x, tmp_from_y, tmp_to_x, tmp_to_y + cdef float lerp_t, bound_y + cdef Py_ssize_t x_set, x_start, x_stop + + cdef cnp.ndarray[dtype=cnp.uint8_t, ndim=2, mode="c"] out = np.zeros((M, N), dtype=np.uint8) + + cdef fpBound[:] fpBounds = np.empty(M, dtype=FP_BOUND_DTYPE) + cdef intBound[:] intBounds = np.empty(M, dtype=INT_BOUND_DTYPE) + + for i in range(M): + fpBounds[i].assigned = False + fpBounds[i].lb = float('inf') + fpBounds[i].rb = -1 + intBounds[i].assigned = False + + prev_x = vx[nr_verts - 1] + prev_y = vy[nr_verts - 1] + + # algorithm relies on vertex validity and counterclockwise orientation of the vertices + for i in range(nr_verts): + curr_x = vx[i] + curr_y = vy[i] + + if prev_x == curr_x: + x_set = (floorf(prev_x) if prev_y < curr_y else ceilf(prev_x)) + + fpBounds[x_set].assigned = True + fpBounds[x_set].lb = min(fpBounds[x_set].lb, prev_y, curr_y) + fpBounds[x_set].rb = max(fpBounds[x_set].rb, prev_y, curr_y) + else: + if prev_x < curr_x: + tmp_from_x = prev_x + tmp_from_y = prev_y + tmp_to_x = curr_x + tmp_to_y = curr_y + else: + tmp_from_x = curr_x + tmp_from_y = curr_y + tmp_to_x = prev_x + tmp_to_y = prev_y + + # vertices are treated as points on image, so include x_stop + x_start = ceilf(tmp_from_x) + x_stop = floorf(tmp_to_x + 1) + + for x_set in range(x_start, x_stop): + lerp_t = (x_set - tmp_from_x) / (tmp_to_x - tmp_from_x) + bound_y = lerp(tmp_from_y, tmp_to_y, lerp_t) + + fpBounds[x_set].assigned = True + fpBounds[x_set].lb = min(fpBounds[x_set].lb, bound_y) + fpBounds[x_set].rb = max(fpBounds[x_set].rb, bound_y) + + prev_x = curr_x + prev_y = curr_y + + # bounds are computed as point interpolation + # so bounds must be valid indices for out array + for m in range(M): + intBounds[m] = intify(fpBounds[m], 0, N - 1) + + for m in range(M): + tmp_int_bound = intBounds[m] + + if tmp_int_bound.assigned: + # Do not forget to fill right bound + for n in range(tmp_int_bound.lb, tmp_int_bound.rb + 1): + out[m, n] = True + + return out + + +# TODO: maybe use round instead of floorf and ceilf? +cdef inline intBound intify(fpBound bound, Py_ssize_t min_idx, Py_ssize_t max_idx): + if bound.assigned: + return intBound( + lb = max(min_idx, ceilf(bound.lb - 0.2)), + rb = min(max_idx, floorf(bound.rb + 0.2)), + assigned=True + ) + + return intBound(lb=0, rb=0, assigned=False) + + +cdef packed struct fpBound: + float lb + float rb + unsigned char assigned + + +cdef packed struct intBound: + int lb + int rb + unsigned char assigned + + +cdef inline float lerp(float y0, float y1, float t): + return y0 * (1 - t) + y1 * t + + +cpdef _left_right_bounds(cnp.uint8_t[:, :] image): + cdef Py_ssize_t i, j, M = image.shape[0], N = image.shape[1], curr_pos = 0, left, right + cdef cnp.ndarray[dtype=int, ndim=2, mode="c"] left_right_bounds = np.zeros((2 * M, 2), dtype=np.int32) + cdef unsigned char found = False + + for i in range(M): + found = False + + for j in range(N): + if image[i, j]: + left = j + found = True + break + + for j in range(N): + if image[i, N - 1 - j]: + right = N - 1 - j + found = True + break + + if found: + left_right_bounds[2 * curr_pos, 0] = i + left_right_bounds[2 * curr_pos, 1] = left + left_right_bounds[2 * curr_pos + 1, 0] = i + left_right_bounds[2 * curr_pos + 1, 1] = right + + curr_pos += 1 + + return np.ascontiguousarray(left_right_bounds[: 2 * curr_pos, :]) + + +cdef inline int set_unique_curr(float* expanded_bounds, int x, int l_b, int r_b): + if l_b == r_b: + expanded_bounds[0] = x + expanded_bounds[1] = l_b - 0.5 + + expanded_bounds[2] = x - 0.5 + expanded_bounds[3] = l_b + + expanded_bounds[4] = x + expanded_bounds[5] = l_b + 0.5 + + return 3 + elif r_b == l_b + 1: + expanded_bounds[0] = x + expanded_bounds[1] = l_b - 0.5 + + expanded_bounds[2] = x - 0.5 + expanded_bounds[3] = l_b + + expanded_bounds[4] = x + expanded_bounds[5] = l_b + 0.5 + + expanded_bounds[6] = x - 0.5 + expanded_bounds[7] = r_b + + expanded_bounds[8] = x + expanded_bounds[9] = r_b + 0.5 + + return 5 + + else: + expanded_bounds[0] = x + expanded_bounds[1] = l_b - 0.5 + + expanded_bounds[2] = x - 0.5 + expanded_bounds[3] = l_b + + expanded_bounds[4] = x + expanded_bounds[5] = l_b + 0.5 + + expanded_bounds[6] = x + expanded_bounds[7] = r_b - 0.5 + + expanded_bounds[8] = x - 0.5 + expanded_bounds[9] = r_b + + expanded_bounds[10] = x + expanded_bounds[11] = r_b + 0.5 + + return 6 + + +cpdef _offset_unique(int[:, :] left_right_bounds): + cdef Py_ssize_t N = left_right_bounds.shape[0], i, curr_pos = 0 + cdef cnp.ndarray[dtype=float, ndim=2, mode="c"] expanded_bounds = np.zeros((4 * N, 2), dtype=np.float32) + + cdef int x_l_prev, y_l_prev, x_r_prev, y_r_prev, x_l_curr, y_l_curr, x_r_curr, y_r_curr + + x_l_prev = left_right_bounds[0, 0] + y_l_prev = left_right_bounds[0, 1] + x_r_prev = left_right_bounds[1, 0] + y_r_prev = left_right_bounds[1, 1] + + curr_pos += set_unique_curr(&expanded_bounds[0, 0], x_l_prev, y_l_prev, y_r_prev) + + for i in range(1, N // 2): + x_l_curr = left_right_bounds[2 * i, 0] + y_l_curr = left_right_bounds[2 * i, 1] + x_r_curr = left_right_bounds[2 * i + 1, 0] + y_r_curr = left_right_bounds[2 * i + 1, 1] + + curr_pos += set_unique_curr(&expanded_bounds[curr_pos, 0], x_l_curr, y_l_curr, y_r_curr) + + if x_l_prev + 1 == x_l_curr and (y_l_prev == y_l_curr or y_l_prev == y_r_curr): + pass + else: + expanded_bounds[curr_pos, 0] = x_l_prev + 0.5 + expanded_bounds[curr_pos, 1] = y_l_prev + curr_pos += 1 + + if x_l_prev + 1 == x_l_curr and (y_r_prev == y_l_curr or y_r_prev == y_r_curr) and (y_r_prev != y_l_prev): + pass + else: + expanded_bounds[curr_pos, 0] = x_l_prev + 0.5 + expanded_bounds[curr_pos, 1] = y_r_prev + curr_pos += 1 + + x_l_prev = x_l_curr + y_l_prev = y_l_curr + x_r_prev = x_r_curr + y_r_prev = y_r_curr + + expanded_bounds[curr_pos, 0] = x_l_prev + 0.5 + expanded_bounds[curr_pos, 1] = y_l_prev + curr_pos += 1 + + if y_r_prev != y_l_prev: + expanded_bounds[curr_pos, 0] = x_l_prev + 0.5 + expanded_bounds[curr_pos, 1] = y_r_prev + curr_pos += 1 + + return np.ascontiguousarray(expanded_bounds[:curr_pos, :]) diff --git a/imops/utils.py b/imops/utils.py index 4fe7ad3e..4b1bd971 100644 --- a/imops/utils.py +++ b/imops/utils.py @@ -1,4 +1,5 @@ import os +import platform from contextlib import contextmanager from itertools import permutations from typing import Callable, Optional, Sequence, Tuple, Union @@ -53,7 +54,7 @@ def normalize_num_threads(num_threads: int, backend: Backend, warn_stacklevel: i env_num_threads = os.environ.get(env_num_threads_var_name, '').strip() env_num_threads = int(env_num_threads) if env_num_threads else None # TODO: maybe let user set the absolute maximum number of threads? - num_available_cpus = len(os.sched_getaffinity(0)) + num_available_cpus = os.cpu_count() if platform.system() == 'Darwin' else len(os.sched_getaffinity(0)) max_num_threads = min(filter(bool, [IMOPS_NUM_THREADS, env_num_threads, num_available_cpus])) @@ -168,7 +169,9 @@ def wrapper( return wrapper -def build_slices(start: Sequence[int], stop: Sequence[int] = None, step: Sequence[int] = None) -> Tuple[slice, ...]: +def build_slices( + start: Sequence[int], stop: Optional[Sequence[int]] = None, step: Optional[Sequence[int]] = None +) -> Tuple[slice, ...]: """ Returns a tuple of slices built from `start` and `stop` with `step`. diff --git a/imops/zoom.py b/imops/zoom.py index a3ceb9f0..b5a0bdeb 100644 --- a/imops/zoom.py +++ b/imops/zoom.py @@ -1,5 +1,5 @@ from platform import python_version -from typing import Callable, Sequence, Union +from typing import Callable, Optional, Sequence, Union from warnings import warn import numpy as np @@ -72,7 +72,7 @@ def _choose_numba_zoom(ndim: int, order: int) -> Callable: def zoom( x: np.ndarray, scale_factor: AxesParams, - axis: AxesLike = None, + axis: Optional[AxesLike] = None, order: int = 1, fill_value: Union[float, Callable] = 0, num_threads: int = -1, @@ -129,7 +129,7 @@ def zoom( def zoom_to_shape( x: np.ndarray, shape: AxesLike, - axis: AxesLike = None, + axis: Optional[AxesLike] = None, order: int = 1, fill_value: Union[float, Callable] = 0, num_threads: int = -1, @@ -191,7 +191,7 @@ def zoom_to_shape( def _zoom( image: np.ndarray, zoom: Sequence[float], - output: np.ndarray = None, + output: Optional[np.ndarray] = None, order: int = 1, mode: str = 'constant', cval: float = 0.0, diff --git a/pyproject.toml b/pyproject.toml index eed6e582..2f7a669b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,17 +4,17 @@ build-backend = 'setuptools.build_meta' [project] name = 'imops' -dynamic = ['version'] +dynamic = ['dependencies', 'version'] description = 'Efficient parallelizable algorithms for multidimensional arrays to speed up your data pipelines' readme = 'README.md' requires-python = '>=3.7' license = { file = 'LICENSE' } keywords = ['image processing', 'fast', 'ndarray', 'data pipelines'] authors = [ - {name = 'maxme1', email = 'max@aumi.ai'}, - {name = 'vovaf709', email = 'vovaf709@yandex.ru'}, - {name = 'talgat', email = 'saparov2130@gmail.com'}, - {name = 'alexeybelkov', email='fpmbelkov@gmail.com'} + { name = 'maxme1', email = 'max@aumi.ai' }, + { name = 'vovaf709', email = 'vovaf709@yandex.ru' }, + { name = 'talgat', email = 'saparov2130@gmail.com' }, + { name = 'alexeybelkov', email = 'fpmbelkov@gmail.com' } ] classifiers = [ 'Development Status :: 5 - Production/Stable', @@ -51,23 +51,21 @@ line_length = 120 lines_after_imports = 2 profile = 'black' combine_as_imports = true -skip_glob=['.asv/*', '.eggs/*'] +skip_glob = ['.asv/*', '.eggs/*'] [tool.cython-lint] max-line-length = 120 -[tool.setuptools] -py-modules = ['_build_utils'] - [tool.setuptools.cmdclass] -build_py = "_build_utils.PyprojectBuild" +build_py = "imops._build_utils.PyprojectBuild" [tool.setuptools.packages.find] include = ['imops'] +exclude = ['tests'] [tool.setuptools.package-data] imops = ['py.typed'] [tool.setuptools.dynamic] -version = {attr = 'imops.__version__.__version__'} +version = { attr = 'imops.__version__.__version__' } dependencies = { file = ['requirements.txt'] } diff --git a/setup.py b/setup.py index e6e3c50e..a2fee0c8 100644 --- a/setup.py +++ b/setup.py @@ -23,13 +23,14 @@ long_description = file.read() version = runpy.run_path(root / name / '__version__.py')['__version__'] -scope = {'__file__': __file__} -exec((root / '_build_utils.py').read_text(), scope) +build_utils = root / name / '_build_utils.py' +scope = {'__file__': str(build_utils)} +exec(build_utils.read_text(), scope) ext_modules = scope['get_ext_modules']() setup( name=name, - packages=find_packages(include=(name,)), + packages=find_packages(include=(name,), exclude=('tests', 'tests.*')), include_package_data=True, version=version, description='Efficient parallelizable algorithms for multidimensional arrays to speed up your data pipelines', @@ -51,8 +52,8 @@ extras_require={'numba': ['numba'], 'all': ['numba']}, setup_requires=[ 'setuptools<69.0.0', - 'Cython>=3.0.0,<4.0.0', 'numpy<3.0.0', + 'Cython>=3.0.0,<4.0.0', 'pybind11', ], ext_modules=ext_modules, diff --git a/tests/test_convex_hull.py b/tests/test_convex_hull.py new file mode 100644 index 00000000..8d83273a --- /dev/null +++ b/tests/test_convex_hull.py @@ -0,0 +1,113 @@ +import numpy as np +import pytest +from skimage import data +from skimage.morphology import convex_hull_image +from skimage.morphology.convex_hull import _offsets_diamond +from skimage.util import invert, unique_rows + +from imops.morphology import convex_hull_image as convex_hull_image_fast +from imops.src._convex_hull import _left_right_bounds, _offset_unique + + +np.random.seed(1337) +N_STRESS = 1000 + + +@pytest.fixture(params=[False, True]) +def offset_coordinates(request): + return request.param + + +def test_bounds(): + for _ in range(N_STRESS): + image = np.zeros((100, 100), dtype=bool) + image[20:70, 20:90] = np.random.randn(50, 70) > 1.5 + + im_any = np.any(image, axis=1) + x_indices = np.arange(0, image.shape[0])[im_any] + y_indices_left = np.argmax(image[im_any], axis=1) + y_indices_right = image.shape[1] - 1 - np.argmax(image[im_any][:, ::-1], axis=1) + left = np.stack((x_indices, y_indices_left), axis=-1) + right = np.stack((x_indices, y_indices_right), axis=-1) + coords_ref = np.vstack((left, right)) + + coords = _left_right_bounds(image) + + # _left_right_bounds has another order + assert len(unique_rows(coords_ref)) == len(unique_rows(np.concatenate((coords, coords_ref), 0))) + + +def test_offset(): + for _ in range(N_STRESS): + image = np.zeros((100, 100), dtype=bool) + image[20:70, 20:90] = np.random.randn(50, 70) > 1.5 + + coords = _left_right_bounds(image) + + offsets = _offsets_diamond(2) + coords_ref = unique_rows((coords[:, None, :] + offsets).reshape(-1, 2)) + + coords = _offset_unique(coords) + + # _left_right_bounds has another order + assert len(coords_ref) == len(unique_rows(np.concatenate((coords, coords_ref), 0))) + + +def test_convex_hull_image(offset_coordinates): + image = invert(data.horse()) + + try: + chull_ref = convex_hull_image(image, offset_coordinates=offset_coordinates, include_borders=True) + except TypeError: + chull_ref = convex_hull_image(image, offset_coordinates=offset_coordinates) + + chull = convex_hull_image_fast(image, offset_coordinates=offset_coordinates) + + assert (chull >= image).all() + assert (chull >= chull_ref).all() + + assert ((chull > chull_ref).sum() / chull_ref.sum()) < 5e-3 + + +def test_convex_hull_image_random(offset_coordinates): + for _ in range(N_STRESS): + image = np.zeros((200, 200), dtype=bool) + + image[15:-15, 5:-25] = np.random.randn(170, 170) > 3 + + try: + chull_ref = convex_hull_image(image, offset_coordinates=offset_coordinates, include_borders=True) + except TypeError: + chull_ref = convex_hull_image(image, offset_coordinates=offset_coordinates) + + chull = convex_hull_image_fast(image, offset_coordinates=offset_coordinates) + + assert (chull >= image).all() + assert (chull >= chull_ref).all() + + assert ((chull > chull_ref).sum() / chull_ref.sum()) < 5e-3 + + +def test_convex_hull_image_non2d(offset_coordinates): + image = np.zeros((3, 3, 3), dtype=bool) + + with pytest.raises(ValueError): + _ = convex_hull_image_fast(image, offset_coordinates=offset_coordinates) + + +def test_convex_hull_image_empty(offset_coordinates): + image = np.zeros((10, 10), dtype=bool) + chull = convex_hull_image_fast(image, offset_coordinates=offset_coordinates) + + assert (chull == np.zeros_like(chull)).all() + + +def test_convex_hull_image_qhullsrc_issues(): + image = np.zeros((10, 10), dtype=bool) + image[1, 1] = True + image[-2, -2] = True + + with pytest.warns(UserWarning): + chull = convex_hull_image_fast(image, offset_coordinates=False) + + assert (chull == np.zeros_like(chull)).all()