diff options
Diffstat (limited to 'src/silx/image/marchingsquares/_mergeimpl.pyx')
| -rw-r--r-- | src/silx/image/marchingsquares/_mergeimpl.pyx | 1319 |
1 files changed, 1319 insertions, 0 deletions
diff --git a/src/silx/image/marchingsquares/_mergeimpl.pyx b/src/silx/image/marchingsquares/_mergeimpl.pyx new file mode 100644 index 0000000..5a7a3b5 --- /dev/null +++ b/src/silx/image/marchingsquares/_mergeimpl.pyx @@ -0,0 +1,1319 @@ +# 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 |