# coding: utf-8 # /*########################################################################## # Copyright (C) 2018 European Synchrotron Radiation Facility # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. # # ############################################################################*/ """ 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 include "../../utils/_have_openmp.pxi" """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 = 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 = 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 = libc.stdlib.malloc(context_size * sizeof(cnumpy.float32_t)) self._max_cache = 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