diff options
Diffstat (limited to 'src/silx/image/marchingsquares')
| -rw-r--r-- | src/silx/image/marchingsquares/__init__.py | 119 | ||||
| -rw-r--r-- | src/silx/image/marchingsquares/_mergeimpl.pyx | 1324 | ||||
| -rw-r--r-- | src/silx/image/marchingsquares/_skimage.py | 138 | ||||
| -rw-r--r-- | src/silx/image/marchingsquares/include/patterns.h | 89 | ||||
| -rw-r--r-- | src/silx/image/marchingsquares/test/__init__.py | 23 | ||||
| -rw-r--r-- | src/silx/image/marchingsquares/test/test_funcapi.py | 90 | ||||
| -rw-r--r-- | src/silx/image/marchingsquares/test/test_mergeimpl.py | 261 |
7 files changed, 2044 insertions, 0 deletions
diff --git a/src/silx/image/marchingsquares/__init__.py b/src/silx/image/marchingsquares/__init__.py new file mode 100644 index 0000000..a310e70 --- /dev/null +++ b/src/silx/image/marchingsquares/__init__.py @@ -0,0 +1,119 @@ +# /*########################################################################## +# 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/src/silx/image/marchingsquares/_mergeimpl.pyx b/src/silx/image/marchingsquares/_mergeimpl.pyx new file mode 100644 index 0000000..84e53bb --- /dev/null +++ b/src/silx/image/marchingsquares/_mergeimpl.pyx @@ -0,0 +1,1324 @@ +#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 + +# /*########################################################################## +# Copyright (C) 2018-2023 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__ = "21/12/2023" + +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() noexcept 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() noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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) noexcept 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/src/silx/image/marchingsquares/_skimage.py b/src/silx/image/marchingsquares/_skimage.py new file mode 100644 index 0000000..2e136f7 --- /dev/null +++ b/src/silx/image/marchingsquares/_skimage.py @@ -0,0 +1,138 @@ +# /*########################################################################## +# 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/src/silx/image/marchingsquares/include/patterns.h b/src/silx/image/marchingsquares/include/patterns.h new file mode 100644 index 0000000..ff86cc1 --- /dev/null +++ b/src/silx/image/marchingsquares/include/patterns.h @@ -0,0 +1,89 @@ +# /*########################################################################## +# +# 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/src/silx/image/marchingsquares/test/__init__.py b/src/silx/image/marchingsquares/test/__init__.py new file mode 100644 index 0000000..0e66f61 --- /dev/null +++ b/src/silx/image/marchingsquares/test/__init__.py @@ -0,0 +1,23 @@ +# +# 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. diff --git a/src/silx/image/marchingsquares/test/test_funcapi.py b/src/silx/image/marchingsquares/test/test_funcapi.py new file mode 100644 index 0000000..0e5471c --- /dev/null +++ b/src/silx/image/marchingsquares/test/test_funcapi.py @@ -0,0 +1,90 @@ +# +# 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) diff --git a/src/silx/image/marchingsquares/test/test_mergeimpl.py b/src/silx/image/marchingsquares/test/test_mergeimpl.py new file mode 100644 index 0000000..db36b54 --- /dev/null +++ b/src/silx/image/marchingsquares/test/test_mergeimpl.py @@ -0,0 +1,261 @@ +# +# 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) |