diff options
Diffstat (limited to 'silx/image')
27 files changed, 0 insertions, 4640 deletions
diff --git a/silx/image/__init__.py b/silx/image/__init__.py deleted file mode 100644 index 12bf320..0000000 --- a/silx/image/__init__.py +++ /dev/null @@ -1,33 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2017-2018 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ -"""This package provides some processing functions for 2D images. - -For more processing functions, see the silx.math and silx.opencl packages. - -See silx documentation: http://www.silx.org/doc/silx/latest/ -""" - -__authors__ = ["H. Payno"] -__license__ = "MIT" -__date__ = "14/02/2017" diff --git a/silx/image/_boundingbox.py b/silx/image/_boundingbox.py deleted file mode 100644 index 1c086b1..0000000 --- a/silx/image/_boundingbox.py +++ /dev/null @@ -1,100 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# -# Copyright (c) 2018 European Synchrotron Radiation Facility -# -# 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. -# -# ###########################################################################*/ -"""offer some generic 2D bounding box features -""" - -__authors__ = ["H. Payno"] -__license__ = "MIT" -__date__ = "07/09/2019" - -from silx.math.combo import min_max -import numpy - - -class _BoundingBox: - """ - Simple 2D bounding box - - :param tuple bottom_left: (y, x) bottom left point - :param tuple top_right: (y, x) top right point - """ - def __init__(self, bottom_left, top_right): - self.bottom_left = bottom_left - self.top_right = top_right - self.min_x = bottom_left[1] - self.min_y = bottom_left[0] - self.max_x = top_right[1] - self.max_y = top_right[0] - - def contains(self, item): - """ - - :param item: bouding box or point. If it is bounding check check if the - bottom_left and top_right corner of the given bounding box - are inside the current bounding box - :type: Union[BoundingBox, tuple] - :return: - """ - if isinstance(item, _BoundingBox): - return self.contains(item.bottom_left) and self.contains(item.top_right) - else: - return ( - (self.min_x <= item[1] <= self.max_x) and - (self.min_y <= item[0] <= self.max_y) - ) - - def collide(self, bb): - """ - Check if the two bounding box collide - - :param bb: bounding box to compare with - :type: :class:BoundingBox - :return: True if the two boxes collides - :rtype: bool - """ - assert isinstance(bb, _BoundingBox) - return ( - (self.min_x < bb.max_x and self.max_x > bb.min_x) and - (self.min_y < bb.max_y and self.max_y > bb.min_y) - ) - - @staticmethod - def from_points(points): - """ - - :param numpy.array tuple points: list of points. Should be 2D: - [(y1, x1), (y2, x2), (y3, x3), ...] - :return: bounding box from two points - :rtype: _BoundingBox - """ - if not isinstance(points, numpy.ndarray): - points_ = numpy.ndarray(points) - else: - points_ = points - x = points_[:, 1] - y = points_[:, 0] - x_min, x_max = min_max(x) - y_min, y_max = min_max(y) - return _BoundingBox(bottom_left=(y_min, x_min), top_right=(y_max, x_max)) diff --git a/silx/image/backprojection.py b/silx/image/backprojection.py deleted file mode 100644 index 63f99ca..0000000 --- a/silx/image/backprojection.py +++ /dev/null @@ -1,25 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2017 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ - -from silx.opencl.backprojection import * diff --git a/silx/image/bilinear.pyx b/silx/image/bilinear.pyx deleted file mode 100644 index 14547f8..0000000 --- a/silx/image/bilinear.pyx +++ /dev/null @@ -1,465 +0,0 @@ -# -*- coding: utf-8 -*- -#cython: embedsignature=True, language_level=3 -## This is for optimisation -#cython: boundscheck=False, wraparound=False, cdivision=True, initializedcheck=False, -## This is for developping: -##cython: profile=True, warn.undeclared=True, warn.unused=True, warn.unused_result=False, warn.unused_arg=True -# -# -# -# Project: silx (originally pyFAI) -# https://github.com/silx-kit/silx -# -# Copyright (C) 2012-2020 European Synchrotron Radiation Facility, Grenoble, France -# -# 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. - -"""Bilinear interpolator, peak finder, line-profile for images""" -__authors__ = ["J. Kieffer"] -__license__ = "MIT" -__date__ = "26/11/2020" - -# C-level imports -from libc.stdint cimport uint8_t -from libc.math cimport floor, ceil, sqrt, NAN, isfinite -from libc.float cimport FLT_MAX - -import cython -import numpy -import logging -logger = logging.getLogger(__name__) - - -#Definition of some constants -# How data are stored -ctypedef float data_t -data_d = numpy.float32 - -#How the mask is stored -ctypedef uint8_t mask_t -mask_d = numpy.uint8 - - -cdef class BilinearImage: - """Bilinear interpolator for images ... or any data on a regular grid - """ - cdef: - readonly data_t[:, ::1] data - readonly mask_t[:, ::1] mask - readonly data_t maxi, mini - readonly Py_ssize_t width, height - readonly bint has_mask - - # C-level declarations - cpdef Py_ssize_t coarse_local_maxi(self, Py_ssize_t) - cdef Py_ssize_t c_local_maxi(self, Py_ssize_t) nogil - cdef data_t c_funct(self, data_t, data_t) nogil - cdef void _init_min_max(self) nogil - - def __cinit__(self, data not None, mask=None): - """Constructor - - :param data: image as a 2D array - """ - assert data.ndim == 2 - self.height = data.shape[0] - self.width = data.shape[1] - self.data = numpy.ascontiguousarray(data, dtype=data_d) - if mask is not None: - self.mask = numpy.ascontiguousarray(mask, dtype=mask_d) - self.has_mask=True - else: - self.mask = None - self.has_mask = False - self._init_min_max() - - def __dealloc__(self): - self.data = None - self.mask = None - - def __call__(self, coord): - """Function f((y, x)) where f is a continuous function - made from the image and (y,x)=(row, column) is the pixel coordinates - in natural C-order - - :param x: 2-tuple of float (row, column) - :return: Interpolated signal from the image - """ - return self.c_funct(coord[1], coord[0]) - - cdef void _init_min_max(self) nogil: - "Calculate the min & max" - cdef: - Py_ssize_t i, j - data_t maxi, mini, value - mini = FLT_MAX - maxi = -FLT_MAX - for i in range(self.height): - for j in range(self.width): - if not (self.has_mask and self.mask[i,j]): - value = self.data[i, j] - maxi = max(value, maxi) - mini = min(value, mini) - self.maxi = maxi - self.mini = mini - - cdef data_t c_funct(self, data_t x, data_t y) nogil: - """Function f(x, y) where f is a continuous function - made from the image. - - :param x (float): column coordinate - :param y (float): row coordinate - :return: Interpolated signal from the image (nearest for outside) - - Cython only function due to NOGIL - """ - cdef: - data_t d0 = min(max(y, 0.0), (self.height - 1.0)) - data_t d1 = min(max(x, 0.0), (self.width - 1.0)) - mask_t m0, m1, m2, m3 - Py_ssize_t i0, i1, j0, j1 - data_t x0, x1, y0, y1, res, scale - - x0 = floor(d0) - x1 = ceil(d0) - y0 = floor(d1) - y1 = ceil(d1) - i0 = < int > x0 - i1 = < int > x1 - j0 = < int > y0 - j1 = < int > y1 - if (i0 == i1) and (j0 == j1): - if not (self.has_mask and self.mask[i0,j0]): - res = self.data[i0, j0] - else: - res = NAN - elif i0 == i1: - if self.has_mask: - m0 = self.mask[i0, j0] - m1 = self.mask[i0, j1] - if m0 and m1: - res = NAN - elif m0: - res = self.data[i0, j1] - elif m1: - res = self.data[i0, j0] - else: - res = (self.data[i0, j0] * (y1 - d1)) + (self.data[i0, j1] * (d1 - y0)) - else: - res = (self.data[i0, j0] * (y1 - d1)) + (self.data[i0, j1] * (d1 - y0)) - elif j0 == j1: - if self.has_mask: - m0 = self.mask[i0, j0] - m1 = self.mask[i1, j0] - if m0 and m1: - res = NAN - elif m0: - res = self.data[i1, j0] - elif m1: - res = self.data[i0, j0] - else: - res = (self.data[i0, j0] * (x1 - d0)) + (self.data[i1, j0] * (d0 - x0)) - else: - res = (self.data[i0, j0] * (x1 - d0)) + (self.data[i1, j0] * (d0 - x0)) - else: - if self.has_mask: - m0 = self.mask[i0, j0] - m1 = self.mask[i1, j0] - m2 = self.mask[i0, j1] - m3 = self.mask[i1, j1] - if m0 and m1 and m2 and m3: - res = NAN - else: - m0 = not m0 - m1 = not m1 - m2 = not m2 - m3 = not m3 - if m0 and m1 and m2 and m3: - res = (self.data[i0, j0] * (x1 - d0) * (y1 - d1)) \ - + (self.data[i1, j0] * (d0 - x0) * (y1 - d1)) \ - + (self.data[i0, j1] * (x1 - d0) * (d1 - y0)) \ - + (self.data[i1, j1] * (d0 - x0) * (d1 - y0)) - else: - res = (m0 * self.data[i0, j0] * (x1 - d0) * (y1 - d1)) \ - + (m1 * self.data[i1, j0] * (d0 - x0) * (y1 - d1)) \ - + (m2 * self.data[i0, j1] * (x1 - d0) * (d1 - y0)) \ - + (m3 * self.data[i1, j1] * (d0 - x0) * (d1 - y0)) - scale = ((m0 * (x1 - d0) * (y1 - d1)) - + (m1 * (d0 - x0) * (y1 - d1)) - + (m2 * (x1 - d0) * (d1 - y0)) - + (m3 * (d0 - x0) * (d1 - y0))) - res /= scale - else: - res = (self.data[i0, j0] * (x1 - d0) * (y1 - d1)) \ - + (self.data[i1, j0] * (d0 - x0) * (y1 - d1)) \ - + (self.data[i0, j1] * (x1 - d0) * (d1 - y0)) \ - + (self.data[i1, j1] * (d0 - x0) * (d1 - y0)) - - return res - - def opp_f(self, coord): - """Function -f((y,x)) for peak finding via minimizer. - - Gives large number outside the boundaries to return into the image - - :param x: 2-tuple of float in natural C order, i.e (row, column) - :return: Negative interpolated signal from the image - """ - cdef: - data_t d0, d1, res - d0, d1 = coord - if d0 < 0: - res = self.mini + d0 - elif d1 < 0: - res = self.mini + d1 - elif d0 > (self.height - 1): - res = self.mini - d0 + self.height - 1 - elif d1 > self.width - 1: - res = self.mini - d1 + self.width - 1 - else: - res = self.c_funct(d1, d0) - return - res - - def local_maxi(self, coord): - """Return the nearest local maximum ... with sub-pixel refinement - - Nearest maximum search: - steepest ascent - - Sub-pixel refinement: - Second order Taylor expansion of the function; - At the maximum, the first derivative is null - delta = x-i = -Inverse[Hessian].gradient - if Hessian is singular or \|delta\|>1: use a center of mass. - - :param coord: 2-tuple of scalar (row, column) - :return: 2-tuple of float with the nearest local maximum - """ - cdef: - int res, current0, current1 - int i0, i1 - data_t tmp, sum0 = 0, sum1 = 0, sum = 0 - data_t a00, a01, a02, a10, a11, a12, a20, a21, a22 - data_t d00, d11, d01, denom, delta0, delta1 - res = self.c_local_maxi(round(coord[0]) * self.width + round(coord[1])) - - current0 = res // self.width - current1 = res % self.width - if (current0 > 0) and (current0 < self.height - 1) and (current1 > 0) and (current1 < self.width - 1): - # Use second order polynomial Taylor expansion - a00 = self.data[current0 - 1, current1 - 1] - a01 = self.data[current0 - 1, current1 ] - a02 = self.data[current0 - 1, current1 + 1] - a10 = self.data[current0 , current1 - 1] - a11 = self.data[current0 , current1 ] - a12 = self.data[current0 , current1 + 1] - a20 = self.data[current0 + 1, current1 - 1] - a21 = self.data[current0 + 1, current1 ] - a22 = self.data[current0 + 1, current1 - 1] - d00 = a12 - 2.0 * a11 + a10 - d11 = a21 - 2.0 * a11 + a01 - d01 = (a00 - a02 - a20 + a22) / 4.0 - denom = 2.0 * (d00 * d11 - d01 * d01) - if abs(denom) < 1e-10: - logger.debug("Singular determinant, Hessian undefined") - else: - delta0 = ((a12 - a10) * d01 + (a01 - a21) * d11) / denom - delta1 = ((a10 - a12) * d00 + (a21 - a01) * d01) / denom - if abs(delta0) <= 1.0 and abs(delta1) <= 1.0: - # Result is OK if lower than 0.5. - return (delta0 + float(current0), delta1 + float(current1)) - else: - logger.debug("Failed to find root using second order expansion") - # refinement of the position by a simple center of mass of the last valid region used - for i0 in range(current0 - 1, current0 + 2): - for i1 in range(current1 - 1, current1 + 2): - tmp = self.data[i0, i1] - sum0 += tmp * i0 - sum1 += tmp * i1 - sum += tmp - if sum > 0: - return (sum0 / sum, sum1 / sum) - - return (float(current0), float(current1)) - - cpdef Py_ssize_t coarse_local_maxi(self, Py_ssize_t x): - """Return the nearest local maximum ... without sub-pixel refinement - - :param idx: start index (=row*width+column) - :return: local maximum index - """ - return self.c_local_maxi(x) - - cdef Py_ssize_t c_local_maxi(self, Py_ssize_t idx) nogil: - """Return the nearest local maximum without sub-pixel refinement - - :param idx: start index (=row*width+column) - :return: local maximum index - - This method is Cython only due to the NOGIL - """ - cdef: - Py_ssize_t current0 = idx // self.width - Py_ssize_t current1 = idx % self.width - Py_ssize_t i0, i1, start0, stop0, start1, stop1, new0, new1, rng, cnt - mask_t m - data_t tmp, value, old_value - - if self.has_mask and self.mask[current0, current1]: - #Start searching for a non masked pixel. - rng = 0 - cnt = 0 - value = self.mini - new0, new1 = current0, current1 - while cnt == 0: - rng += 1 - cnt = 0 - start0 = max(0, current0 - rng) - stop0 = min(self.height, current0 + rng + 1) - start1 = max(0, current1 - rng) - stop1 = min(self.width, current1 + rng + 1) - for i0 in range(start0, stop0): - for i1 in range(start1, stop1): - m = not self.mask[i0, i1] - cnt += m - if m: - tmp = self.data[i0, i1] - if tmp > value: - new0, new1 = i0, i1 - value = tmp - current0, current1 = new0, new1 - else: - value = self.data[current0, current1] - - old_value = value -1 - new0, new1 = current0, current1 - - while value > old_value: - old_value = value - start0 = max(0, current0 - 1) - stop0 = min(self.height, current0 + 2) - start1 = max(0, current1 - 1) - stop1 = min(self.width, current1 + 2) - for i0 in range(start0, stop0): - for i1 in range(start1, stop1): - if self.has_mask and self.mask[current0, current1]: - continue - tmp = self.data[i0, i1] - if tmp > value: - new0, new1 = i0, i1 - value = tmp - current0, current1 = new0, new1 - return self.width * current0 + current1 - - def map_coordinates(self, coordinates): - """Map coordinates of the array on the image - - :param coordinates: 2-tuple of array of the same size (row_array, column_array) - :return: array of values at given coordinates - """ - cdef: - data_t[:] d0, d1, res - Py_ssize_t size, i - shape = coordinates[0].shape - size = coordinates[0].size - d0 = numpy.ascontiguousarray(coordinates[0].ravel(), dtype=data_d) - d1 = numpy.ascontiguousarray(coordinates[1].ravel(), dtype=data_d) - assert size == d1.size - res = numpy.empty(size, dtype=data_d) - with nogil: - for i in range(size): - res[i] = self.c_funct(d1[i], d0[i]) - return numpy.asarray(res).reshape(shape) - - def profile_line(self, src, dst, int linewidth=1, method='mean'): - """Return the mean or sum of intensity profile of an image measured - along a scan line. - - :param src: The start point of the scan line. - :type src: 2-tuple of numeric scalar - :param dst: The end point of the scan line. - The destination point is included in the profile, - in contrast to standard numpy indexing. - :type dst: 2-tuple of numeric scalar - :param int linewidth: Width of the scanline (unit image pixel). - :param str method: 'mean' or 'sum' depending if we want to compute the - mean intensity along the line or the sum. - :return: The intensity profile along the scan line. - The length of the profile is the ceil of the computed length - of the scan line. - :rtype: 1d array - - Inspired from skimage - """ - cdef: - data_t src_row, src_col, dst_row, dst_col, d_row, d_col - data_t length, col_width, row_width, sum, row, col, new_row, new_col, val - Py_ssize_t lengt, i, j, cnt - bint compute_mean - data_t[::1] result - src_row, src_col = src - dst_row, dst_col = dst - if (src_row == dst_row) and (src_col == dst_col): - logger.warning("Source and destination points are the same") - return numpy.array([self.c_funct(src_col, src_row)]) - d_row = dst_row - src_row - d_col = dst_col - src_col - - # Offsets to deal with linewidth - length = sqrt(d_row * d_row + d_col * d_col) - row_width = d_col / length - col_width = - d_row / length - - lengt = <int> ceil(length + 1) - d_row /= <data_t> (lengt -1) - d_col /= <data_t> (lengt -1) - - result = numpy.zeros(lengt, dtype=data_d) - - # Offset position to the center of the bottom pixels of the profile - src_row -= row_width * (linewidth - 1) / 2. - src_col -= col_width * (linewidth - 1) / 2. - - compute_mean = (method == 'mean') - with nogil: - for i in range(lengt): - sum = 0 - cnt = 0 - - row = src_row + i * d_row - col = src_col + i * d_col - - for j in range(linewidth): - new_row = row + j * row_width - new_col = col + j * col_width - if ((new_col >= 0) and (new_col < self.width) and - (new_row >= 0) and (new_row < self.height)): - val = self.c_funct(new_col, new_row) - if isfinite(val): - cnt += 1 - sum += val - if cnt: - if compute_mean: - result[i] += sum / cnt - else: - result[i] += sum - elif compute_mean: - result[i] += NAN - # Ensures the result is exported as numpy array and not memory view. - return numpy.asarray(result) diff --git a/silx/image/marchingsquares/__init__.py b/silx/image/marchingsquares/__init__.py deleted file mode 100644 index a47a7f6..0000000 --- a/silx/image/marchingsquares/__init__.py +++ /dev/null @@ -1,117 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2018 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ -""" -This module provides implementations based on marching squares algorithms. - -The main implementation is done by :class:`MarchingSquaresMergeImpl`. It was -designed to speed up the computation of iso surface using Cython and OpenMP. -It also provides features like support of mask, and cache of min/max per tiles -which is very efficient to find many iso contours from image gradient. - -Utilitary functions are provided as facade for simple use. -:meth:`find_contours` to find iso contours from an image and using the same -main signature as `find_contours` from `skimage`, but supporting mask. -And :meth:`find_pixels` which returns a set of pixel coords containing the -points of the iso contours. -""" - -__authors__ = ["V. Valls"] -__license__ = "MIT" -__date__ = "02/07/2018" - - -from ._mergeimpl import MarchingSquaresMergeImpl - - -def _factory(engine, image, mask): - """Factory to create the marching square implementation from the engine - name""" - if engine == "merge": - return MarchingSquaresMergeImpl(image, mask) - elif engine == "skimage": - from _skimage import MarchingSquaresSciKitImage - return MarchingSquaresSciKitImage(image, mask) - else: - raise ValueError("Engine '%s' is not supported ('merge' or 'skimage' expected).") - - -def find_pixels(image, level, mask=None): - """ - Find the pixels following the iso contours at the given `level`. - - These pixels are localized by the bound of the segment generated by the - iso contour algorithm. - - The result is returned as a numpy array storing a list of coordinates y/x. - - .. code-block:: python - - # Example using a mask - shape = 100, 100 - image = numpy.random.random(shape) - mask = numpy.random.random(shape) < 0.01 - pixels = silx.image.marchingsquares.find_pixels(image, 0.5, mask=mask) - - :param numpy.ndarray image: Image to process - :param float level: Level of the requested iso contours. - :param numpy.ndarray mask: An optional mask (a non-zero value invalidate - the pixels of the image) - :returns: An array of coordinates in y/x - :rtype: numpy.ndarray - """ - assert(image is not None) - if mask is not None: - assert(image.shape == mask.shape) - engine = "merge" - impl = _factory(engine, image, mask) - return impl.find_pixels(level) - - -def find_contours(image, level, mask=None): - """ - Find the iso contours at the given `level`. - - The result is returned as a list of polygons. - - .. code-block:: python - - # Example using a mask - shape = 100, 100 - image = numpy.random.random(shape) - mask = numpy.random.random(shape) < 0.01 - polygons = silx.image.marchingsquares.find_contours(image, 0.5, mask=mask) - - :param numpy.ndarray image: Image to process - :param float level: Level of the requested iso contours. - :param numpy.ndarray mask: An optional mask (a non-zero value invalidate - the pixels of the image) - :returns: A list of array containing y-x coordinates of points - :rtype: List[numpy.ndarray] - """ - assert(image is not None) - if mask is not None: - assert(image.shape == mask.shape) - engine = "merge" - impl = _factory(engine, image, mask) - return impl.find_contours(level) diff --git a/silx/image/marchingsquares/_mergeimpl.pyx b/silx/image/marchingsquares/_mergeimpl.pyx deleted file mode 100644 index 5a7a3b5..0000000 --- a/silx/image/marchingsquares/_mergeimpl.pyx +++ /dev/null @@ -1,1319 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2018-2020 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ -""" -Marching squares implementation based on a merge of segements and polygons. -""" - -__authors__ = ["Almar Klein", "Jerome Kieffer", "Valentin Valls"] -__license__ = "MIT" -__date__ = "23/04/2018" - -import numpy -cimport numpy as cnumpy - -from libcpp.vector cimport vector -from libcpp.list cimport list as clist -from libcpp.set cimport set as cset -from libcpp.map cimport map -from libcpp cimport bool -from libc.math cimport fabs -from libc.math cimport floor - -from cython.parallel import prange -from cython.operator cimport dereference -from cython.operator cimport preincrement -cimport libc.stdlib -cimport libc.string - -cimport cython - -from ...utils._have_openmp cimport COMPILED_WITH_OPENMP -"""Store in the module if it was compiled with OpenMP""" - -cdef double EPSILON = numpy.finfo(numpy.float64).eps - -# Windows compatibility: Cross-platform INFINITY -from libc.float cimport DBL_MAX -cdef double INFINITY = DBL_MAX + DBL_MAX -# from libc.math cimport INFINITY - -cdef extern from "include/patterns.h": - cdef unsigned char EDGE_TO_POINT[][2] - cdef unsigned char CELL_TO_EDGE[][5] - cdef struct coord_t: - short x - short y - -ctypedef cnumpy.uint32_t point_index_t -"""Type of the unique index identifying a connection for the polygons.""" - -"""Define a point of a polygon.""" -cdef struct point_t: - cnumpy.float32_t x - cnumpy.float32_t y - -"""Description of a non-final polygon.""" -cdef cppclass PolygonDescription: - point_index_t begin - point_index_t end - clist[point_t] points - - PolygonDescription() nogil: - pass - -"""Description of a tile context. - -It contains structure to store intermediate and final data of a thread. -Pixels and contours structures are merged together as it looks to have -mostly no cost. -""" -cdef cppclass TileContext: - int pos_x - int pos_y - int dim_x - int dim_y - - # Only used to find contours - clist[PolygonDescription*] final_polygons - map[point_index_t, PolygonDescription*] polygons - - # Only used to find pixels - clist[coord_t] final_pixels - cset[coord_t] pixels - - TileContext() nogil: - pass - - -cdef class _MarchingSquaresAlgorithm(object): - """Abstract class managing a marching squares algorithm. - - It provides common methods to execute the process, with the support of - OpenMP, plus some hooks. Mostly created to be able to reuse part of the - logic between `_MarchingSquaresContours` and `_MarchingSquaresPixels`. - """ - - cdef cnumpy.float32_t *_image_ptr - cdef cnumpy.int8_t *_mask_ptr - cdef int _dim_x - cdef int _dim_y - cdef int _group_size - cdef bool _use_minmax_cache - cdef bool _force_sequencial_reduction - - cdef TileContext* _final_context - - cdef cnumpy.float32_t *_min_cache - cdef cnumpy.float32_t *_max_cache - - def __cinit__(self): - self._use_minmax_cache = False - self._force_sequencial_reduction = False - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void marching_squares(self, cnumpy.float64_t level) nogil: - """ - Main method to execute the marching squares. - - :param level: The level expected. - """ - cdef: - TileContext** contexts - TileContext** valid_contexts - int nb_contexts, nb_valid_contexts - int i, j - TileContext* context - int dim_x, dim_y - - contexts = self.create_contexts(level, &dim_x, &dim_y, &nb_valid_contexts) - nb_contexts = dim_x * dim_y - - if nb_valid_contexts == 0: - # shortcut - self._final_context = new TileContext() - libc.stdlib.free(contexts) - return - - j = 0 - valid_contexts = <TileContext **>libc.stdlib.malloc(nb_valid_contexts * sizeof(TileContext*)) - for i in xrange(nb_contexts): - if contexts[i] != NULL: - valid_contexts[j] = contexts[i] - j += 1 - - # openmp - for i in prange(nb_valid_contexts, nogil=True): - self.marching_squares_mp(valid_contexts[i], level) - - if nb_valid_contexts == 1: - # shortcut - self._final_context = valid_contexts[0] - libc.stdlib.free(valid_contexts) - libc.stdlib.free(contexts) - return - - if self._force_sequencial_reduction: - self.sequencial_reduction(nb_valid_contexts, valid_contexts) - # FIXME can only be used if compiled with openmp - # elif copenmp.omp_get_num_threads() <= 1: - # self._sequencial_reduction(nb_valid_contexts, valid_contexts) - else: - self.reduction_2d(dim_x, dim_y, contexts) - - libc.stdlib.free(valid_contexts) - libc.stdlib.free(contexts) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void reduction_2d(self, int dim_x, int dim_y, TileContext **contexts) nogil: - """ - Reduce the problem merging first neighbours together in a recursive - process. Optimized with OpenMP. - - :param dim_x: Number of contexts in the x dimension - :param dim_y: Number of contexts in the y dimension - :param contexts: Array of contexts - """ - cdef: - int x1, y1, x2, y2, i1, i2 - int delta = 1 - - while True: - if delta >= dim_x and delta >= dim_y: - break - # NOTE: Cython 0.21.1 is buggy with prange + steps - # It is needed to add a delta and the 'to' - # Here is what we can use with Cython 0.28: - # for i in prange(0, dim_x, (delta + delta)): - for i1 in prange(0, dim_x + (delta + delta - 1), delta + delta, nogil=True): - x1 = i1 - if x1 + delta < dim_x: - y1 = 0 - while y1 < dim_y: - self.merge_array_contexts(contexts, y1 * dim_x + x1, y1 * dim_x + x1 + delta) - y1 = y1 + delta - - # NOTE: Cython 0.21.1 is buggy with prange + steps - # It is needed to add a delta and the 'to' - # Here is what we can use with Cython 0.28: - # for i in prange(0, dim_y, (delta + delta)): - for i2 in prange(0, dim_y + (delta + delta - 1), delta + delta, nogil=True): - y2 = i2 - if y2 + delta < dim_y: - x2 = 0 - while x2 < dim_x: - self.merge_array_contexts(contexts, y2 * dim_x + x2, (y2 + delta) * dim_x + x2) - x2 = x2 + delta + delta - delta <<= 1 - - self._final_context = contexts[0] - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef inline void merge_array_contexts(self, - TileContext **contexts, - int index1, - int index2) nogil: - """ - Merge contexts from `index2` to `index1` and delete the one from index2. - If the one from index1 was NULL, the one from index2 is moved to index1 - and is not deleted. - - This intermediate function was needed to avoid compilation problem of - Cython + OpenMP. - """ - cdef: - TileContext *context1 - TileContext *context2 - - context1 = contexts[index1] - context2 = contexts[index2] - if context1 != NULL and context2 != NULL: - self.merge_context(context1, context2) - del context2 - elif context2 != NULL: - contexts[index1] = context2 - # for sanity - # contexts[index2] = NULL - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void sequencial_reduction(self, - int nb_contexts, - TileContext **contexts) nogil: - """ - Reduce the problem sequencially without taking care of the topology - - :param nb_contexts: Number of contexts - :param contexts: Array of contexts - """ - cdef: - int i - # merge - self._final_context = new TileContext() - for i in xrange(nb_contexts): - if contexts[i] != NULL: - self.merge_context(self._final_context, contexts[i]) - del contexts[i] - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void marching_squares_mp(self, - TileContext *context, - cnumpy.float64_t level) nogil: - """ - Main entry of the marching squares algorithm for each threads. - - :param context: Context used by the thread to store data - :param level: The requested level - """ - cdef: - int x, y, pattern - cnumpy.float64_t tmpf - cnumpy.float32_t *image_ptr - cnumpy.int8_t *mask_ptr - - image_ptr = self._image_ptr + (context.pos_y * self._dim_x + context.pos_x) - if self._mask_ptr != NULL: - mask_ptr = self._mask_ptr + (context.pos_y * self._dim_x + context.pos_x) - else: - mask_ptr = NULL - - for y in range(context.pos_y, context.pos_y + context.dim_y): - for x in range(context.pos_x, context.pos_x + context.dim_x): - # Calculate index. - pattern = 0 - if image_ptr[0] > level: - pattern += 1 - if image_ptr[1] > level: - pattern += 2 - if image_ptr[self._dim_x] > level: - pattern += 8 - if image_ptr[self._dim_x + 1] > level: - pattern += 4 - - # Resolve ambiguity - if pattern == 5 or pattern == 10: - # Calculate value of cell center (i.e. average of corners) - tmpf = 0.25 * (image_ptr[0] + - image_ptr[1] + - image_ptr[self._dim_x] + - image_ptr[self._dim_x + 1]) - # If below level, swap - if tmpf <= level: - if pattern == 5: - pattern = 10 - else: - pattern = 5 - - # Cache mask information - if mask_ptr != NULL: - # Note: Store the mask in the index. It could be usefull to - # generate accurate segments in some cases, but yet it - # is not used - if mask_ptr[0] > 0: - pattern += 16 - if mask_ptr[1] > 0: - pattern += 32 - if mask_ptr[self._dim_x] > 0: - pattern += 128 - if mask_ptr[self._dim_x + 1] > 0: - pattern += 64 - mask_ptr += 1 - - if pattern < 16 and pattern != 0 and pattern != 15: - self.insert_pattern(context, x, y, pattern, level) - - image_ptr += 1 - - # There is a missing pixel at the end of each rows - image_ptr += self._dim_x - context.dim_x - if mask_ptr != NULL: - mask_ptr += self._dim_x - context.dim_x - - self.after_marching_squares(context) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void after_marching_squares(self, TileContext *context) nogil: - """ - Called by each threads after execution of the marching squares - algorithm. Called before merging together the contextes. - - :param context: Context used by the thread to store data - """ - pass - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void insert_pattern(self, - TileContext *context, - int x, - int y, - int pattern, - cnumpy.float64_t level) nogil: - """ - Called by the marching squares algorithm each time a pattern is found. - - :param context: Context used by the thread to store data - :param x: X location of the pattern - :param y: Y location of the pattern - :param pattern: Binary-field identifying lower and higher pixel levels - :param level: The requested level - """ - pass - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void merge_context(self, - TileContext *context, - TileContext *other) nogil: - """ - Merge into a context another context. - - :param context: Context which will contains the merge result - :param other: Context to merge into the other one. The merging process - is destructive. The context may returns empty. - """ - pass - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef TileContext** create_contexts(self, - cnumpy.float64_t level, - int* dim_x, - int* dim_y, - int* nb_valid_contexts) nogil: - """ - Create and initialize a 2d-array of contexts. - - If the minmax cache is used, only useful context will be created. - Thous with the minmax range excluding the level will not be created and - will have a `NULL` reference in the context array. - - :param level: The requested level - :param dim_x: Resulting X dimension of context array - :param dim_x: Resulting Y dimension of context array - :param nb_valid_contexts: Resulting number of created contexts - :return: The context array - """ - cdef: - int context_dim_x, context_dim_y - int context_size, valid_contexts - int x, y - int icontext - TileContext* context - TileContext** contexts - - context_dim_x = self._dim_x // self._group_size + (self._dim_x % self._group_size > 0) - context_dim_y = self._dim_y // self._group_size + (self._dim_y % self._group_size > 0) - context_size = context_dim_x * context_dim_y - contexts = <TileContext **>libc.stdlib.malloc(context_size * sizeof(TileContext*)) - libc.string.memset(contexts, 0, context_size * sizeof(TileContext*)) - - valid_contexts = 0 - icontext = 0 - y = 0 - while y < self._dim_y - 1: - x = 0 - while x < self._dim_x - 1: - if self._use_minmax_cache: - if level < self._min_cache[icontext] or level > self._max_cache[icontext]: - icontext += 1 - x += self._group_size - continue - context = self.create_context(x, y, self._group_size, self._group_size) - contexts[icontext] = context - icontext += 1 - valid_contexts += 1 - x += self._group_size - y += self._group_size - - # dereference is not working here... then we uses array index but - # it is not the proper way - dim_x[0] = context_dim_x - dim_y[0] = context_dim_y - nb_valid_contexts[0] = valid_contexts - return contexts - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef TileContext *create_context(self, - int x, - int y, - int dim_x, - int dim_y) nogil: - """ - Allocate and initialize a context. - - :param x: Left location of the context into the image - :param y: Top location of the context into the image - :param dim_x: Size of the context in the X dimension of the image - :param dim_y: Size of the context in the Y dimension of the image - :return: The context - """ - cdef: - TileContext *context - context = new TileContext() - context.pos_x = x - context.pos_y = y - context.dim_x = dim_x - context.dim_y = dim_y - if x + context.dim_x > self._dim_x - 1: - context.dim_x = self._dim_x - 1 - x - if y + context.dim_y > self._dim_y - 1: - context.dim_y = self._dim_y - 1 - y - if context.dim_x <= 0 or context.dim_y <= 0: - del context - return NULL - return context - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void compute_point(self, - cnumpy.uint32_t x, - cnumpy.uint32_t y, - cnumpy.uint8_t edge, - cnumpy.float64_t level, - point_t *result_point) nogil: - """ - Compute the location of a point of the polygons according to the level - and the neighbours. - - :param x: X location of the 4-pixels - :param y: Y location of the 4-pixels - :param edge: Enumeration identifying the 2-pixels to process - :param level: The requested level - :param result_point: Resulting value of the point - """ - cdef: - int dx1, dy1, index1 - int dx2, dy2, index2 - cnumpy.float64_t fx, fy, ff, weight1, weight2 - # Use these to look up the relative positions of the pixels to interpolate - dx1, dy1 = EDGE_TO_POINT[edge][0], EDGE_TO_POINT[edge][1] - dx2, dy2 = EDGE_TO_POINT[edge + 1][0], EDGE_TO_POINT[edge + 1][1] - # Define "strength" of each corner of the cube that we need - index1 = (y + dy1) * self._dim_x + x + dx1 - index2 = (y + dy2) * self._dim_x + x + dx2 - weight1 = 1.0 / (EPSILON + fabs(self._image_ptr[index1] - level)) - weight2 = 1.0 / (EPSILON + fabs(self._image_ptr[index2] - level)) - # Apply a kind of center-of-mass method - fx, fy, ff = 0.0, 0.0, 0.0 - fx += dx1 * weight1 - fy += dy1 * weight1 - ff += weight1 - fx += dx2 * weight2 - fy += dy2 * weight2 - ff += weight2 - fx /= ff - fy /= ff - result_point.x = x + fx - result_point.y = y + fy - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void compute_ipoint(self, - cnumpy.uint32_t x, - cnumpy.uint32_t y, - cnumpy.uint8_t edge, - cnumpy.float64_t level, - coord_t *result_coord) nogil: - """ - Compute the location of pixel which contains the point of the polygons - according to the level and the neighbours. - - This implementation is supposed to be faster than `compute_point` when - we only request the location of the pixel. - - :param x: X location of the 4-pixels - :param y: Y location of the 4-pixels - :param edge: Enumeration identifying the 2-pixels to process - :param level: The requested level - :param result_coord: Resulting location of the pixel - """ - cdef: - int dx1, dy1, index1 - int dx2, dy2, index2 - cnumpy.float64_t fx, fy, ff, weight1, weight2 - # Use these to look up the relative positions of the pixels to interpolate - dx1, dy1 = EDGE_TO_POINT[edge][0], EDGE_TO_POINT[edge][1] - dx2, dy2 = EDGE_TO_POINT[edge + 1][0], EDGE_TO_POINT[edge + 1][1] - # Define "strength" of each corner of the cube that we need - index1 = (y + dy1) * self._dim_x + x + dx1 - index2 = (y + dy2) * self._dim_x + x + dx2 - weight1 = EPSILON + fabs(self._image_ptr[index1] - level) - weight2 = EPSILON + fabs(self._image_ptr[index2] - level) - # Apply a kind of center-of-mass method - if edge == 0: - result_coord.x = x + (weight1 > weight2) - result_coord.y = y - elif edge == 1: - result_coord.x = x + 1 - result_coord.y = y + (weight1 > weight2) - elif edge == 2: - result_coord.x = x + (weight1 < weight2) - result_coord.y = y + 1 - elif edge == 3: - result_coord.x = x - result_coord.y = y + (weight1 < weight2) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef point_index_t create_point_index(self, int yx, cnumpy.uint8_t edge) nogil: - """ - Create a unique identifier for a point of a polygon based on the - pattern location and the edge. - - A index can be shared by different pixel coordinates. For example, - the index of the tuple (x=0, y=0, edge=2) is equal to the one of - (x=1, y=0, edge=0). - - :param yx: Index of the location of the pattern in the image - :param edge: Enumeration identifying the edge of the pixel - :return: An index - """ - if edge == 2: - yx += self._dim_x - edge = 0 - elif edge == 1: - yx += 1 - elif edge == 3: - edge = 1 - - # Reserve the zero value - yx += 1 - - return edge + (yx << 1) - - -cdef class _MarchingSquaresContours(_MarchingSquaresAlgorithm): - """Implementation of the marching squares algorithm to find iso contours. - """ - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void insert_pattern(self, - TileContext *context, - int x, - int y, - int pattern, - cnumpy.float64_t level) nogil: - cdef: - int segment - for segment in range(CELL_TO_EDGE[pattern][0]): - begin_edge = CELL_TO_EDGE[pattern][1 + segment * 2 + 0] - end_edge = CELL_TO_EDGE[pattern][1 + segment * 2 + 1] - self.insert_segment(context, x, y, begin_edge, end_edge, level) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void insert_segment(self, TileContext *context, - int x, int y, - cnumpy.uint8_t begin_edge, - cnumpy.uint8_t end_edge, - cnumpy.float64_t level) nogil: - cdef: - int i, yx - point_t point - point_index_t begin, end - PolygonDescription *description - PolygonDescription *description_begin - PolygonDescription *description_end - map[point_index_t, PolygonDescription*].iterator it_begin - map[point_index_t, PolygonDescription*].iterator it_end - - yx = self._dim_x * y + x - begin = self.create_point_index(yx, begin_edge) - end = self.create_point_index(yx, end_edge) - - it_begin = context.polygons.find(begin) - it_end = context.polygons.find(end) - if it_begin == context.polygons.end() and it_end == context.polygons.end(): - # insert a new polygon - description = new PolygonDescription() - description.begin = begin - description.end = end - self.compute_point(x, y, begin_edge, level, &point) - description.points.push_back(point) - self.compute_point(x, y, end_edge, level, &point) - description.points.push_back(point) - context.polygons[begin] = description - context.polygons[end] = description - elif it_begin == context.polygons.end(): - # insert the beginning point to an existing polygon - self.compute_point(x, y, begin_edge, level, &point) - description = dereference(it_end).second - context.polygons.erase(it_end) - if end == description.begin: - # insert at start - description.points.push_front(point) - description.begin = begin - context.polygons[begin] = description - else: - # insert on tail - description.points.push_back(point) - description.end = begin - context.polygons[begin] = description - elif it_end == context.polygons.end(): - # insert the ending point to an existing polygon - self.compute_point(x, y, end_edge, level, &point) - description = dereference(it_begin).second - context.polygons.erase(it_begin) - if begin == description.begin: - # insert at start - description.points.push_front(point) - description.begin = end - context.polygons[end] = description - else: - # insert on tail - description.points.push_back(point) - description.end = end - context.polygons[end] = description - else: - # merge 2 polygons using this segment - description_begin = dereference(it_begin).second - description_end = dereference(it_end).second - if description_begin == description_end: - # The segment closes a polygon - # FIXME: this intermediate assign is not needed - point = description_begin.points.front() - description_begin.points.push_back(point) - context.polygons.erase(begin) - context.polygons.erase(end) - context.final_polygons.push_back(description_begin) - else: - if ((begin == description_begin.begin or end == description_begin.begin) and - (begin == description_end.end or end == description_end.end)): - # worst case, let's make it faster - description = description_end - description_end = description_begin - description_begin = description - - # FIXME: We can recycle a description instead of creating a new one - description = new PolygonDescription() - - # Make sure the last element of the list is the one to connect - if description_begin.begin == begin or description_begin.begin == end: - # O(n) - description_begin.points.reverse() - description.begin = description_begin.end - else: - description.begin = description_begin.begin - - # O(1) - description.points.splice(description.points.end(), description_begin.points) - - # Make sure the first element of the list is the one to connect - if description_end.end == begin or description_end.end == end: - description_end.points.reverse() - description.end = description_end.begin - else: - description.end = description_end.end - - description.points.splice(description.points.end(), description_end.points) - - context.polygons.erase(it_begin) - context.polygons.erase(it_end) - context.polygons[description.begin] = description - context.polygons[description.end] = description - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void merge_context(self, TileContext *context, TileContext *other) nogil: - cdef: - map[point_index_t, PolygonDescription*].iterator it_begin - map[point_index_t, PolygonDescription*].iterator it_end - map[point_index_t, PolygonDescription*].iterator it - PolygonDescription *description_other - PolygonDescription *description - PolygonDescription *description2 - point_index_t point_index - vector[PolygonDescription*] mergeable_polygons - size_t i - - # merge final polygons - context.final_polygons.splice(context.final_polygons.end(), other.final_polygons) - - # mergeable_polygons.reserve(other.polygons.size() / 2) - it = other.polygons.begin() - while it != other.polygons.end(): - point_index = dereference(it).first - description_other = dereference(it).second - if description_other.begin == point_index: - mergeable_polygons.push_back(description_other) - preincrement(it) - - for i in range(mergeable_polygons.size()): - description_other = mergeable_polygons[i] - it_begin = context.polygons.find(description_other.begin) - it_end = context.polygons.find(description_other.end) - - if it_begin == context.polygons.end() and it_end == context.polygons.end(): - # It's a new polygon - context.polygons[description_other.begin] = description_other - context.polygons[description_other.end] = description_other - elif it_end == context.polygons.end(): - # The head of the polygon have to be merged - description = dereference(it_begin).second - context.polygons.erase(description.begin) - context.polygons.erase(description.end) - if description.begin == description_other.begin: - description.begin = description.end - description.points.reverse() - description.end = description_other.end - # remove the dup element - description_other.points.pop_front() - description.points.splice(description.points.end(), description_other.points) - context.polygons[description.begin] = description - context.polygons[description.end] = description - del description_other - elif it_begin == context.polygons.end(): - # The tail of the polygon have to be merged - description = dereference(it_end).second - context.polygons.erase(description.begin) - context.polygons.erase(description.end) - if description.begin == description_other.end: - description.begin = description.end - description.points.reverse() - description.end = description_other.begin - description_other.points.reverse() - # remove the dup element - description_other.points.pop_front() - description.points.splice(description.points.end(), description_other.points) - context.polygons[description.begin] = description - context.polygons[description.end] = description - del description_other - else: - # Both sides have to be merged - description = dereference(it_begin).second - description2 = dereference(it_end).second - if description == description2: - # It became a closed polygon - context.polygons.erase(description.begin) - context.polygons.erase(description.end) - if description.begin == description_other.begin: - description.begin = description.end - description.points.reverse() - description.end = description_other.end - # remove the dup element - description_other.points.pop_front() - description.points.splice(description.points.end(), description_other.points) - context.final_polygons.push_back(description) - del description_other - else: - context.polygons.erase(description.begin) - context.polygons.erase(description.end) - context.polygons.erase(description2.begin) - context.polygons.erase(description2.end) - if description.begin == description_other.begin: - description.begin = description.end - description.points.reverse() - if description2.end == description_other.end: - description.end = description2.begin - description2.points.reverse() - else: - description.end = description2.end - description_other.points.pop_front() - description2.points.pop_front() - description.points.splice(description.points.end(), description_other.points) - description.points.splice(description.points.end(), description2.points) - context.polygons[description.begin] = description - context.polygons[description.end] = description - del description_other - del description2 - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef extract_polygons(self): - cdef: - size_t i - int i_pixel - cnumpy.uint8_t index - map[point_index_t, PolygonDescription*].iterator it - vector[PolygonDescription*] descriptions - clist[point_t].iterator it_points - PolygonDescription *description - cnumpy.float32_t[:, ::1] polygon - - if self._final_context == NULL: - return [] - - # move all the polygons in a final structure - with nogil: - it = self._final_context.polygons.begin() - while it != self._final_context.polygons.end(): - description = dereference(it).second - if dereference(it).first == description.begin: - # polygones are stored 2 times - # only use one - descriptions.push_back(description) - preincrement(it) - self._final_context.polygons.clear() - - descriptions.insert(descriptions.end(), - self._final_context.final_polygons.begin(), - self._final_context.final_polygons.end()) - self._final_context.final_polygons.clear() - - del self._final_context - self._final_context = NULL - - # create result and clean up allocated memory - polygons = [] - for i in range(descriptions.size()): - description = descriptions[i] - polygon = numpy.empty((description.points.size(), 2), dtype=numpy.float32) - it_points = description.points.begin() - i_pixel = 0 - while it_points != description.points.end(): - polygon[i_pixel, 0] = dereference(it_points).y - polygon[i_pixel, 1] = dereference(it_points).x - i_pixel += 1 - preincrement(it_points) - polygons.append(numpy.asarray(polygon)) - del description - - return polygons - - -cdef class _MarchingSquaresPixels(_MarchingSquaresAlgorithm): - """Implementation of the marching squares algorithm to find pixels of the - image containing points of the polygons of the iso contours. - """ - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void insert_pattern(self, - TileContext *context, - int x, - int y, - int pattern, - cnumpy.float64_t level) nogil: - cdef: - int segment - for segment in range(CELL_TO_EDGE[pattern][0]): - begin_edge = CELL_TO_EDGE[pattern][1 + segment * 2 + 0] - end_edge = CELL_TO_EDGE[pattern][1 + segment * 2 + 1] - self.insert_segment(context, x, y, begin_edge, end_edge, level) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void insert_segment(self, TileContext *context, - int x, int y, - cnumpy.uint8_t begin_edge, - cnumpy.uint8_t end_edge, - cnumpy.float64_t level) nogil: - cdef: - coord_t coord - self.compute_ipoint(x, y, begin_edge, level, &coord) - context.pixels.insert(coord) - self.compute_ipoint(x, y, end_edge, level, &coord) - context.pixels.insert(coord) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void after_marching_squares(self, TileContext *context) nogil: - cdef: - coord_t coord - cset[coord_t].iterator it_coord - cset[coord_t].iterator it_coord_erase - pass - - it_coord = context.pixels.begin() - while it_coord != context.pixels.end(): - coord = dereference(it_coord) - if (coord.x > context.pos_x and coord.x < context.pos_x + context.dim_x - 1 and - coord.y > context.pos_y and coord.y < context.pos_y + context.dim_y - 1): - it_coord_erase = it_coord - preincrement(it_coord) - context.pixels.erase(it_coord_erase) - context.final_pixels.push_back(coord) - else: - preincrement(it_coord) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void merge_context(self, TileContext *context, TileContext *other) nogil: - cdef: - cset[coord_t].iterator it_coord - - # merge final pixels - context.final_pixels.splice(context.final_pixels.end(), other.final_pixels) - - # merge final pixels - # NOTE: This is not declared in Cython - # context.final_pixels.insert(other.final_pixels.begin(), other.final_pixels.end()) - it_coord = other.pixels.begin() - while it_coord != other.pixels.end(): - context.pixels.insert(dereference(it_coord)) - preincrement(it_coord) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef extract_pixels(self): - cdef: - int i, x, y - point_index_t index - cset[coord_t].iterator it - clist[coord_t].iterator it_coord - coord_t coord - cnumpy.int32_t[:, ::1] pixels - - if self._final_context == NULL: - return numpy.empty((0, 2), dtype=numpy.int32) - - # create result - it = self._final_context.pixels.begin() - while it != self._final_context.pixels.end(): - coord = dereference(it) - self._final_context.final_pixels.push_back(coord) - preincrement(it) - - pixels = numpy.empty((self._final_context.final_pixels.size(), 2), dtype=numpy.int32) - i = 0 - - it_coord = self._final_context.final_pixels.begin() - while it_coord != self._final_context.final_pixels.end(): - coord = dereference(it_coord) - pixels[i, 0] = coord.y - pixels[i, 1] = coord.x - i += 1 - preincrement(it_coord) - - del self._final_context - self._final_context = NULL - - return numpy.asarray(pixels) - - -cdef class MarchingSquaresMergeImpl(object): - """ - Marching squares implementation based on a merge of segements and polygons. - - The main logic is based on the common marching squares algorithms. - Segments of the iso-valued contours are identified using a pattern based - on blocks of 2*2 pixels. The image is read sequencially and when a segment - is identified it is inserted at a right place is a set of valid polygons. - This process can grow up polygons on bounds, or merge polygons together. - - The algorithm can take care of a mask. If a pixel is invalidated by a - non-zero value of the mask at it's location, the computation of the pattern - cancelled and no segements are generated. - - This implementation based on merge allow to use divide and conquer - implementation in multi process using OpenMP.The image is subdivised into - many tiles, each one is processed independantly. The result is finally - reduced by consecutives polygon merges. - - The OpenMP group size can also by used to skip part of the image using - pre-computed informations. `use_minmax_cache` can enable the computation of - minimum and maximum pixel levels available on each tile groups. It was - designed to improve the efficiency of the extraction of many contour levels - from the same gradient image. - - Finally the implementation provides an implementation to reach polygons - (:meth:`find_contours`) or pixels (:meth:`find_pixels`) from the iso-valued - data. - - .. code-block:: python - - # Example using a mask - shape = 100, 100 - image = numpy.random.random(shape) - mask = numpy.random.random(shape) < 0.01 - ms = MarchingSquaresMergeImpl(image, mask) - polygons = ms.find_contours(level=0.5) - for polygon in polygons: - print(polygon) - - .. code-block:: python - - # Example using multi requests - shape = 1000, 1000 - image = numpy.random.random(shape) - ms = MarchingSquaresMergeImpl(image) - levels = numpy.arange(0, 1, 0.05) - for level in levels: - polygons = ms.find_contours(level=level) - - .. code-block:: python - - # Efficient cache using multi requests - shape = 1000, 1000 - image = numpy.arange(shape[0] * shape[1]) / (shape[0] * shape[1]) - image.shape = shape - ms = MarchingSquaresMergeImpl(image, use_minmax_cache=True) - levels = numpy.arange(0, 1, 0.05) - for level in levels: - polygons = ms.find_contours(level=level) - - :param numpy.ndarray image: Image to process. - If the image is not a continuous array of native float 32bits, the data - will be first normalized. This can reduce efficiency. - :param numpy.ndarray mask: An optional mask (a non-zero value invalidate - the pixels of the image) - If the image is not a continuous array of signed integer 8bits, the - data will be first normalized. This can reduce efficiency. - :param int group_size: Specify the size of the tile to split the - computation with OpenMP. It is also used as tile size to compute the - min/max cache - :param bool use_minmax_cache: If true the min/max cache is enabled. - """ - - cdef cnumpy.float32_t[:, ::1] _image - cdef cnumpy.int8_t[:, ::1] _mask - - cdef cnumpy.float32_t *_image_ptr - cdef cnumpy.int8_t *_mask_ptr - cdef int _dim_x - cdef int _dim_y - cdef int _group_size - cdef bool _use_minmax_cache - - cdef cnumpy.float32_t *_min_cache - cdef cnumpy.float32_t *_max_cache - - cdef _MarchingSquaresContours _contours_algo - cdef _MarchingSquaresPixels _pixels_algo - - def __init__(self, - image, mask=None, - group_size=256, - use_minmax_cache=False): - if not isinstance(image, numpy.ndarray) or len(image.shape) != 2: - raise ValueError("Only 2D arrays are supported.") - if image.shape[0] < 2 or image.shape[1] < 2: - raise ValueError("Input array must be at least 2x2.") - # Force contiguous native array - self._image = numpy.ascontiguousarray(image, dtype='=f4') - self._image_ptr = &self._image[0][0] - if mask is not None: - if not isinstance(mask, numpy.ndarray): - raise ValueError("Only 2D arrays are supported.") - if image.shape != mask.shape: - raise ValueError("Mask size and image size must be the same.") - # Force contiguous native array - self._mask = numpy.ascontiguousarray(mask, dtype='=i1') - self._mask_ptr = &self._mask[0][0] - else: - self._mask = None - self._mask_ptr = NULL - self._group_size = group_size - self._use_minmax_cache = use_minmax_cache - self._min_cache = NULL - self._max_cache = NULL - with nogil: - self._dim_y = self._image.shape[0] - self._dim_x = self._image.shape[1] - self._contours_algo = None - self._pixels_algo = None - - def __dealloc__(self): - if self._min_cache != NULL: - libc.stdlib.free(self._min_cache) - if self._max_cache != NULL: - libc.stdlib.free(self._max_cache) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void _compute_minmax_on_block(self, int block_x, int block_y, int block_index) nogil: - """ - Initialize the minmax cache. - - The cache is computed for each tiles of the image. It reuses the OpenMP - group size for the size of the tile, which allow to skip a full OpenMP - context in case the requested level do not match the cache. - - The minmax is compuded with an overlap of 1 pixel, in order to match - the marching squares algorithm. - - The mask is taking into accound. As result if a tile is fully masked, - the minmax cache result for this tile will have infinit values. - - :param block_x: X location of tile in block unit - :param block_y: Y location of tile in block unit - :param block_index: Index of the tile in the minmax cache structure - """ - cdef: - int x, y - int pos_x, end_x, pos_y, end_y - cnumpy.float32_t minimum, maximum, value - cnumpy.float32_t *image_ptr - cnumpy.int8_t *mask_ptr - - pos_x = block_x * self._group_size - end_x = pos_x + self._group_size + 1 - if end_x > self._dim_x: - end_x = self._dim_x - pos_y = block_y * self._group_size - end_y = pos_y + self._group_size + 1 - if end_y > self._dim_y: - end_y = self._dim_y - - image_ptr = self._image_ptr + (pos_y * self._dim_x + pos_x) - if self._mask_ptr != NULL: - mask_ptr = self._mask_ptr + (pos_y * self._dim_x + pos_x) - else: - mask_ptr = NULL - minimum = INFINITY - maximum = -INFINITY - - for y in range(pos_y, end_y): - for x in range(pos_x, end_x): - if mask_ptr != NULL: - if mask_ptr[0] != 0: - image_ptr += 1 - mask_ptr += 1 - continue - value = image_ptr[0] - if value < minimum: - minimum = value - if value > maximum: - maximum = value - image_ptr += 1 - if mask_ptr != NULL: - mask_ptr += 1 - image_ptr += self._dim_x + pos_x - end_x - if mask_ptr != NULL: - mask_ptr += self._dim_x + pos_x - end_x - - self._min_cache[block_index] = minimum - self._max_cache[block_index] = maximum - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - cdef void _create_minmax_cache(self) nogil: - """ - Create and initialize minmax cache. - """ - cdef: - int icontext, context_x, context_y - int context_dim_x, context_dim_y, context_size - - context_dim_x = self._dim_x // self._group_size + (self._dim_x % self._group_size > 0) - context_dim_y = self._dim_y // self._group_size + (self._dim_y % self._group_size > 0) - context_size = context_dim_x * context_dim_y - - self._min_cache = <cnumpy.float32_t *>libc.stdlib.malloc(context_size * sizeof(cnumpy.float32_t)) - self._max_cache = <cnumpy.float32_t *>libc.stdlib.malloc(context_size * sizeof(cnumpy.float32_t)) - - for icontext in prange(context_size, nogil=True): - context_x = icontext % context_dim_x - context_y = icontext // context_dim_x - self._compute_minmax_on_block(context_x, context_y, icontext) - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - def find_pixels(self, level): - """ - Compute the pixels from the image over the requested iso contours - at this `level`. Pixels are those over the bound of the segments. - - :param float level: Level of the requested iso contours. - :returns: An array of y-x coordinates. - :rtype: numpy.ndarray - """ - if self._use_minmax_cache and self._min_cache == NULL: - self._create_minmax_cache() - - if self._pixels_algo is None: - algo = _MarchingSquaresPixels() - algo._image_ptr = self._image_ptr - algo._mask_ptr = self._mask_ptr - algo._dim_x = self._dim_x - algo._dim_y = self._dim_y - algo._group_size = self._group_size - algo._use_minmax_cache = self._use_minmax_cache - algo._force_sequencial_reduction = COMPILED_WITH_OPENMP == 0 - if self._use_minmax_cache: - algo._min_cache = self._min_cache - algo._max_cache = self._max_cache - self._pixels_algo = algo - else: - algo = self._pixels_algo - - algo.marching_squares(level) - pixels = algo.extract_pixels() - return pixels - - @cython.boundscheck(False) - @cython.wraparound(False) - @cython.cdivision(True) - def find_contours(self, level=None): - """ - Compute the list of polygons of the iso contours at this `level`. - - :param float level: Level of the requested iso contours. - :returns: A list of array containg y-x coordinates of points - :rtype: List[numpy.ndarray] - """ - if self._use_minmax_cache and self._min_cache == NULL: - self._create_minmax_cache() - - if self._contours_algo is None: - algo = _MarchingSquaresContours() - algo._image_ptr = self._image_ptr - algo._mask_ptr = self._mask_ptr - algo._dim_x = self._dim_x - algo._dim_y = self._dim_y - algo._group_size = self._group_size - algo._use_minmax_cache = self._use_minmax_cache - algo._force_sequencial_reduction = COMPILED_WITH_OPENMP == 0 - if self._use_minmax_cache: - algo._min_cache = self._min_cache - algo._max_cache = self._max_cache - self._contours_algo = algo - else: - algo = self._contours_algo - - algo.marching_squares(level) - polygons = algo.extract_polygons() - return polygons diff --git a/silx/image/marchingsquares/_skimage.py b/silx/image/marchingsquares/_skimage.py deleted file mode 100644 index d49eeb0..0000000 --- a/silx/image/marchingsquares/_skimage.py +++ /dev/null @@ -1,139 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2018 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ - -__authors__ = ["V. Valls"] -__license__ = "MIT" -__date__ = "05/04/2018" - - -import numpy -import skimage.measure - - -class MarchingSquaresSciKitImage(object): - """Reference implementation of a marching squares using sci-kit image. - - It uses `skimage.measure.find_contours` to find iso contours taking care of - an optional mask. As result the computation is not accurate but can be used - as reference for benchmark or for testing the API without compiling the - cython part of silx. - - :param numpy.ndarray image: 2d-image containing the values - :param numpy.ndarray mask: Optional 2d-image containing mask to cancel - mask on part of the image. A `0` means the pixel at this location is - valid, else the pixel from the image will not be used. - """ - - def __init__(self, image, mask=None): - self._image = image - self._mask = mask - - _deltas = [(0.0, 0.0), (0.99, 0.0), (0.0, 0.99), (0.99, 0.99)] - - def _flag_coord_over_mask(self, coord): - """Flag coord over the mask as NaN""" - for dx, dy in self._deltas: - if self._mask[int(coord[0] + dx), int(coord[1] + dy)] != 0: - return float("nan"), float("nan") - return coord - - def find_pixels(self, level): - """ - Compute the pixels from the image over the requested iso contours - at this `level`. - - This implementation have to use `skimage.measure.find_contours` then - it is not accurate nor efficient. - - :param float level: Level of the requested iso contours. - :returns: An array of y-x coordinates. - :rtype: numpy.ndarray - """ - polylines = skimage.measure.find_contours(self._image, level=level) - size = 0 - for polyline in polylines: - size += len(polyline) - result = numpy.empty((size, 2), dtype=numpy.int32) - size = 0 - delta = numpy.array([0.5, 0.5]) - for polyline in polylines: - if len(polyline) == 0: - continue - integer_polyline = numpy.floor(polyline + delta) - result[size:size + len(polyline)] = integer_polyline - size += len(polyline) - - if len(result) == 0: - return result - - if self._mask is not None: - # filter out pixels over the mask - x_dim = self._image.shape[1] - indexes = result[:, 0] * x_dim + result[:, 1] - indexes = indexes.ravel() - mask = self._mask.ravel() - indexes = numpy.unique(indexes) - indexes = indexes[mask[indexes] == 0] - pixels = numpy.concatenate((indexes // x_dim, indexes % x_dim)) - pixels.shape = 2, -1 - pixels = pixels.T - result = pixels - else: - # Note: Cound be done using a single line numpy.unique(result, axis=0) - # But here it supports Debian 8 - x_dim = self._image.shape[1] - indexes = result[:, 0] * x_dim + result[:, 1] - indexes = indexes.ravel() - indexes = numpy.unique(indexes) - pixels = numpy.concatenate((indexes // x_dim, indexes % x_dim)) - pixels.shape = 2, -1 - pixels = pixels.T - result = pixels - return result - - def find_contours(self, level): - """ - Compute the list of polygons of the iso contours at this `level`. - - If no mask is involved, the result is the same as - `skimage.measure.find_contours`. - - If the result have to be filtered with a mask, the result is not - accurate nor efficient. Polygons are not splited, but only points are - filtered out using NaN coordinates. This could create artefacts. - - :param float level: Level of the requested iso contours. - :returns: A list of array containg y-x coordinates of points - :rtype: List[numpy.ndarray] - """ - polylines = skimage.measure.find_contours(self._image, level=level) - if self._mask is None: - return polylines - result = [] - for polyline in polylines: - polyline = map(self._flag_coord_over_mask, polyline) - polyline = list(polyline) - polyline = numpy.array(polyline) - result.append(polyline) - return result diff --git a/silx/image/marchingsquares/include/patterns.h b/silx/image/marchingsquares/include/patterns.h deleted file mode 100644 index ff86cc1..0000000 --- a/silx/image/marchingsquares/include/patterns.h +++ /dev/null @@ -1,89 +0,0 @@ -# /*########################################################################## -# -# Copyright (c) 2018 European Synchrotron Radiation Facility -# -# 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. -# -# ###########################################################################*/ - -/** This header provides static lookup table used by marching square - * algorithms. - */ - -#ifndef __MARCHINGSQUARE_PATTERNS_H__ -#define __MARCHINGSQUARE_PATTERNS_H__ - -/** - * Array containing pixel's coordinate linked together by an edge index. - */ -const unsigned char EDGE_TO_POINT[][2] = { - {0, 0}, - {1, 0}, - {1, 1}, - {0, 1}, - {0, 0} -}; - -/** - * Array of index containing 5 values: {nb, seg1b, seg2e, seg2b, seg2e} - * nb: number of segments (up to 2) - * seg1b: index of the start of the 1st edge - * seg1e: index of the end of the 1st edge - * seg2b: index of the start of the 2nd edge - * seg2e: index of the end of the 2nd edge - */ -const unsigned char CELL_TO_EDGE[][5] = { - {0, 0, 0, 0, 0}, // Case 0: 0000: nothing - {1, 0, 3, 0, 0}, // Case 1: 0001 - {1, 0, 1, 0, 0}, // Case 2: 0010 - {1, 1, 3, 0, 0}, // Case 3: 0011 - - {1, 1, 2, 0, 0}, // Case 4: 0100 - {2, 0, 1, 2, 3}, // Case 5: 0101 > ambiguous - {1, 0, 2, 0, 0}, // Case 6: 0110 - {1, 2, 3, 0, 0}, // Case 7: 0111 - - {1, 2, 3, 0, 0}, // Case 8: 1000 - {1, 0, 2, 0, 0}, // Case 9: 1001 - {2, 0, 3, 1, 2}, // Case 10: 1010 > ambiguous - {1, 1, 2, 0, 0}, // Case 11: 1011 - - {1, 1, 3, 0, 0}, // Case 12: 1100 - {1, 0, 1, 0, 0}, // Case 13: 1101 - {1, 0, 3, 0, 0}, // Case 14: 1110 - {0, 0, 0, 0, 0}, // Case 15: 1111 -}; - - -typedef struct coord_t { - short x; - short y; - - bool operator<(const coord_t& other) const { - if (y < other.y) { - return true; - } else if (y == other.y) { - return x < other.x; - } else { - return false; - } - } -} coord_t; - -#endif /*__MARCHINGSQUARE_PATTERNS_H__*/ diff --git a/silx/image/marchingsquares/setup.py b/silx/image/marchingsquares/setup.py deleted file mode 100644 index 902f297..0000000 --- a/silx/image/marchingsquares/setup.py +++ /dev/null @@ -1,51 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2016-2017 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ - -__authors__ = ["V. Valls"] -__license__ = "MIT" -__date__ = "23/04/2018" - -import os -import numpy -from numpy.distutils.misc_util import Configuration - - -def configuration(parent_package='', top_path=None): - config = Configuration('marchingsquares', parent_package, top_path) - config.add_subpackage('test') - - silx_include = os.path.join(top_path, "silx", "utils", "include") - config.add_extension('_mergeimpl', - sources=['_mergeimpl.pyx'], - include_dirs=[numpy.get_include(), silx_include], - language='c++', - extra_link_args=['-fopenmp'], - extra_compile_args=['-fopenmp']) - - return config - - -if __name__ == "__main__": - from numpy.distutils.core import setup - setup(configuration=configuration) diff --git a/silx/image/marchingsquares/test/__init__.py b/silx/image/marchingsquares/test/__init__.py deleted file mode 100644 index 5351a28..0000000 --- a/silx/image/marchingsquares/test/__init__.py +++ /dev/null @@ -1,40 +0,0 @@ -# -*- coding: utf-8 -*- -# -# Project: silx -# https://github.com/silx-kit/silx -# -# Copyright (C) 2012-2016 European Synchrotron Radiation Facility, Grenoble, France -# -# 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. - -__authors__ = ["V. Valls"] -__license__ = "MIT" -__date__ = "17/04/2018" - -import unittest -from . import test_funcapi -from . import test_mergeimpl - - -def suite(): - """Test suite for module silx.image.test""" - test_suite = unittest.TestSuite() - test_suite.addTest(test_funcapi.suite()) - test_suite.addTest(test_mergeimpl.suite()) - return test_suite diff --git a/silx/image/marchingsquares/test/test_funcapi.py b/silx/image/marchingsquares/test/test_funcapi.py deleted file mode 100644 index a84a493..0000000 --- a/silx/image/marchingsquares/test/test_funcapi.py +++ /dev/null @@ -1,99 +0,0 @@ -# -*- coding: utf-8 -*- -# -# Project: silx -# https://github.com/silx-kit/silx -# -# Copyright (C) 2012-2016 European Synchrotron Radiation Facility, Grenoble, France -# -# 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. - -__authors__ = ["V. Valls"] -__license__ = "MIT" -__date__ = "17/04/2018" - -import unittest -import numpy -import silx.image.marchingsquares - - -class MockMarchingSquares(object): - - last = None - - def __init__(self, image, mask=None): - MockMarchingSquares.last = self - self.events = [] - self.events.append(("image", image)) - self.events.append(("mask", mask)) - - def find_pixels(self, level): - self.events.append(("find_pixels", level)) - return None - - def find_contours(self, level): - self.events.append(("find_contours", level)) - return None - - -class TestFunctionalApi(unittest.TestCase): - """Test that the default functional API is called using the right - parameters to the right location.""" - - def setUp(self): - self.old_impl = silx.image.marchingsquares.MarchingSquaresMergeImpl - silx.image.marchingsquares.MarchingSquaresMergeImpl = MockMarchingSquares - - def tearDown(self): - silx.image.marchingsquares.MarchingSquaresMergeImpl = self.old_impl - del self.old_impl - - def test_default_find_contours(self): - image = numpy.ones((2, 2), dtype=numpy.float32) - mask = numpy.zeros((2, 2), dtype=numpy.int32) - level = 2.5 - silx.image.marchingsquares.find_contours(image=image, level=level, mask=mask) - events = MockMarchingSquares.last.events - self.assertEqual(len(events), 3) - self.assertEqual(events[0][0], "image") - self.assertEqual(events[0][1][0, 0], 1) - self.assertEqual(events[1][0], "mask") - self.assertEqual(events[1][1][0, 0], 0) - self.assertEqual(events[2][0], "find_contours") - self.assertEqual(events[2][1], level) - - def test_default_find_pixels(self): - image = numpy.ones((2, 2), dtype=numpy.float32) - mask = numpy.zeros((2, 2), dtype=numpy.int32) - level = 3.5 - silx.image.marchingsquares.find_pixels(image=image, level=level, mask=mask) - events = MockMarchingSquares.last.events - self.assertEqual(len(events), 3) - self.assertEqual(events[0][0], "image") - self.assertEqual(events[0][1][0, 0], 1) - self.assertEqual(events[1][0], "mask") - self.assertEqual(events[1][1][0, 0], 0) - self.assertEqual(events[2][0], "find_pixels") - self.assertEqual(events[2][1], level) - - -def suite(): - test_suite = unittest.TestSuite() - loadTests = unittest.defaultTestLoader.loadTestsFromTestCase - test_suite.addTest(loadTests(TestFunctionalApi)) - return test_suite diff --git a/silx/image/marchingsquares/test/test_mergeimpl.py b/silx/image/marchingsquares/test/test_mergeimpl.py deleted file mode 100644 index 1c14f33..0000000 --- a/silx/image/marchingsquares/test/test_mergeimpl.py +++ /dev/null @@ -1,272 +0,0 @@ -# -*- coding: utf-8 -*- -# -# Project: silx -# https://github.com/silx-kit/silx -# -# Copyright (C) 2012-2016 European Synchrotron Radiation Facility, Grenoble, France -# -# 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. - -__authors__ = ["V. Valls"] -__license__ = "MIT" -__date__ = "18/04/2018" - -import unittest -import numpy -from .._mergeimpl import MarchingSquaresMergeImpl - - -class TestMergeImplApi(unittest.TestCase): - - def test_image_not_an_array(self): - bad_image = 1 - self.assertRaises(ValueError, MarchingSquaresMergeImpl, bad_image) - - def test_image_bad_dim(self): - bad_image = numpy.array([[[1.0]]]) - self.assertRaises(ValueError, MarchingSquaresMergeImpl, bad_image) - - def test_image_not_big_enough(self): - bad_image = numpy.array([[1.0, 1.0, 1.0, 1.0]]) - self.assertRaises(ValueError, MarchingSquaresMergeImpl, bad_image) - - def test_mask_not_an_array(self): - image = numpy.array([[1.0, 1.0], [1.0, 1.0]]) - bad_mask = 1 - self.assertRaises(ValueError, MarchingSquaresMergeImpl, image, bad_mask) - - def test_mask_not_match(self): - image = numpy.array([[1.0, 1.0], [1.0, 1.0]]) - bad_mask = numpy.array([[1.0, 1.0]]) - self.assertRaises(ValueError, MarchingSquaresMergeImpl, image, bad_mask) - - def test_ok_anyway_bad_type(self): - image = numpy.array([[1.0, 1.0], [1.0, 1.0]], dtype=numpy.int32) - mask = numpy.array([[1.0, 1.0], [1.0, 1.0]], dtype=numpy.float32) - MarchingSquaresMergeImpl(image, mask) - - def test_find_contours_result(self): - image = numpy.zeros((2, 2)) - image[0, 0] = 1 - ms = MarchingSquaresMergeImpl(image) - polygons = ms.find_contours(0.5) - self.assertIsInstance(polygons, list) - self.assertTrue(len(polygons), 1) - self.assertIsInstance(polygons[0], numpy.ndarray) - self.assertEqual(polygons[0].shape[1], 2) - self.assertEqual(polygons[0].dtype.kind, "f") - - def test_find_pixels_result(self): - image = numpy.zeros((2, 2)) - image[0, 0] = 1 - ms = MarchingSquaresMergeImpl(image) - pixels = ms.find_pixels(0.5) - self.assertIsInstance(pixels, numpy.ndarray) - self.assertEqual(pixels.shape[1], 2) - self.assertEqual(pixels.dtype.kind, "i") - - def test_find_contours_empty_result(self): - image = numpy.zeros((2, 2)) - ms = MarchingSquaresMergeImpl(image) - polygons = ms.find_contours(0.5) - self.assertIsInstance(polygons, list) - self.assertEqual(len(polygons), 0) - - def test_find_pixels_empty_result(self): - image = numpy.zeros((2, 2)) - ms = MarchingSquaresMergeImpl(image) - pixels = ms.find_pixels(0.5) - self.assertIsInstance(pixels, numpy.ndarray) - self.assertEqual(pixels.shape[1], 2) - self.assertEqual(pixels.shape[0], 0) - self.assertEqual(pixels.dtype.kind, "i") - - def test_find_contours_yx_result(self): - image = numpy.zeros((2, 2)) - image[1, 0] = 1 - ms = MarchingSquaresMergeImpl(image) - polygons = ms.find_contours(0.5) - polygon = polygons[0] - self.assertTrue((polygon == (0.5, 0)).any()) - self.assertTrue((polygon == (1, 0.5)).any()) - - def test_find_pixels_yx_result(self): - image = numpy.zeros((2, 2)) - image[1, 0] = 1 - ms = MarchingSquaresMergeImpl(image) - pixels = ms.find_pixels(0.5) - self.assertTrue((pixels == (1, 0)).any()) - - -class TestMergeImplContours(unittest.TestCase): - - def test_merge_segments(self): - image = numpy.zeros((4, 4)) - image[(2, 3), :] = 1 - ms = MarchingSquaresMergeImpl(image) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 1) - - def test_merge_segments_2(self): - image = numpy.zeros((4, 4)) - image[(2, 3), :] = 1 - image[2, 2] = 0 - ms = MarchingSquaresMergeImpl(image) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 1) - - def test_merge_tiles(self): - image = numpy.zeros((4, 4)) - image[(2, 3), :] = 1 - ms = MarchingSquaresMergeImpl(image, group_size=2) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 1) - - def test_fully_masked(self): - image = numpy.zeros((5, 5)) - image[(2, 3), :] = 1 - mask = numpy.ones((5, 5)) - ms = MarchingSquaresMergeImpl(image, mask) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 0) - - def test_fully_masked_minmax(self): - """This invalidates all the tiles. The route is not the same.""" - image = numpy.zeros((5, 5)) - image[(2, 3), :] = 1 - mask = numpy.ones((5, 5)) - ms = MarchingSquaresMergeImpl(image, mask, group_size=2, use_minmax_cache=True) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 0) - - def test_masked_segments(self): - image = numpy.zeros((5, 5)) - image[(2, 3, 4), :] = 1 - mask = numpy.zeros((5, 5)) - mask[:, 2] = 1 - ms = MarchingSquaresMergeImpl(image, mask) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 2) - - def test_closed_polygon(self): - image = numpy.zeros((5, 5)) - image[2, 2] = 1 - image[1, 2] = 1 - image[3, 2] = 1 - image[2, 1] = 1 - image[2, 3] = 1 - mask = None - ms = MarchingSquaresMergeImpl(image, mask) - polygons = ms.find_contours(0.9) - self.assertEqual(len(polygons), 1) - self.assertEqual(list(polygons[0][0]), list(polygons[0][-1])) - - def test_closed_polygon_between_tiles(self): - image = numpy.zeros((5, 5)) - image[2, 2] = 1 - image[1, 2] = 1 - image[3, 2] = 1 - image[2, 1] = 1 - image[2, 3] = 1 - mask = None - ms = MarchingSquaresMergeImpl(image, mask, group_size=2) - polygons = ms.find_contours(0.9) - self.assertEqual(len(polygons), 1) - self.assertEqual(list(polygons[0][0]), list(polygons[0][-1])) - - def test_open_polygon(self): - image = numpy.zeros((5, 5)) - image[2, 2] = 1 - image[1, 2] = 1 - image[3, 2] = 1 - image[2, 1] = 1 - image[2, 3] = 1 - mask = numpy.zeros((5, 5)) - mask[1, 1] = 1 - ms = MarchingSquaresMergeImpl(image, mask) - polygons = ms.find_contours(0.9) - self.assertEqual(len(polygons), 1) - self.assertNotEqual(list(polygons[0][0]), list(polygons[0][-1])) - - def test_ambiguous_pattern(self): - image = numpy.zeros((6, 8)) - image[(3, 4), :] = 1 - image[:, (0, -1)] = 0 - image[3, 3] = -0.001 - image[4, 4] = 0.0 - mask = None - ms = MarchingSquaresMergeImpl(image, mask) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 2) - - def test_ambiguous_pattern_2(self): - image = numpy.zeros((6, 8)) - image[(3, 4), :] = 1 - image[:, (0, -1)] = 0 - image[3, 3] = +0.001 - image[4, 4] = 0.0 - mask = None - ms = MarchingSquaresMergeImpl(image, mask) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 1) - - def count_closed_polygons(self, polygons): - closed = 0 - for polygon in polygons: - if list(polygon[0]) == list(polygon[-1]): - closed += 1 - return closed - - def test_image(self): - # example from skimage - x, y = numpy.ogrid[-numpy.pi:numpy.pi:100j, -numpy.pi:numpy.pi:100j] - image = numpy.sin(numpy.exp((numpy.sin(x)**3 + numpy.cos(y)**2))) - mask = None - ms = MarchingSquaresMergeImpl(image, mask) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 11) - self.assertEqual(self.count_closed_polygons(polygons), 3) - - def test_image_tiled(self): - # example from skimage - x, y = numpy.ogrid[-numpy.pi:numpy.pi:100j, -numpy.pi:numpy.pi:100j] - image = numpy.sin(numpy.exp((numpy.sin(x)**3 + numpy.cos(y)**2))) - mask = None - ms = MarchingSquaresMergeImpl(image, mask, group_size=50) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 11) - self.assertEqual(self.count_closed_polygons(polygons), 3) - - def test_image_tiled_minmax(self): - # example from skimage - x, y = numpy.ogrid[-numpy.pi:numpy.pi:100j, -numpy.pi:numpy.pi:100j] - image = numpy.sin(numpy.exp((numpy.sin(x)**3 + numpy.cos(y)**2))) - mask = None - ms = MarchingSquaresMergeImpl(image, mask, group_size=50, use_minmax_cache=True) - polygons = ms.find_contours(0.5) - self.assertEqual(len(polygons), 11) - self.assertEqual(self.count_closed_polygons(polygons), 3) - - -def suite(): - test_suite = unittest.TestSuite() - loadTests = unittest.defaultTestLoader.loadTestsFromTestCase - test_suite.addTest(loadTests(TestMergeImplApi)) - test_suite.addTest(loadTests(TestMergeImplContours)) - return test_suite diff --git a/silx/image/medianfilter.py b/silx/image/medianfilter.py deleted file mode 100644 index 857f73d..0000000 --- a/silx/image/medianfilter.py +++ /dev/null @@ -1,114 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2017-2018 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ - -""" -This module provides :func:`medfilt2d`, a 2D median filter function -with the choice between 2 implementations: 'cpp' and 'opencl'. -""" - -__authors__ = ["H. Payno"] -__license__ = "MIT" -__date__ = "04/05/2017" - - -import logging - -from silx.math import medianfilter as medianfilter_cpp -from silx.opencl import ocl as _ocl -if _ocl is not None: - from silx.opencl import medfilt as medfilt_opencl -else: # No OpenCL device or pyopencl not installed - medfilt_opencl = None - - -_logger = logging.getLogger(__name__) - - -MEDFILT_ENGINES = ['cpp', 'opencl'] - - -def medfilt2d(image, kernel_size=3, engine='cpp'): - """Apply a median filter on an image. - - This median filter is using a 'nearest' padding for values - past the array edges. If you want more padding options or - functionalities for the median filter (conditional filter - for example) please have a look at - :mod:`silx.math.medianfilter`. - - :param numpy.ndarray image: the 2D array for which we want to apply - the median filter. - :param kernel_size: the dimension of the kernel. - Kernel size must be odd. - If a scalar is given, then it is used as the size in both dimension. - Default: (3, 3) - :type kernel_size: A int or a list of 2 int (kernel_height, kernel_width) - :param engine: the type of implementation to use. - Valid values are: 'cpp' (default) and 'opencl' - - :returns: the array with the median value for each pixel. - - .. note:: if the opencl implementation is requested but - is not present or fails, the cpp implementation is called. - - """ - if engine not in MEDFILT_ENGINES: - err = 'silx doesn\'t have an implementation for the requested engine: ' - err += '%s' % engine - raise ValueError(err) - - if len(image.shape) != 2: - raise ValueError('medfilt2d deals with arrays of dimension 2 only') - - if engine == 'cpp': - return medianfilter_cpp.medfilt(data=image, - kernel_size=kernel_size, - conditional=False) - elif engine == 'opencl': - if medfilt_opencl is None: - wrn = 'opencl median filter not available. ' - wrn += 'Launching cpp implementation.' - _logger.warning(wrn) - # instead call the cpp implementation - return medianfilter_cpp.medfilt(data=image, - kernel_size=kernel_size, - conditional=False) - else: - try: - medianfilter = medfilt_opencl.MedianFilter2D(image.shape, - devicetype="gpu") - res = medianfilter.medfilt2d(image, kernel_size) - except(RuntimeError, MemoryError, ImportError): - wrn = 'Exception occured in opencl median filter. ' - wrn += 'To get more information see debug log.' - wrn += 'Launching cpp implementation.' - _logger.warning(wrn) - _logger.debug("median filter - openCL implementation issue.", - exc_info=True) - # instead call the cpp implementation - res = medianfilter_cpp.medfilt(data=image, - kernel_size=kernel_size, - conditional=False) - - return res diff --git a/silx/image/phantomgenerator.py b/silx/image/phantomgenerator.py deleted file mode 100644 index 10b249b..0000000 --- a/silx/image/phantomgenerator.py +++ /dev/null @@ -1,160 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# -# Copyright (c) 2016 European Synchrotron Radiation Facility -# -# 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. -# -# ###########################################################################*/ - -__authors__ = ["N. Vigano", "H. Payno", "P. Paleo"] -__license__ = "MIT" -__date__ = "19/09/2017" - -import numpy - - -class PhantomGenerator(object): - """ - Class for generating Phantoms - """ - - class _Ellipsoid: - def __init__(self, a, b, c, x0, y0, z0, alpha, mu): - self.a = a - self.b = b - self.c = c - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.alpha = alpha * numpy.pi / 180.0 - self.mu = mu - self.cosAlpha = numpy.cos(self.alpha) - self.sinAlpha = numpy.sin(self.alpha) - - SHEPP_LOGAN = [ - # a b c x0 y0 z0 alpha mu - _Ellipsoid(0.69, 0.92, 0.90, 0.0, 0.0, 0.0, 0.0, 0.10), - _Ellipsoid(0.6624, 0.874, 0.88, 0.0, -0.02, 0.0, 0.0, -0.08), - _Ellipsoid(0.11, 0.31, 0.21, 0.22, -0.0, 0.0, -18.0, -0.02), - _Ellipsoid(0.16, 0.41, 0.22, -0.22, 0.0, -0.25, 18.0, -0.02), - _Ellipsoid(0.21, 0.25, 0.35, 0.0, 0.35, -0.25, 0.0, 0.03), - _Ellipsoid(0.046, 0.046, 0.046, 0.0, 0.10, -0.25, 0.0, 0.01), - _Ellipsoid(0.046, 0.046, 0.02, 0.0, -0.10, -0.25, 0.0, 0.01), - _Ellipsoid(0.046, 0.023, 0.02, -0.08, -0.605, -0.25, 0.0, 0.01), - _Ellipsoid(0.023, 0.023, 0.10, 0.0, -0.605, -0.25, 0.0, 0.01), - _Ellipsoid(0.023, 0.046, 0.10, 0.06, -0.605, -0.25, 0.0, 0.01) - ] - - @staticmethod - def get2DPhantomSheppLogan(n, ellipsoidID=None): - """ - generate a classical 2D shepp logan phantom. - - :param n: The width (and height) of the phantom to generate - :param ellipsoidID: The Id of the ellipsoid to pick. If None will - produce every ellipsoid - :return numpy.ndarray: shepp logan phantom - """ - assert(ellipsoidID is None or (ellipsoidID >= 0 and ellipsoidID < len(PhantomGenerator.SHEPP_LOGAN))) - if ellipsoidID is None: - area = PhantomGenerator._get2DPhantom(n, - PhantomGenerator.SHEPP_LOGAN) - else: - area = PhantomGenerator._get2DPhantom(n, - [PhantomGenerator.SHEPP_LOGAN[ellipsoidID]]) - - indices = numpy.abs(area) > 0 - area[indices] = numpy.multiply(area[indices] + 0.1, 5) - return area / 100.0 - - @staticmethod - def _get2DPhantom(n, phantomSpec): - area = numpy.ndarray(shape=(n, n)) - area.fill(0.) - - count = 0 - for ell in phantomSpec: - count = count+1 - for x in range(n): - sumSquareXandY = PhantomGenerator._getSquareXandYsum(n, x, ell) - indices = sumSquareXandY <= 1 - area[indices, x] = ell.mu - return area - - @staticmethod - def _getSquareXandYsum(n, x, ell): - supportX1 = numpy.ndarray(shape=(n, )) - supportX2 = numpy.ndarray(shape=(n, )) - support_consts = numpy.ndarray(shape=(n, )) - - xScaled = float(2*x-n)/float(n) - xCos = xScaled * ell.cosAlpha - xSin = -xScaled * ell.sinAlpha - supportX1.fill(xCos) - supportX2.fill(xSin) - - supportY1 = numpy.arange(n) - support_consts.fill(2.) - supportY1 = numpy.multiply(support_consts, supportY1) - support_consts.fill(n) - supportY1 = numpy.subtract(supportY1, support_consts) - support_consts.fill(n) - supportY1 = numpy.divide(supportY1, support_consts) - supportY2 = numpy.array(supportY1) - - support_consts.fill(ell.sinAlpha) - supportY1 = numpy.add(supportX1, - numpy.multiply(supportY1, support_consts)) - support_consts.fill(ell.cosAlpha) - supportY2 = numpy.add(supportX2, - numpy.multiply(supportY2, support_consts)) - - support_consts.fill(ell.x0) - supportY1 = numpy.subtract(supportY1, support_consts) - support_consts.fill(ell.y0) - supportY2 = numpy.subtract(supportY2, support_consts) - - support_consts.fill(ell.a) - supportY1 = numpy.power((numpy.divide(supportY1, support_consts)), - 2) - support_consts.fill(ell.b) - supportY2 = numpy.power(numpy.divide(supportY2, support_consts), - 2) - - return numpy.add(supportY1, supportY2) - - @staticmethod - def _getSquareZ(n, ell): - supportZ1 = numpy.arange(n) - support_consts = numpy.ndarray(shape=(n, )) - support_consts.fill(2.) - supportZ1 = numpy.multiply(support_consts, supportZ1) - support_consts.fill(n) - supportZ1 = numpy.subtract(supportZ1, support_consts) - support_consts.fill(n) - supportZ1 = numpy.divide(supportZ1, support_consts) - - support_consts.fill(ell.z0) - supportZ1 = numpy.subtract(supportZ1, ell.z0) - - support_consts.fill(ell.c) - return numpy.power(numpy.divide(supportZ1, support_consts), - 2) - diff --git a/silx/image/projection.py b/silx/image/projection.py deleted file mode 100644 index 5c76c35..0000000 --- a/silx/image/projection.py +++ /dev/null @@ -1,25 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2017 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ - -from silx.opencl.projection import * diff --git a/silx/image/reconstruction.py b/silx/image/reconstruction.py deleted file mode 100644 index 875b66b..0000000 --- a/silx/image/reconstruction.py +++ /dev/null @@ -1,25 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2017 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ - -from silx.opencl.reconstruction import * diff --git a/silx/image/setup.py b/silx/image/setup.py deleted file mode 100644 index 69d5b1b..0000000 --- a/silx/image/setup.py +++ /dev/null @@ -1,47 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2016-2018 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ - -__authors__ = ["J. Kieffer"] -__license__ = "MIT" -__date__ = "05/04/2018" - -from numpy.distutils.misc_util import Configuration - - -def configuration(parent_package='', top_path=None): - config = Configuration('image', parent_package, top_path) - config.add_subpackage('test') - config.add_extension('bilinear', - sources=["bilinear.pyx"], - language='c') - config.add_extension('shapes', - sources=["shapes.pyx"], - language='c') - config.add_subpackage('marchingsquares') - return config - - -if __name__ == "__main__": - from numpy.distutils.core import setup - setup(configuration=configuration) diff --git a/silx/image/shapes.pyx b/silx/image/shapes.pyx deleted file mode 100644 index 9284811..0000000 --- a/silx/image/shapes.pyx +++ /dev/null @@ -1,321 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# -# Copyright (c) 2015-2017 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ -"""This module provides functions making masks on an image. - -- :func:`circle_fill` function generates coordinates of a circle in an image. -- :func:`draw_line` function generates coordinates of a line in an image. -- :func:`polygon_fill_mask` function generates a mask from a set of points - defining a polygon. - -The :class:`Polygon` class provides checking if a point is inside a polygon. - -The whole module uses the (row, col) (i.e., (y, x))) convention -for 2D coordinates. -""" - -__authors__ = ["Jérôme Kieffer", "T. Vincent"] -__license__ = "MIT" -__date__ = "15/02/2019" -__status__ = "dev" - - -cimport cython -import numpy -from libc.math cimport ceil, fabs - - -cdef class Polygon(object): - """Define a polygon that provides inside check and mask generation. - - :param vertices: corners of the polygon - :type vertices: Nx2 array of floats of (row, col) - """ - - cdef float[:,:] vertices - cdef int nvert - - def __init__(self, vertices): - self.vertices = numpy.ascontiguousarray(vertices, dtype=numpy.float32) - if self.vertices.ndim != 2 or self.vertices.shape[1] != 2: - raise ValueError("A list of 2d vertices is expected (n,2 dimensional array)") - self.nvert = self.vertices.shape[0] - - def is_inside(self, row, col): - """Check if (row, col) is inside or outside the polygon - - :param float row: - :param float col: - :return: True if position is inside polygon, False otherwise - """ - return self.c_is_inside(row, col) - - @cython.cdivision(True) - @cython.wraparound(False) - @cython.boundscheck(False) - cdef bint c_is_inside(self, float row, float col) nogil: - """Check if (row, col) is inside or outside the polygon - - Pure C_Cython class implementation. - - :param float row: - :param float col: - :return: True if position is inside polygon, False otherwise - """ - cdef int index, is_inside - cdef float pt1x, pt1y, pt2x, pt2y, xinters - is_inside = 0 - - pt1x = self.vertices[self.nvert-1, 1] - pt1y = self.vertices[self.nvert-1, 0] - for index in range(self.nvert): - pt2x = self.vertices[index, 1] - pt2y = self.vertices[index, 0] - - if (((pt1y <= row and row < pt2y) or - (pt2y <= row and row < pt1y)) and - # Extra (optional) condition to avoid some computation - (col <= pt1x or col <= pt2x)): - xinters = (row - pt1y) * (pt2x - pt1x) / (pt2y - pt1y) + pt1x - is_inside ^= col < xinters - pt1x, pt1y = pt2x, pt2y - return is_inside - - @cython.cdivision(True) - @cython.wraparound(False) - @cython.boundscheck(False) - def make_mask(self, int height, int width): - """Create a mask array representing the filled polygon - - :param int height: Height of the mask array - :param int width: Width of the mask array - :return: 2D array (height, width) - """ - cdef unsigned char[:, :] mask = numpy.zeros((height, width), - dtype=numpy.uint8) - cdef int row_min, row_max, col_min, col_max # mask subpart to update - cdef int row, col, index # Loop indixes - cdef float pt1x, pt1y, pt2x, pt2y # segment end points - cdef int xinters, is_inside, current - - row_min = max(int(min(self.vertices[:, 0])), 0) - row_max = min(int(max(self.vertices[:, 0])) + 1, height) - - # Can be replaced by prange(row_min, row_max, nogil=True) - with nogil: - for row in range(row_min, row_max): - # For each line of the image, mark intersection of all segments - # in the line and then run a xor scan to fill inner parts - # Adapted from http://alienryderflex.com/polygon_fill/ - pt1x = self.vertices[self.nvert-1, 1] - pt1y = self.vertices[self.nvert-1, 0] - col_min = width - 1 - col_max = 0 - is_inside = 0 # Init with whether first col is inside or not - - for index in range(self.nvert): - pt2x = self.vertices[index, 1] - pt2y = self.vertices[index, 0] - - if ((pt1y <= row and row < pt2y) or - (pt2y <= row and row < pt1y)): - # Intersection casted to int so that ]x, x+1] => x - xinters = (<int>ceil(pt1x + (row - pt1y) * - (pt2x - pt1x) / (pt2y - pt1y))) - 1 - - # Update column range to patch - if xinters < col_min: - col_min = xinters - if xinters > col_max: - col_max = xinters - - if xinters < 0: - # Add an intersection to init value of xor scan - is_inside ^= 1 - elif xinters < width: - # Mark intersection in mask - mask[row, xinters] ^= 1 - # else: do not consider intersection on the right - - pt1x, pt1y = pt2x, pt2y - - if col_min < col_max: - # Clip column range to mask - if col_min < 0: - col_min = 0 - if col_max > width - 1: - col_max = width - 1 - - # xor exclusive scan - for col in range(col_min, col_max + 1): - current = mask[row, col] - mask[row, col] = is_inside - is_inside = current ^ is_inside - - # Ensures the result is exported as numpy array and not memory view. - return numpy.asarray(mask) - - -def polygon_fill_mask(vertices, shape): - """Return a mask of boolean, True for pixels inside a polygon. - - :param vertices: Strip of segments end points (row, column) or (y, x) - :type vertices: numpy.ndarray like container of dimension Nx2 - :param shape: size of the mask as (height, width) - :type shape: 2-tuple of int - :return: Mask corresponding to the polygon - :rtype: numpy.ndarray of dimension shape - """ - return Polygon(vertices).make_mask(shape[0], shape[1]) - - -@cython.wraparound(False) -@cython.boundscheck(False) -def draw_line(int row0, int col0, int row1, int col1, int width=1): - """Line includes both end points. - Width is handled by drawing parallel lines, so junctions of lines belonging - to different octant with width > 1 will not look nice. - - Using Bresenham line algorithm: - Bresenham, J. E. - Algorithm for computer control of a digital plotter. - IBM Systems Journal. Vol 4 No 1. 1965. pp 25-30 - - :param int row0: Start point row - :param int col0: Start point col - :param int row1: End point row - :param int col1: End point col - :param int width: Thickness of the line in pixels (default 1) - Width must be at least 1. - :return: Array coordinates of points inside the line (might be negative) - :rtype: 2-tuple of numpy.ndarray (rows, cols) - """ - cdef int drow, dcol, invert_coords - cdef int db, da, delta, b, a, step_a, step_b - cdef int index, offset # Loop indices - - # Store coordinates of points of the line - cdef int[:, :] b_coords - cdef int[:, :] a_coords - - dcol = abs(col1 - col0) - drow = abs(row1 - row0) - invert_coords = dcol < drow - - if dcol == 0 and drow == 0: - return (numpy.array((row0,), dtype=numpy.int32), - numpy.array((col0,), dtype=numpy.int32)) - - if width < 1: - width = 1 - - # Set a and b according to segment octant - if not invert_coords: - da = dcol - db = drow - step_a = 1 if col1 > col0 else -1 - step_b = 1 if row1 > row0 else -1 - a = col0 - b = row0 - - else: - da = drow - db = dcol - step_a = 1 if row1 > row0 else -1 - step_b = 1 if col1 > col0 else -1 - a = row0 - b = col0 - - b_coords = numpy.empty((da + 1, width), dtype=numpy.int32) - a_coords = numpy.empty((da + 1, width), dtype=numpy.int32) - - with nogil: - b -= (width - 1) // 2 - delta = 2 * db - da - for index in range(da + 1): - for offset in range(width): - b_coords[index, offset] = b + offset - a_coords[index, offset] = a - - if delta >= 0: # M2: Move by step_a + step_b - b += step_b - delta -= 2 * da - # else M1: Move by step_a - - a += step_a - delta += 2 * db - - if not invert_coords: - return (numpy.asarray(b_coords).reshape(-1), - numpy.asarray(a_coords).reshape(-1)) - else: - return (numpy.asarray(a_coords).reshape(-1), - numpy.asarray(b_coords).reshape(-1)) - - -def circle_fill(int crow, int ccol, float radius): - """Generates coordinate of image points lying in a disk. - - :param int crow: Row of the center of the disk - :param int ccol: Column of the center of the disk - :param float radius: Radius of the disk - :return: Array coordinates of points inside the disk (might be negative) - :rtype: 2-tuple of numpy.ndarray (rows, cols) - """ - cdef int i_radius, len_coords - - radius = fabs(radius) - i_radius = <int>radius - - coords = numpy.arange(-i_radius, ceil(radius) + 1, - dtype=numpy.float32) ** 2 - len_coords = len(coords) - # rows, cols = where(row**2 + col**2 < radius**2) - rows, cols = numpy.where(coords.reshape(1, len_coords) + - coords.reshape(len_coords, 1) < radius ** 2) - return rows + crow - i_radius, cols + ccol - i_radius - - -def ellipse_fill(int crow, int ccol, float radius_r, float radius_c): - """Generates coordinate of image points lying in a ellipse. - - :param int crow: Row of the center of the ellipse - :param int ccol: Column of the center of the ellipse - :param float radius_r: Radius of the ellipse in the row - :param float radius_c: Radius of the ellipse in the column - :return: Array coordinates of points inside the ellipse (might be negative) - :rtype: 2-tuple of numpy.ndarray (rows, cols) - """ - cdef int i_radius_r - cdef int i_radius_c - - i_radius_r = <int>fabs(radius_r) - i_radius_c = <int>fabs(radius_c) - - x_coords = numpy.arange(-i_radius_r, ceil(radius_r) + 1, dtype=numpy.float32).reshape(-1, 1) - y_coords = numpy.arange(-i_radius_c, ceil(radius_c) + 1, dtype=numpy.float32).reshape(1, -1) - - # rows, cols = where(x**2 + col**2 < 1) - rows, cols = numpy.where(x_coords**2 / radius_r**2 + y_coords**2 / radius_c**2 < 1.0) - return rows + crow - i_radius_r, cols + ccol - i_radius_c diff --git a/silx/image/sift.py b/silx/image/sift.py deleted file mode 100644 index cb1e6bd..0000000 --- a/silx/image/sift.py +++ /dev/null @@ -1,25 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2017 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ - -from silx.opencl.sift import * diff --git a/silx/image/test/__init__.py b/silx/image/test/__init__.py deleted file mode 100644 index f469edc..0000000 --- a/silx/image/test/__init__.py +++ /dev/null @@ -1,48 +0,0 @@ -# -*- coding: utf-8 -*- -# -# Project: silx -# https://github.com/silx-kit/silx -# -# Copyright (C) 2012-2018 European Synchrotron Radiation Facility, Grenoble, France -# -# 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. - -__authors__ = ["J. Kieffer"] -__license__ = "MIT" -__date__ = "17/04/2018" - -import unittest -from . import test_bilinear -from . import test_shapes -from . import test_medianfilter -from . import test_tomography -from . import test_bb -from ..marchingsquares.test import suite as marchingsquares_suite - - -def suite(): - """Test suite for module silx.image.test""" - test_suite = unittest.TestSuite() - test_suite.addTest(test_bilinear.suite()) - test_suite.addTest(test_medianfilter.suite()) - test_suite.addTest(test_shapes.suite()) - test_suite.addTest(test_tomography.suite()) - test_suite.addTest(marchingsquares_suite()) - test_suite.addTest(test_bb.suite()) - return test_suite diff --git a/silx/image/test/test_bb.py b/silx/image/test/test_bb.py deleted file mode 100644 index 3f33e80..0000000 --- a/silx/image/test/test_bb.py +++ /dev/null @@ -1,86 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# -# Copyright (c) 2016 European Synchrotron Radiation Facility -# -# 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. -# -# ###########################################################################*/ -"""Basic tests for Bounding box""" - -__authors__ = ["H. Payno"] -__license__ = "MIT" -__date__ = "27/09/2019" - - -import unittest -import numpy -from silx.image._boundingbox import _BoundingBox - - -class TestBB(unittest.TestCase): - """Some simple test on the bounding box class""" - def test_creation(self): - """test some constructors""" - pts = numpy.array([(0, 0), (10, 20), (20, 0)]) - bb = _BoundingBox.from_points(pts) - self.assertTrue(bb.bottom_left == (0, 0)) - self.assertTrue(bb.top_right == (20, 20)) - pts = numpy.array([(0, 10), (10, 20), (45, 30), (35, 0)]) - bb = _BoundingBox.from_points(pts) - self.assertTrue(bb.bottom_left == (0, 0)) - print(bb.top_right) - self.assertTrue(bb.top_right == (45, 30)) - - def test_isIn_pt(self): - """test the isIn function with points""" - bb = _BoundingBox(bottom_left=(6, 2), top_right=(12, 6)) - self.assertTrue(bb.contains((10, 4))) - self.assertTrue(bb.contains((6, 2))) - self.assertTrue(bb.contains((12, 2))) - self.assertFalse(bb.contains((0, 0))) - self.assertFalse(bb.contains((20, 0))) - self.assertFalse(bb.contains((10, 0))) - - def test_collide(self): - """test the collide function""" - bb1 = _BoundingBox(bottom_left=(6, 2), top_right=(12, 6)) - self.assertTrue(bb1.collide(_BoundingBox(bottom_left=(6, 2), top_right=(12, 6)))) - bb1 = _BoundingBox(bottom_left=(6, 2), top_right=(12, 6)) - self.assertFalse(bb1.collide(_BoundingBox(bottom_left=(12, 2), top_right=(12, 2)))) - - def test_isIn_bb(self): - """test the isIn function with other bounding box""" - bb1 = _BoundingBox(bottom_left=(6, 2), top_right=(12, 6)) - self.assertTrue(bb1.contains(_BoundingBox(bottom_left=(6, 2), top_right=(12, 6)))) - bb1 = _BoundingBox(bottom_left=(6, 2), top_right=(12, 6)) - self.assertTrue(bb1.contains(_BoundingBox(bottom_left=(12, 2), top_right=(12, 2)))) - self.assertFalse(_BoundingBox(bottom_left=(12, 2), top_right=(12, 2)).contains(bb1)) - - -def suite(): - test_suite = unittest.TestSuite() - for TestClass in (TestBB,): - test_suite.addTest( - unittest.defaultTestLoader.loadTestsFromTestCase(TestClass)) - return test_suite - - -if __name__ == '__main__': - unittest.main(defaultTest='suite') diff --git a/silx/image/test/test_bilinear.py b/silx/image/test/test_bilinear.py deleted file mode 100644 index 55eaccb..0000000 --- a/silx/image/test/test_bilinear.py +++ /dev/null @@ -1,178 +0,0 @@ -# -*- coding: utf-8 -*- -# -# Project: silx (originally pyFAI) -# https://github.com/silx-kit/silx -# -# Copyright (C) 2012-2017 European Synchrotron Radiation Facility, Grenoble, France -# 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. - -__authors__ = ["J. Kieffer"] -__license__ = "MIT" -__date__ = "25/11/2020" - -import unittest -import numpy -import logging -logger = logging.getLogger(__name__) -from ..bilinear import BilinearImage - - -class TestBilinear(unittest.TestCase): - """basic maximum search test""" - N = 1000 - - def test_max_search_round(self): - """test maximum search using random points: maximum is at the pixel center""" - a = numpy.arange(100) - 40. - b = numpy.arange(100) - 60. - ga = numpy.exp(-a * a / 4000) - gb = numpy.exp(-b * b / 6000) - gg = numpy.outer(ga, gb) - b = BilinearImage(gg) - - self.assertAlmostEqual(b.maxi, 1, 2, "maxi is almost 1") - self.assertLess(b.mini, 0.3, "mini should be around 0.23") - - ok = 0 - for s in range(self.N): - i, j = numpy.random.randint(100), numpy.random.randint(100) - k, l = b.local_maxi((i, j)) - if abs(k - 40) > 1e-4 or abs(l - 60) > 1e-4: - logger.warning("Wrong guess maximum (%i,%i) -> (%.1f,%.1f)", i, j, k, l) - else: - logger.debug("Good guess maximum (%i,%i) -> (%.1f,%.1f)", i, j, k, l) - ok += 1 - logger.debug("Success rate: %.1f", 100. * ok / self.N) - self.assertEqual(ok, self.N, "Maximum is always found") - - def test_max_search_half(self): - """test maximum search using random points: maximum is at a pixel edge""" - a = numpy.arange(100) - 40.5 - b = numpy.arange(100) - 60.5 - ga = numpy.exp(-a * a / 4000) - gb = numpy.exp(-b * b / 6000) - gg = numpy.outer(ga, gb) - b = BilinearImage(gg) - ok = 0 - for s in range(self.N): - i, j = numpy.random.randint(100), numpy.random.randint(100) - k, l = b.local_maxi((i, j)) - if abs(k - 40.5) > 0.5 or abs(l - 60.5) > 0.5: - logger.warning("Wrong guess maximum (%i,%i) -> (%.1f,%.1f)", i, j, k, l) - else: - logger.debug("Good guess maximum (%i,%i) -> (%.1f,%.1f)", i, j, k, l) - ok += 1 - logger.debug("Success rate: %.1f", 100. * ok / self.N) - self.assertEqual(ok, self.N, "Maximum is always found") - - def test_map(self): - N = 6 - y, x = numpy.ogrid[:N,:N + 10] - img = x + y - b = BilinearImage(img) - x2d = numpy.zeros_like(y) + x - y2d = numpy.zeros_like(x) + y - res1 = b.map_coordinates((y2d, x2d)) - self.assertEqual(abs(res1 - img).max(), 0, "images are the same (corners)") - - x2d = numpy.zeros_like(y) + (x[:,:-1] + 0.5) - y2d = numpy.zeros_like(x[:,:-1]) + y - res1 = b.map_coordinates((y2d, x2d)) - self.assertEqual(abs(res1 - img[:,:-1] - 0.5).max(), 0, "images are the same (middle)") - - x2d = numpy.zeros_like(y[:-1,:]) + (x[:,:-1] + 0.5) - y2d = numpy.zeros_like(x[:,:-1]) + (y[:-1,:] + 0.5) - res1 = b.map_coordinates((y2d, x2d)) - self.assertEqual(abs(res1 - img[:-1, 1:]).max(), 0, "images are the same (center)") - - def test_mask_grad(self): - N = 100 - img = numpy.arange(N * N).reshape(N, N) - # No mask on the boundaries, makes the test complicated, pixel always separated - masked = 2 * numpy.random.randint(0, int((N - 1) / 2), size=(2, N)) + 1 - mask = numpy.zeros((N, N), dtype=numpy.uint8) - mask[(masked[0], masked[1])] = 1 - self.assertLessEqual(mask.sum(), N, "At most N pixels are masked") - - b = BilinearImage(img, mask=mask) - self.assertEqual(b.has_mask, True, "interpolator has mask") - self.assertEqual(b.maxi, N * N - 1, "maxi is N²-1") - self.assertEqual(b.mini, 0, "mini is 0") - - y, x = numpy.ogrid[:N,:N] - x2d = numpy.zeros_like(y) + x - y2d = numpy.zeros_like(x) + y - res1 = b.map_coordinates((y2d, x2d)) - self.assertEqual(numpy.nanmax(abs(res1 - img)), 0, "images are the same (corners), or Nan ") - - x2d = numpy.zeros_like(y) + (x[:,:-1] + 0.5) - y2d = numpy.zeros_like(x[:,:-1]) + y - res1 = b.map_coordinates((y2d, x2d)) - self.assertLessEqual(numpy.max(abs(res1 - img[:, 1:] + 1 / 2.)), 0.5, "images are the same (middle) +/- 0.5") - - x2d = numpy.zeros_like(y[:-1]) + (x[:,:-1] + 0.5) - y2d = numpy.zeros_like(x[:,:-1]) + (y[:-1] + 0.5) - res1 = b.map_coordinates((y2d, x2d)) - exp = 0.25 * (img[:-1,:-1] + img[:-1, 1:] + img[1:,:-1] + img[1:, 1:]) - self.assertLessEqual(abs(res1 - exp).max(), N / 4, "images are almost the same (center)") - - def test_profile_grad(self): - N = 100 - img = numpy.arange(N * N).reshape(N, N) - b = BilinearImage(img) - res1 = b.profile_line((0, 0), (N - 1, N - 1)) - l = numpy.ceil(numpy.sqrt(2) * N) - self.assertEqual(len(res1), l, "Profile has correct length") - self.assertLess((res1[:-2] - res1[1:-1]).std(), 1e-3, "profile is linear (excluding last point)") - - def test_profile_gaus(self): - N = 100 - x = numpy.arange(N) - N // 2.0 - g = numpy.exp(-x * x / (N * N)) - img = numpy.outer(g, g) - b = BilinearImage(img) - res_hor = b.profile_line((N // 2, 0), (N // 2, N - 1)) - res_ver = b.profile_line((0, N // 2), (N - 1, N // 2)) - self.assertEqual(len(res_hor), N, "Profile has correct length") - self.assertEqual(len(res_ver), N, "Profile has correct length") - self.assertLess(abs(res_hor - g).max(), 1e-5, "correct horizontal profile") - self.assertLess(abs(res_ver - g).max(), 1e-5, "correct vertical profile") - - # Profile with linewidth=3 - expected_profile = img[:, N // 2 - 1:N // 2 + 2].mean(axis=1) - res_hor = b.profile_line((N // 2, 0), (N // 2, N - 1), linewidth=3) - res_ver = b.profile_line((0, N // 2), (N - 1, N // 2), linewidth=3) - - self.assertEqual(len(res_hor), N, "Profile has correct length") - self.assertEqual(len(res_ver), N, "Profile has correct length") - self.assertLess(abs(res_hor - expected_profile).max(), 1e-5, - "correct horizontal profile") - self.assertLess(abs(res_ver - expected_profile).max(), 1e-5, - "correct vertical profile") - - -def suite(): - testsuite = unittest.TestSuite() - testsuite.addTest(TestBilinear("test_max_search_round")) - testsuite.addTest(TestBilinear("test_max_search_half")) - testsuite.addTest(TestBilinear("test_map")) - testsuite.addTest(TestBilinear("test_profile_grad")) - testsuite.addTest(TestBilinear("test_profile_gaus")) - testsuite.addTest(TestBilinear("test_mask_grad")) - return testsuite diff --git a/silx/image/test/test_medianfilter.py b/silx/image/test/test_medianfilter.py deleted file mode 100644 index 5b062d9..0000000 --- a/silx/image/test/test_medianfilter.py +++ /dev/null @@ -1,76 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# -# Copyright (c) 2016-2017 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ -"""Tests that the different implementation of opencl (cpp, opencl) are - accessible -""" - -__authors__ = ["H. Payno"] -__license__ = "MIT" -__date__ = "11/05/2017" - -import unittest -from silx.image import medianfilter -import numpy - -from silx.opencl.common import ocl - - -class TestMedianFilterEngines(unittest.TestCase): - """Make sure we have access to all the different implementation of - median filter from image medfilt""" - - - IMG = numpy.arange(10000.).reshape(100, 100) - - KERNEL = (1, 1) - - def testCppMedFilt2d(self): - """test cpp engine for medfilt2d""" - res = medianfilter.medfilt2d( - image=TestMedianFilterEngines.IMG, - kernel_size=TestMedianFilterEngines.KERNEL, - engine='cpp') - self.assertTrue(numpy.array_equal(res, TestMedianFilterEngines.IMG)) - - @unittest.skipUnless(ocl, "PyOpenCl is missing") - def testOpenCLMedFilt2d(self): - """test cpp engine for medfilt2d""" - res = medianfilter.medfilt2d( - image=TestMedianFilterEngines.IMG, - kernel_size=TestMedianFilterEngines.KERNEL, - engine='opencl') - self.assertTrue(numpy.array_equal(res, TestMedianFilterEngines.IMG)) - - -def suite(): - test_suite = unittest.TestSuite() - for testClass in (TestMedianFilterEngines, ): - test_suite.addTest( - unittest.defaultTestLoader.loadTestsFromTestCase(testClass)) - return test_suite - - -if __name__ == '__main__': - unittest.main(defaultTest='suite') diff --git a/silx/image/test/test_shapes.py b/silx/image/test/test_shapes.py deleted file mode 100644 index 6539bba..0000000 --- a/silx/image/test/test_shapes.py +++ /dev/null @@ -1,366 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# -# Copyright (c) 2016 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ -"""Tests for polygon functions -""" - -__authors__ = ["T. Vincent"] -__license__ = "MIT" -__date__ = "15/02/2019" - - -import logging -import unittest -import numpy - -from silx.utils.testutils import ParametricTestCase -from silx.image import shapes - -_logger = logging.getLogger(__name__) - - -class TestPolygonFill(ParametricTestCase): - """basic poylgon test""" - - def test_squares(self): - """Test polygon fill for a square polygons""" - mask_shape = 4, 4 - tests = { - # test name: [(row min, row max), (col min, col max)] - 'square in': [(1, 3), (1, 3)], - 'square out': [(1, 3), (1, 10)], - 'square around': [(-1, 5), (-1, 5)], - } - - for test_name, (rows, cols) in tests.items(): - with self.subTest(msg=test_name, rows=rows, cols=cols, - mask_shape=mask_shape): - ref_mask = numpy.zeros(mask_shape, dtype=numpy.uint8) - ref_mask[max(0, rows[0]):rows[1], - max(0, cols[0]):cols[1]] = True - - vertices = [(rows[0], cols[0]), (rows[1], cols[0]), - (rows[1], cols[1]), (rows[0], cols[1])] - mask = shapes.polygon_fill_mask(vertices, ref_mask.shape) - is_equal = numpy.all(numpy.equal(ref_mask, mask)) - if not is_equal: - _logger.debug('%s failed with mask != ref_mask:', - test_name) - _logger.debug('result:\n%s', str(mask)) - _logger.debug('ref:\n%s', str(ref_mask)) - self.assertTrue(is_equal) - - def test_eight(self): - """Tests with eight shape with different rotation and direction""" - ref_mask = numpy.array(( - (1, 1, 1, 1, 1, 0), - (0, 1, 1, 1, 0, 0), - (0, 0, 1, 0, 0, 0), - (0, 0, 1, 0, 0, 0), - (0, 1, 1, 1, 0, 0), - (0, 0, 0, 0, 0, 0)), dtype=numpy.uint8) - ref_mask_rot = numpy.asarray(numpy.logical_not(ref_mask), - dtype=numpy.uint8) - ref_mask_rot[:, -1] = 0 - ref_mask_rot[-1, :] = 0 - - tests = { - 'dir 1': ([(0, 0), (5, 5), (5, 0), (0, 5)], ref_mask), - 'dir 1, rot 90': ([(5, 0), (0, 5), (5, 5), (0, 0)], ref_mask_rot), - 'dir 1, rot 180': ([(5, 5), (0, 0), (0, 5), (5, 0)], ref_mask), - 'dir 1, rot -90': ([(0, 5), (5, 0), (0, 0), (5, 5)], ref_mask_rot), - 'dir 2': ([(0, 0), (0, 5), (5, 0), (5, 5)], ref_mask), - 'dir 2, rot 90': ([(5, 0), (0, 0), (5, 5), (0, 5)], ref_mask_rot), - 'dir 2, rot 180': ([(5, 5), (5, 0), (0, 5), (0, 0)], ref_mask), - 'dir 2, rot -90': ([(0, 5), (5, 5), (0, 0), (5, 0)], ref_mask_rot), - } - - for test_name, (vertices, ref_mask) in tests.items(): - with self.subTest(msg=test_name): - mask = shapes.polygon_fill_mask(vertices, ref_mask.shape) - is_equal = numpy.all(numpy.equal(ref_mask, mask)) - if not is_equal: - _logger.debug('%s failed with mask != ref_mask:', - test_name) - _logger.debug('result:\n%s', str(mask)) - _logger.debug('ref:\n%s', str(ref_mask)) - self.assertTrue(is_equal) - - def test_shapes(self): - """Tests with shapes and reference mask""" - tests = { - # name: ( - # polygon corners as a list of (row, col), - # ref_mask) - 'concave polygon': ( - [(1, 1), (4, 3), (1, 5), (2, 3)], - numpy.array(( - (0, 0, 0, 0, 0, 0, 0, 0), - (0, 0, 0, 0, 0, 0, 0, 0), - (0, 0, 1, 1, 1, 0, 0, 0), - (0, 0, 0, 1, 0, 0, 0, 0), - (0, 0, 0, 0, 0, 0, 0, 0), - (0, 0, 0, 0, 0, 0, 0, 0)), dtype=numpy.uint8)), - 'concave polygon partly outside mask': ( - [(-1, -1), (4, 3), (1, 5), (2, 3)], - numpy.array(( - (1, 0, 0, 0, 0, 0), - (0, 1, 0, 0, 0, 0), - (0, 0, 1, 1, 1, 0), - (0, 0, 0, 1, 0, 0), - (0, 0, 0, 0, 0, 0), - (0, 0, 0, 0, 0, 0), - (0, 0, 0, 0, 0, 0), - (0, 0, 0, 0, 0, 0)), dtype=numpy.uint8)), - 'polygon surrounding mask': ( - [(-1, -1), (-1, 7), (7, 7), (7, -1), (0, -1), - (8, -2), (8, 8), (-2, 8)], - numpy.zeros((6, 6), dtype=numpy.uint8)) - } - - for test_name, (vertices, ref_mask) in tests.items(): - with self.subTest(msg=test_name): - mask = shapes.polygon_fill_mask(vertices, ref_mask.shape) - is_equal = numpy.all(numpy.equal(ref_mask, mask)) - if not is_equal: - _logger.debug('%s failed with mask != ref_mask:', - test_name) - _logger.debug('result:\n%s', str(mask)) - _logger.debug('ref:\n%s', str(ref_mask)) - self.assertTrue(is_equal) - - -class TestDrawLine(ParametricTestCase): - """basic draw line test""" - - def test_aligned_lines(self): - """Test drawing horizontal, vertical and diagonal lines""" - - lines = { # test_name: (drow, dcol) - 'Horizontal line, col0 < col1': (0, 10), - 'Horizontal line, col0 > col1': (0, -10), - 'Vertical line, row0 < row1': (10, 0), - 'Vertical line, row0 > row1': (-10, 0), - 'Diagonal col0 < col1 and row0 < row1': (10, 10), - 'Diagonal col0 < col1 and row0 > row1': (-10, 10), - 'Diagonal col0 > col1 and row0 < row1': (10, -10), - 'Diagonal col0 > col1 and row0 > row1': (-10, -10), - } - row0, col0 = 1, 2 # Start point - - for test_name, (drow, dcol) in lines.items(): - row1 = row0 + drow - col1 = col0 + dcol - with self.subTest(msg=test_name, drow=drow, dcol=dcol): - # Build reference coordinates from drow and dcol - if drow == 0: - rows = row0 * numpy.ones(abs(dcol) + 1) - else: - step = 1 if drow > 0 else -1 - rows = row0 + numpy.arange(0, drow + step, step) - - if dcol == 0: - cols = col0 * numpy.ones(abs(drow) + 1) - else: - step = 1 if dcol > 0 else -1 - cols = col0 + numpy.arange(0, dcol + step, step) - ref_coords = rows, cols - - result = shapes.draw_line(row0, col0, row1, col1) - self.assertTrue(self.isEqual(test_name, result, ref_coords)) - - def test_noline(self): - """Test pt0 == pt1""" - for width in range(4): - with self.subTest(width=width): - result = shapes.draw_line(1, 2, 1, 2, width) - self.assertTrue(numpy.all(numpy.equal(result, [(1,), (2,)]))) - - def test_lines(self): - """Test lines not aligned with axes for 8 slopes and directions""" - row0, col0 = 1, 1 - - dy, dx = 3, 5 - ref_coords = numpy.array( - [(0, 0), (1, 1), (1, 2), (2, 3), (2, 4), (3, 5)]) - - # Build lines for the 8 octants from this coordinantes - lines = { # name: (drow, dcol, ref_coords) - '1st octant': (dy, dx, ref_coords), - '2nd octant': (dx, dy, ref_coords[:, (1, 0)]), # invert x and y - '3rd octant': (dx, -dy, ref_coords[:, (1, 0)] * (1, -1)), - '4th octant': (dy, -dx, ref_coords * (1, -1)), - '5th octant': (-dy, -dx, ref_coords * (-1, -1)), - '6th octant': (-dx, -dy, ref_coords[:, (1, 0)] * (-1, -1)), - '7th octant': (-dx, dy, ref_coords[:, (1, 0)] * (-1, 1)), - '8th octant': (-dy, dx, ref_coords * (-1, 1)) - } - - # Test with different starting points with positive and negative coords - for row0, col0 in ((0, 0), (2, 3), (-4, 1), (-5, -6), (8, -7)): - for name, (drow, dcol, ref_coords) in lines.items(): - row1 = row0 + drow - col1 = col0 + dcol - # Transpose from ((row0, col0), ...) to (rows, cols) - ref_coords = numpy.transpose(ref_coords + (row0, col0)) - - with self.subTest(msg=name, - pt0=(row0, col0), pt1=(row1, col1)): - result = shapes.draw_line(row0, col0, row1, col1) - self.assertTrue(self.isEqual(name, result, ref_coords)) - - def test_width(self): - """Test of line width""" - - lines = { # test_name: row0, col0, row1, col1, width, ref - 'horizontal w=2': - (0, 0, 0, 1, 2, ((0, 1, 0, 1), - (0, 0, 1, 1))), - 'horizontal w=3': - (0, 0, 0, 1, 3, ((-1, 0, 1, -1, 0, 1), - (0, 0, 0, 1, 1, 1))), - 'vertical w=2': - (0, 0, 1, 0, 2, ((0, 0, 1, 1), - (0, 1, 0, 1))), - 'vertical w=3': - (0, 0, 1, 0, 3, ((0, 0, 0, 1, 1, 1), - (-1, 0, 1, -1, 0, 1))), - 'diagonal w=3': - (0, 0, 1, 1, 3, ((-1, 0, 1, 0, 1, 2), - (0, 0, 0, 1, 1, 1))), - '1st octant w=3': - (0, 0, 1, 2, 3, - numpy.array(((-1, 0), (0, 0), (1, 0), - (0, 1), (1, 1), (2, 1), - (0, 2), (1, 2), (2, 2))).T), - '2nd octant w=3': - (0, 0, 2, 1, 3, - numpy.array(((0, -1), (0, 0), (0, 1), - (1, 0), (1, 1), (1, 2), - (2, 0), (2, 1), (2, 2))).T), - } - - for test_name, (row0, col0, row1, col1, width, ref) in lines.items(): - with self.subTest(msg=test_name, - pt0=(row0, col0), pt1=(row1, col1), width=width): - result = shapes.draw_line(row0, col0, row1, col1, width) - self.assertTrue(self.isEqual(test_name, result, ref)) - - def isEqual(self, test_name, result, ref): - """Test equality of two numpy arrays and log them if different""" - is_equal = numpy.all(numpy.equal(result, ref)) - if not is_equal: - _logger.debug('%s failed with result != ref:', - test_name) - _logger.debug('result:\n%s', str(result)) - _logger.debug('ref:\n%s', str(ref)) - return is_equal - - -class TestCircleFill(ParametricTestCase): - """Tests for circle filling""" - - def testCircle(self): - """Test circle_fill with different input parameters""" - - square3x3 = numpy.array(((-1, -1, -1, 0, 0, 0, 1, 1, 1), - (-1, 0, 1, -1, 0, 1, -1, 0, 1))) - - tests = [ - # crow, ccol, radius, ref_coords = (ref_rows, ref_cols) - (0, 0, 1, ((0,), (0,))), - (10, 15, 1, ((10,), (15,))), - (0, 0, 1.5, square3x3), - (5, 10, 2, (5 + square3x3[0], 10 + square3x3[1])), - (10, 20, 3.5, ( - 10 + numpy.array((-3, -3, -3, - -2, -2, -2, -2, -2, - -1, -1, -1, -1, -1, -1, -1, - 0, 0, 0, 0, 0, 0, 0, - 1, 1, 1, 1, 1, 1, 1, - 2, 2, 2, 2, 2, - 3, 3, 3)), - 20 + numpy.array((-1, 0, 1, - -2, -1, 0, 1, 2, - -3, -2, -1, 0, 1, 2, 3, - -3, -2, -1, 0, 1, 2, 3, - -3, -2, -1, 0, 1, 2, 3, - -2, -1, 0, 1, 2, - -1, 0, 1)))), - ] - - for crow, ccol, radius, ref_coords in tests: - with self.subTest(crow=crow, ccol=ccol, radius=radius): - coords = shapes.circle_fill(crow, ccol, radius) - is_equal = numpy.all(numpy.equal(coords, ref_coords)) - if not is_equal: - _logger.debug('result:\n%s', str(coords)) - _logger.debug('ref:\n%s', str(ref_coords)) - self.assertTrue(is_equal) - - -class TestEllipseFill(unittest.TestCase): - """Tests for ellipse filling""" - - def testPoint(self): - args = [1, 1, 1, 1] - result = shapes.ellipse_fill(*args) - expected = numpy.array(([1], [1])) - numpy.testing.assert_array_equal(result, expected) - - def testTranslatedPoint(self): - args = [10, 10, 1, 1] - result = shapes.ellipse_fill(*args) - expected = numpy.array(([10], [10])) - numpy.testing.assert_array_equal(result, expected) - - def testEllipse(self): - args = [0, 0, 20, 10] - rows, cols = shapes.ellipse_fill(*args) - self.assertEqual(len(rows), 617) - self.assertEqual(rows.mean(), 0) - self.assertAlmostEqual(rows.std(), 10.025575, places=3) - self.assertEqual(len(cols), 617) - self.assertEqual(cols.mean(), 0) - self.assertAlmostEqual(cols.std(), 4.897325, places=3) - - def testTranslatedEllipse(self): - args = [0, 0, 20, 10] - expected_rows, expected_cols = shapes.ellipse_fill(*args) - args = [10, 50, 20, 10] - rows, cols = shapes.ellipse_fill(*args) - numpy.testing.assert_allclose(rows, expected_rows + 10) - numpy.testing.assert_allclose(cols, expected_cols + 50) - - -def suite(): - test_suite = unittest.TestSuite() - for testClass in (TestPolygonFill, TestDrawLine, TestCircleFill, TestEllipseFill): - test_suite.addTest( - unittest.defaultTestLoader.loadTestsFromTestCase(testClass)) - return test_suite - - -if __name__ == '__main__': - unittest.main(defaultTest='suite') diff --git a/silx/image/test/test_tomography.py b/silx/image/test/test_tomography.py deleted file mode 100644 index 2a6a33c..0000000 --- a/silx/image/test/test_tomography.py +++ /dev/null @@ -1,66 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# -# Copyright (c) 2016-2017 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ -""" -Tests that the functions of tomography are valid -""" - -__authors__ = ["H. Payno"] -__license__ = "MIT" -__date__ = "12/09/2017" - -import unittest -import numpy -from silx.test.utils import utilstest -from silx.image import tomography - -class TestTomography(unittest.TestCase): - """ - - """ - - def setUp(self): - self.sinoTrueData = numpy.load(utilstest.getfile("sino500.npz"))["data"] - - def testCalcCenterCentroid(self): - centerTD = tomography.calc_center_centroid(self.sinoTrueData) - self.assertTrue(numpy.isclose(centerTD, 256, rtol=0.01)) - - def testCalcCenterCorr(self): - centerTrueData = tomography.calc_center_corr(self.sinoTrueData, - fullrot=False, - props=1) - self.assertTrue(numpy.isclose(centerTrueData, 256, rtol=0.01)) - - -def suite(): - test_suite = unittest.TestSuite() - for testClass in (TestTomography, ): - test_suite.addTest( - unittest.defaultTestLoader.loadTestsFromTestCase(testClass)) - return test_suite - - -if __name__ == '__main__': - unittest.main(defaultTest='suite') diff --git a/silx/image/tomography.py b/silx/image/tomography.py deleted file mode 100644 index 53855c1..0000000 --- a/silx/image/tomography.py +++ /dev/null @@ -1,301 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# Copyright (C) 2017 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ -""" -This module contains utilitary functions for tomography -""" - -__author__ = ["P. Paleo"] -__license__ = "MIT" -__date__ = "12/09/2017" - - -import numpy as np -from math import pi -from functools import lru_cache -from itertools import product -from bisect import bisect -from silx.math.fit import leastsq - -# ------------------------------------------------------------------------------ -# -------------------- Filtering-related functions ----------------------------- -# ------------------------------------------------------------------------------ - -def compute_ramlak_filter(dwidth_padded, dtype=np.float32): - """ - Compute the Ramachandran-Lakshminarayanan (Ram-Lak) filter, used in - filtered backprojection. - - :param dwidth_padded: width of the 2D sinogram after padding - :param dtype: data type - """ - L = dwidth_padded - h = np.zeros(L, dtype=dtype) - L2 = L//2+1 - h[0] = 1/4. - j = np.linspace(1, L2, L2//2, False).astype(dtype) # np < 1.9.0 - h[1:L2:2] = -1./(pi**2 * j**2) - h[L2:] = np.copy(h[1:L2-1][::-1]) - return h - - -def tukey(N, alpha=0.5): - """ - Compute the Tukey apodization window. - - :param int N: Number of points. - :param float alpha: - """ - apod = np.zeros(N) - x = np.arange(N)/(N-1) - r = alpha - M1 = (0 <= x) * (x < r/2) - M2 = (r/2 <= x) * (x <= 1 - r/2) - M3 = (1 - r/2 < x) * (x <= 1) - apod[M1] = (1 + np.cos(2*pi/r * (x[M1] - r/2)))/2. - apod[M2] = 1. - apod[M3] = (1 + np.cos(2*pi/r * (x[M3] - 1 + r/2)))/2. - return apod - - -def lanczos(N): - """ - Compute the Lanczos window (truncated sinc) of width N. - - :param int N: window width - """ - x = np.arange(N)/(N-1) - return np.sin(pi*(2*x-1))/(pi*(2*x-1)) - - -def compute_fourier_filter(dwidth_padded, filter_name, cutoff=1.): - """ - Compute the filter used for FBP. - - :param dwidth_padded: padded detector width. As the filtering is done by the - Fourier convolution theorem, dwidth_padded should be at least 2*dwidth. - :param filter_name: Name of the filter. Available filters are: - Ram-Lak, Shepp-Logan, Cosine, Hamming, Hann, Tukey, Lanczos. - :param cutoff: Cut-off frequency, if relevant. - """ - Nf = dwidth_padded - #~ filt_f = np.abs(np.fft.fftfreq(Nf)) - rl = compute_ramlak_filter(Nf, dtype=np.float64) - filt_f = np.fft.fft(rl) - - filter_name = filter_name.lower() - if filter_name in ["ram-lak", "ramlak"]: - return filt_f - - w = 2 * pi * np.fft.fftfreq(dwidth_padded) - d = cutoff - apodization = { - # ~OK - "shepp-logan": np.sin(w[1:Nf]/(2*d))/(w[1:Nf]/(2*d)), - # ~OK - "cosine": np.cos(w[1:Nf]/(2*d)), - # OK - "hamming": 0.54*np.ones_like(filt_f)[1:Nf] + .46 * np.cos(w[1:Nf]/d), - # OK - "hann": (np.ones_like(filt_f)[1:Nf] + np.cos(w[1:Nf]/d))/2., - # These one is not compatible with Astra - TODO investigate why - "tukey": np.fft.fftshift(tukey(dwidth_padded, alpha=d/2.))[1:Nf], - "lanczos": np.fft.fftshift(lanczos(dwidth_padded))[1:Nf], - } - if filter_name not in apodization: - raise ValueError("Unknown filter %s. Available filters are %s" % - (filter_name, str(apodization.keys()))) - filt_f[1:Nf] *= apodization[filter_name] - return filt_f - - -@lru_cache(maxsize=1) -def generate_powers(): - """ - Generate a list of powers of [2, 3, 5, 7], - up to (2**15)*(3**9)*(5**6)*(7**5). - """ - primes = [2, 3, 5, 7] - maxpow = {2: 15, 3: 9, 5: 6, 7: 5} - valuations = [] - for prime in primes: - # disallow any odd number (for R2C transform), and any number - # not multiple of 4 (Ram-Lak filter behaves strangely when - # dwidth_padded/2 is not even) - minval = 2 if prime == 2 else 0 - valuations.append(range(minval, maxpow[prime]+1)) - powers = product(*valuations) - res = [] - for pw in powers: - res.append(np.prod(list(map(lambda x : x[0]**x[1], zip(primes, pw))))) - return np.unique(res) - - -def get_next_power(n, powers=None): - """ - Given a number, get the closest (upper) number p such that - p is a power of 2, 3, 5 and 7. - """ - if powers is None: - powers = generate_powers() - idx = bisect(powers, n) - if powers[idx-1] == n: - return n - return powers[idx] - - -# ------------------------------------------------------------------------------ -# ------------- Functions for determining the center of rotation -------------- -# ------------------------------------------------------------------------------ - - - -def calc_center_corr(sino, fullrot=False, props=1): - """ - Compute a guess of the Center of Rotation (CoR) of a given sinogram. - The computation is based on the correlation between the line projections at - angle (theta = 0) and at angle (theta = 180). - - Note that for most scans, the (theta=180) angle is not included, - so the CoR might be underestimated. - In a [0, 360[ scan, the projection angle at (theta=180) is exactly in the - middle for odd number of projections. - - :param numpy.ndarray sino: Sinogram - :param bool fullrot: optional. If False (default), the scan is assumed to - be [0, 180). - If True, the scan is assumed to be [0, 380). - :param int props: optional. Number of propositions for the CoR - """ - - n_a, n_d = sino.shape - first = 0 - last = -1 if not(fullrot) else n_a // 2 - proj1 = sino[first, :] - proj2 = sino[last, :][::-1] - - # Compute the correlation in the Fourier domain - proj1_f = np.fft.fft(proj1, 2 * n_d) - proj2_f = np.fft.fft(proj2, 2 * n_d) - corr = np.abs(np.fft.ifft(proj1_f * proj2_f.conj())) - - if props == 1: - pos = np.argmax(corr) - if pos > n_d // 2: - pos -= n_d - return (n_d + pos) / 2. - else: - corr_argsorted = np.argsort(corr)[:props] - corr_argsorted[corr_argsorted > n_d // 2] -= n_d - return (n_d + corr_argsorted) / 2. - - -def _sine_function(t, offset, amplitude, phase): - """ - Helper function for calc_center_centroid - """ - n_angles = t.shape[0] - res = amplitude * np.sin(2 * pi * (1. / (2 * n_angles)) * t + phase) - return offset + res - - -def _sine_function_derivative(t, params, eval_idx): - """ - Helper function for calc_center_centroid - """ - offset, amplitude, phase = params - n_angles = t.shape[0] - w = 2.0 * pi * (1. / (2.0 * n_angles)) * t + phase - grad = (1.0, np.sin(w), amplitude*np.cos(w)) - return grad[eval_idx] - - -def calc_center_centroid(sino): - """ - Compute a guess of the Center of Rotation (CoR) of a given sinogram. - The computation is based on the computation of the centroid of each - projection line, which should be a sine function according to the - Helgason-Ludwig condition. - This method is unlikely to work in local tomography. - - :param numpy.ndarray sino: Sinogram - """ - - n_a, n_d = sino.shape - # Compute the vector of centroids of the sinogram - i = np.arange(n_d) - centroids = np.sum(sino*i, axis=1)/np.sum(sino, axis=1) - - # Fit with a sine function : phase, amplitude, offset - # Using non-linear Levenberg–Marquardt algorithm - angles = np.linspace(0, n_a, n_a, True) - # Initial parameter vector - cmax, cmin = centroids.max(), centroids.min() - offs = (cmax + cmin) / 2. - amp = (cmax - cmin) / 2. - phi = 1.1 - p0 = (offs, amp, phi) - - constraints = np.zeros((3, 3)) - - popt, _ = leastsq(model=_sine_function, - xdata=angles, - ydata=centroids, - p0=p0, - sigma=None, - constraints=constraints, - model_deriv=None, - epsfcn=None, - deltachi=None, - full_output=0, - check_finite=True, - left_derivative=False, - max_iter=100) - return popt[0] - - - -# ------------------------------------------------------------------------------ -# -------------------- Visualization-related functions ------------------------- -# ------------------------------------------------------------------------------ - - -def rescale_intensity(img, from_subimg=None, percentiles=None): - """ - clamp intensity into the [2, 98] percentiles - - :param img: - :param from_subimg: - :param percentiles: - :return: the rescale intensity - """ - if percentiles is None: - percentiles = [2, 98] - else: - assert type(percentiles) in (tuple, list) - assert(len(percentiles) == 2) - data = from_subimg if from_subimg is not None else img - imin, imax = np.percentile(data, percentiles) - res = np.clip(img, imin, imax) - return res - diff --git a/silx/image/utils.py b/silx/image/utils.py deleted file mode 100644 index 996d010..0000000 --- a/silx/image/utils.py +++ /dev/null @@ -1,53 +0,0 @@ -# -*- coding: utf-8 -*- -# /*########################################################################## -# Copyright (C) 2019 European Synchrotron Radiation Facility -# -# 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. -# -# ############################################################################*/ - -import numpy as np -from math import ceil - -def gaussian_kernel(sigma, cutoff=4, force_odd_size=False): - """ - Generates a Gaussian convolution kernel. - - :param sigma: Standard Deviation of the Gaussian curve. - :param cutoff: Parameter tuning the truncation of the Gaussian. - The higher cutoff, the biggest the array will be (and the closest to - a "true" Gaussian function). - :param force_odd_size: when set to True, the resulting array will always - have an odd size, regardless of the values of "sigma" and "cutoff". - :return: a numpy.ndarray containing the truncated Gaussian function. - The array size is 2*c*s+1 where c=cutoff, s=sigma. - - Nota: due to the quick decay of the Gaussian function, small values of the - "cutoff" parameter are usually fine. The energy difference between a - Gaussian truncated to [-c, c] and a "true" one is - erfc(c/(sqrt(2)*s)) - so choosing cutoff=4*sigma keeps the truncation error below 1e-4. - """ - size = int(ceil(2 * cutoff * sigma + 1)) - if force_odd_size and size % 2 == 0: - size += 1 - x = np.arange(size) - (size - 1.0) / 2.0 - g = np.exp(-(x / sigma) ** 2 / 2.0) - g /= g.sum() - return g |