path: root/silx/image/marchingsquares/_mergeimpl.pyx

# 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 = <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