New upstream version 1.0.0+dfsg

27 files changed, 4526 insertions, 0 deletions

diff --git a/src/silx/image/__init__.py b/src/silx/image/__init__.py
new file mode 100644
index 0000000..12bf320
--- /dev/null
+++ b/src/silx/image/__init__.py
@@ -0,0 +1,33 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2017-2018 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+"""This package provides some processing functions for 2D images.
+
+For more processing functions, see the silx.math and silx.opencl packages.
+
+See silx documentation: http://www.silx.org/doc/silx/latest/
+"""
+
+__authors__ = ["H. Payno"]
+__license__ = "MIT"
+__date__ = "14/02/2017"
diff --git a/src/silx/image/_boundingbox.py b/src/silx/image/_boundingbox.py
new file mode 100644
index 0000000..1c086b1
--- /dev/null
+++ b/src/silx/image/_boundingbox.py
@@ -0,0 +1,100 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2018 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""offer some generic 2D bounding box features
+"""
+
+__authors__ = ["H. Payno"]
+__license__ = "MIT"
+__date__ = "07/09/2019"
+
+from silx.math.combo import min_max
+import numpy
+
+
+class _BoundingBox:
+    """
+    Simple 2D bounding box
+
+    :param tuple bottom_left: (y, x) bottom left point
+    :param tuple top_right: (y, x) top right point
+    """
+    def __init__(self, bottom_left, top_right):
+        self.bottom_left = bottom_left
+        self.top_right = top_right
+        self.min_x = bottom_left[1]
+        self.min_y = bottom_left[0]
+        self.max_x = top_right[1]
+        self.max_y = top_right[0]
+
+    def contains(self, item):
+        """
+
+        :param item: bouding box or point. If it is bounding check check if the
+                     bottom_left and top_right corner of the given bounding box
+                     are inside the current bounding box
+        :type: Union[BoundingBox, tuple]
+        :return:
+        """
+        if isinstance(item, _BoundingBox):
+            return self.contains(item.bottom_left) and self.contains(item.top_right)
+        else:
+            return (
+                    (self.min_x <= item[1] <= self.max_x) and
+                    (self.min_y <= item[0] <= self.max_y)
+            )
+
+    def collide(self, bb):
+        """
+        Check if the two bounding box collide
+
+        :param bb: bounding box to compare with
+        :type: :class:BoundingBox
+        :return: True if the two boxes collides
+        :rtype: bool
+        """
+        assert isinstance(bb, _BoundingBox)
+        return (
+            (self.min_x < bb.max_x and self.max_x > bb.min_x) and
+            (self.min_y < bb.max_y and self.max_y > bb.min_y)
+        )
+
+    @staticmethod
+    def from_points(points):
+        """
+
+        :param numpy.array tuple points: list of points. Should be 2D:
+                                         [(y1, x1), (y2, x2), (y3, x3), ...]
+        :return: bounding box from two points
+        :rtype: _BoundingBox
+        """
+        if not isinstance(points, numpy.ndarray):
+            points_ = numpy.ndarray(points)
+        else:
+            points_ = points
+        x = points_[:, 1]
+        y = points_[:, 0]
+        x_min, x_max = min_max(x)
+        y_min, y_max = min_max(y)
+        return _BoundingBox(bottom_left=(y_min, x_min), top_right=(y_max, x_max))
diff --git a/src/silx/image/backprojection.py b/src/silx/image/backprojection.py
new file mode 100644
index 0000000..63f99ca
--- /dev/null
+++ b/src/silx/image/backprojection.py
@@ -0,0 +1,25 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+
+from silx.opencl.backprojection import *
diff --git a/src/silx/image/bilinear.pyx b/src/silx/image/bilinear.pyx
new file mode 100644
index 0000000..14547f8
--- /dev/null
+++ b/src/silx/image/bilinear.pyx
@@ -0,0 +1,465 @@
+# -*- coding: utf-8 -*-
+#cython: embedsignature=True, language_level=3
+## This is for optimisation
+#cython: boundscheck=False, wraparound=False, cdivision=True, initializedcheck=False,
+## This is for developping:
+##cython: profile=True, warn.undeclared=True, warn.unused=True, warn.unused_result=False, warn.unused_arg=True
+#
+#
+#
+#    Project: silx (originally pyFAI)
+#             https://github.com/silx-kit/silx
+#
+#    Copyright (C) 2012-2020  European Synchrotron Radiation Facility, Grenoble, France
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+
+"""Bilinear interpolator, peak finder, line-profile for images"""
+__authors__ = ["J. Kieffer"]
+__license__ = "MIT"
+__date__ = "26/11/2020"
+
+# C-level imports
+from libc.stdint cimport uint8_t
+from libc.math cimport floor, ceil, sqrt, NAN, isfinite
+from libc.float cimport FLT_MAX
+
+import cython
+import numpy
+import logging
+logger = logging.getLogger(__name__)
+
+
+#Definition of some constants
+# How data are stored 
+ctypedef float data_t
+data_d = numpy.float32
+
+#How the mask is stored
+ctypedef uint8_t mask_t
+mask_d = numpy.uint8
+
+
+cdef class BilinearImage:
+    """Bilinear interpolator for images ... or any data on a regular grid
+    """
+    cdef:
+        readonly data_t[:, ::1] data
+        readonly mask_t[:, ::1] mask
+        readonly data_t maxi, mini
+        readonly Py_ssize_t width, height
+        readonly bint has_mask 
+    
+    # C-level declarations
+    cpdef Py_ssize_t coarse_local_maxi(self, Py_ssize_t)
+    cdef Py_ssize_t c_local_maxi(self, Py_ssize_t) nogil
+    cdef data_t c_funct(self, data_t, data_t) nogil
+    cdef void _init_min_max(self) nogil
+    
+    def __cinit__(self, data not None, mask=None):
+        """Constructor
+
+        :param data: image as a 2D array
+        """
+        assert data.ndim == 2
+        self.height = data.shape[0]
+        self.width = data.shape[1]
+        self.data = numpy.ascontiguousarray(data, dtype=data_d)
+        if mask is not None:
+            self.mask = numpy.ascontiguousarray(mask, dtype=mask_d)
+            self.has_mask=True
+        else:
+            self.mask = None
+            self.has_mask = False
+        self._init_min_max()
+
+    def __dealloc__(self):
+        self.data = None
+        self.mask = None
+
+    def __call__(self, coord):
+        """Function f((y, x)) where f is a continuous function
+        made from the image and (y,x)=(row, column) is the pixel coordinates
+        in natural C-order
+
+        :param x: 2-tuple of float (row, column)
+        :return: Interpolated signal from the image
+        """
+        return self.c_funct(coord[1], coord[0])
+    
+    cdef void _init_min_max(self) nogil:
+        "Calculate the min & max"
+        cdef:
+            Py_ssize_t i, j
+            data_t maxi, mini, value
+        mini = FLT_MAX
+        maxi = -FLT_MAX
+        for i in range(self.height):
+            for j in range(self.width):
+                if not (self.has_mask and self.mask[i,j]):
+                    value = self.data[i, j] 
+                    maxi = max(value, maxi)
+                    mini = min(value, mini)
+        self.maxi = maxi
+        self.mini = mini 
+
+    cdef data_t c_funct(self, data_t x, data_t y) nogil:
+        """Function f(x, y) where f is a continuous function
+        made from the image.
+
+        :param x (float): column coordinate
+        :param y (float): row coordinate
+        :return: Interpolated signal from the image (nearest for outside)
+
+        Cython only function due to NOGIL
+        """
+        cdef:
+            data_t d0 = min(max(y, 0.0), (self.height - 1.0))
+            data_t d1 = min(max(x, 0.0), (self.width - 1.0))
+            mask_t m0, m1, m2, m3
+            Py_ssize_t i0, i1, j0, j1
+            data_t x0, x1, y0, y1, res, scale
+
+        x0 = floor(d0)
+        x1 = ceil(d0)
+        y0 = floor(d1)
+        y1 = ceil(d1)
+        i0 = < int > x0
+        i1 = < int > x1
+        j0 = < int > y0
+        j1 = < int > y1
+        if (i0 == i1) and (j0 == j1):
+            if not (self.has_mask and self.mask[i0,j0]):
+                res = self.data[i0, j0]
+            else:
+                res = NAN
+        elif i0 == i1:
+            if self.has_mask:
+                m0 = self.mask[i0, j0]
+                m1 = self.mask[i0, j1]
+                if m0 and m1:
+                    res = NAN
+                elif m0:
+                    res = self.data[i0, j1]
+                elif m1:
+                    res = self.data[i0, j0]
+                else:
+                    res = (self.data[i0, j0] * (y1 - d1)) + (self.data[i0, j1] * (d1 - y0))                
+            else:
+                res = (self.data[i0, j0] * (y1 - d1)) + (self.data[i0, j1] * (d1 - y0))
+        elif j0 == j1:
+            if self.has_mask:
+                m0 = self.mask[i0, j0]
+                m1 = self.mask[i1, j0]
+                if m0 and m1:
+                    res = NAN
+                elif m0:
+                    res = self.data[i1, j0]
+                elif m1:
+                    res = self.data[i0, j0]
+                else:
+                    res = (self.data[i0, j0] * (x1 - d0)) + (self.data[i1, j0] * (d0 - x0))                
+            else:
+                res = (self.data[i0, j0] * (x1 - d0)) + (self.data[i1, j0] * (d0 - x0))
+        else:
+            if self.has_mask:
+                m0 = self.mask[i0, j0]
+                m1 = self.mask[i1, j0]
+                m2 = self.mask[i0, j1]
+                m3 = self.mask[i1, j1]
+                if m0 and m1 and m2 and m3:
+                    res = NAN
+                else:
+                    m0 = not m0
+                    m1 = not m1
+                    m2 = not m2
+                    m3 = not m3
+                    if m0 and m1 and m2 and m3:
+                        res = (self.data[i0, j0] * (x1 - d0) * (y1 - d1))  \
+                            + (self.data[i1, j0] * (d0 - x0) * (y1 - d1))  \
+                            + (self.data[i0, j1] * (x1 - d0) * (d1 - y0))  \
+                            + (self.data[i1, j1] * (d0 - x0) * (d1 - y0))                    
+                    else:
+                        res = (m0 * self.data[i0, j0] * (x1 - d0) * (y1 - d1))  \
+                            + (m1 * self.data[i1, j0] * (d0 - x0) * (y1 - d1))  \
+                            + (m2 * self.data[i0, j1] * (x1 - d0) * (d1 - y0))  \
+                            + (m3 * self.data[i1, j1] * (d0 - x0) * (d1 - y0))
+                        scale = ((m0 * (x1 - d0) * (y1 - d1)) 
+                               + (m1 * (d0 - x0) * (y1 - d1))
+                               + (m2 * (x1 - d0) * (d1 - y0))
+                               + (m3 * (d0 - x0) * (d1 - y0)))
+                        res /= scale
+            else:
+                res = (self.data[i0, j0] * (x1 - d0) * (y1 - d1))  \
+                    + (self.data[i1, j0] * (d0 - x0) * (y1 - d1))  \
+                    + (self.data[i0, j1] * (x1 - d0) * (d1 - y0))  \
+                    + (self.data[i1, j1] * (d0 - x0) * (d1 - y0))                    
+                
+        return res
+
+    def opp_f(self, coord):
+        """Function -f((y,x)) for peak finding via minimizer.
+
+        Gives large number outside the boundaries to return into the image
+
+        :param x: 2-tuple of float in natural C order, i.e (row, column)
+        :return: Negative interpolated signal from the image
+        """
+        cdef:
+            data_t d0, d1, res
+        d0, d1 = coord
+        if d0 < 0:
+            res = self.mini + d0
+        elif d1 < 0:
+            res = self.mini + d1
+        elif d0 > (self.height - 1):
+            res = self.mini - d0 + self.height - 1
+        elif d1 > self.width - 1:
+            res = self.mini - d1 + self.width - 1
+        else:
+            res = self.c_funct(d1, d0)
+        return - res
+
+    def local_maxi(self, coord):
+        """Return the nearest local maximum ... with sub-pixel refinement
+
+        Nearest maximum search:
+            steepest ascent
+
+        Sub-pixel refinement:
+            Second order Taylor expansion of the function;
+            At the maximum, the first derivative is null
+            delta = x-i = -Inverse[Hessian].gradient
+            if Hessian is singular or \|delta\|>1: use a center of mass.
+
+        :param coord: 2-tuple of scalar (row, column)
+        :return: 2-tuple of float with the nearest local maximum
+        """
+        cdef:
+            int res, current0, current1
+            int i0, i1
+            data_t tmp, sum0 = 0, sum1 = 0, sum = 0
+            data_t a00, a01, a02, a10, a11, a12, a20, a21, a22
+            data_t d00, d11, d01, denom, delta0, delta1
+        res = self.c_local_maxi(round(coord[0]) * self.width + round(coord[1]))
+
+        current0 = res // self.width
+        current1 = res % self.width
+        if (current0 > 0) and (current0 < self.height - 1) and (current1 > 0) and (current1 < self.width - 1):
+            # Use second order polynomial Taylor expansion
+            a00 = self.data[current0 - 1, current1 - 1]
+            a01 = self.data[current0 - 1, current1    ]
+            a02 = self.data[current0 - 1, current1 + 1]
+            a10 = self.data[current0    , current1 - 1]
+            a11 = self.data[current0    , current1    ]
+            a12 = self.data[current0    , current1 + 1]
+            a20 = self.data[current0 + 1, current1 - 1]
+            a21 = self.data[current0 + 1, current1    ]
+            a22 = self.data[current0 + 1, current1 - 1]
+            d00 = a12 - 2.0 * a11 + a10
+            d11 = a21 - 2.0 * a11 + a01
+            d01 = (a00 - a02 - a20 + a22) / 4.0
+            denom = 2.0 * (d00 * d11 - d01 * d01)
+            if abs(denom) < 1e-10:
+                logger.debug("Singular determinant, Hessian undefined")
+            else:
+                delta0 = ((a12 - a10) * d01 + (a01 - a21) * d11) / denom
+                delta1 = ((a10 - a12) * d00 + (a21 - a01) * d01) / denom
+                if abs(delta0) <= 1.0 and abs(delta1) <= 1.0:
+                    # Result is OK if lower than 0.5.
+                    return (delta0 + float(current0), delta1 + float(current1))
+                else:
+                    logger.debug("Failed to find root using second order expansion")
+            # refinement of the position by a simple center of mass of the last valid region used
+            for i0 in range(current0 - 1, current0 + 2):
+                for i1 in range(current1 - 1, current1 + 2):
+                    tmp = self.data[i0, i1]
+                    sum0 += tmp * i0
+                    sum1 += tmp * i1
+                    sum += tmp
+            if sum > 0:
+                return (sum0 / sum, sum1 / sum)
+
+        return (float(current0), float(current1))
+
+    cpdef Py_ssize_t coarse_local_maxi(self, Py_ssize_t x):
+        """Return the nearest local maximum ... without sub-pixel refinement
+
+        :param idx: start index (=row*width+column)
+        :return: local maximum index
+        """
+        return self.c_local_maxi(x)
+
+    cdef Py_ssize_t c_local_maxi(self, Py_ssize_t idx) nogil:
+        """Return the nearest local maximum without sub-pixel refinement
+
+        :param idx: start index (=row*width+column)
+        :return: local maximum index
+
+        This method is Cython only due to the NOGIL
+        """
+        cdef:
+            Py_ssize_t current0 = idx // self.width
+            Py_ssize_t current1 = idx % self.width
+            Py_ssize_t i0, i1, start0, stop0, start1, stop1, new0, new1, rng, cnt
+            mask_t m
+            data_t tmp, value, old_value
+
+        if self.has_mask and self.mask[current0, current1]:
+            #Start searching for a non masked pixel.
+            rng = 0
+            cnt = 0
+            value = self.mini
+            new0, new1 = current0, current1
+            while cnt == 0:
+                rng += 1
+                cnt = 0
+                start0 = max(0, current0 - rng)
+                stop0 = min(self.height, current0 + rng + 1)
+                start1 = max(0, current1 - rng)
+                stop1 = min(self.width, current1 + rng + 1)
+                for i0 in range(start0, stop0):
+                    for i1 in range(start1, stop1):
+                        m = not self.mask[i0, i1] 
+                        cnt += m
+                        if m:
+                            tmp = self.data[i0, i1]
+                            if tmp > value:
+                                new0, new1 = i0, i1
+                                value = tmp
+            current0, current1 = new0, new1
+        else:
+            value = self.data[current0, current1]
+        
+        old_value = value -1
+        new0, new1 = current0, current1
+
+        while value > old_value:
+            old_value = value
+            start0 = max(0, current0 - 1)
+            stop0 = min(self.height, current0 + 2)
+            start1 = max(0, current1 - 1)
+            stop1 = min(self.width, current1 + 2)
+            for i0 in range(start0, stop0):
+                for i1 in range(start1, stop1):
+                    if self.has_mask and self.mask[current0, current1]:
+                        continue
+                    tmp = self.data[i0, i1]
+                    if tmp > value:
+                        new0, new1 = i0, i1
+                        value = tmp
+            current0, current1 = new0, new1
+        return self.width * current0 + current1
+
+    def map_coordinates(self, coordinates):
+        """Map coordinates of the array on the image
+
+        :param coordinates: 2-tuple of array of the same size (row_array, column_array)
+        :return: array of values at given coordinates
+        """
+        cdef:
+            data_t[:] d0, d1, res
+            Py_ssize_t size, i
+        shape = coordinates[0].shape
+        size = coordinates[0].size
+        d0 = numpy.ascontiguousarray(coordinates[0].ravel(), dtype=data_d)
+        d1 = numpy.ascontiguousarray(coordinates[1].ravel(), dtype=data_d)
+        assert size == d1.size
+        res = numpy.empty(size, dtype=data_d)
+        with nogil:
+            for i in range(size):
+                res[i] = self.c_funct(d1[i], d0[i])
+        return numpy.asarray(res).reshape(shape)
+
+    def profile_line(self, src, dst, int linewidth=1, method='mean'):
+        """Return the mean or sum of intensity profile of an image measured
+        along a scan line.
+
+        :param src: The start point of the scan line.
+        :type src: 2-tuple of numeric scalar
+        :param dst: The end point of the scan line.
+            The destination point is included in the profile,
+            in contrast to standard numpy indexing.
+        :type dst: 2-tuple of numeric scalar
+        :param int linewidth: Width of the scanline (unit image pixel).
+        :param str method: 'mean' or 'sum' depending if we want to compute the
+            mean intensity along the line or the sum.
+        :return: The intensity profile along the scan line.
+            The length of the profile is the ceil of the computed length
+            of the scan line.
+        :rtype: 1d array
+
+        Inspired from skimage
+        """
+        cdef:
+            data_t src_row, src_col, dst_row, dst_col, d_row, d_col
+            data_t length, col_width, row_width, sum, row, col, new_row, new_col, val
+            Py_ssize_t lengt, i, j, cnt
+            bint compute_mean
+            data_t[::1] result
+        src_row, src_col = src
+        dst_row, dst_col = dst
+        if (src_row == dst_row) and (src_col == dst_col):
+            logger.warning("Source and destination points are the same")
+            return numpy.array([self.c_funct(src_col, src_row)])
+        d_row = dst_row - src_row
+        d_col = dst_col - src_col
+
+        # Offsets to deal with linewidth
+        length = sqrt(d_row * d_row + d_col * d_col)
+        row_width = d_col / length
+        col_width = - d_row / length
+
+        lengt = <int> ceil(length + 1)
+        d_row /= <data_t> (lengt -1)
+        d_col /= <data_t> (lengt -1)
+
+        result = numpy.zeros(lengt, dtype=data_d)
+
+        # Offset position to the center of the bottom pixels of the profile
+        src_row -= row_width * (linewidth - 1) / 2.
+        src_col -= col_width * (linewidth - 1) / 2.
+
+        compute_mean = (method == 'mean')
+        with nogil:
+            for i in range(lengt):
+                sum = 0
+                cnt = 0
+
+                row = src_row + i * d_row
+                col = src_col + i * d_col
+
+                for j in range(linewidth):
+                    new_row = row + j * row_width
+                    new_col = col + j * col_width
+                    if ((new_col >= 0) and (new_col < self.width) and
+                            (new_row >= 0) and (new_row < self.height)):
+                        val = self.c_funct(new_col, new_row)
+                        if isfinite(val):
+                            cnt += 1
+                            sum += val  
+                if cnt:
+                        if compute_mean:
+                            result[i] += sum / cnt
+                        else:
+                            result[i] += sum
+                elif compute_mean:
+                    result[i] += NAN
+        # Ensures the result is exported as numpy array and not memory view.
+        return numpy.asarray(result)
diff --git a/src/silx/image/marchingsquares/__init__.py b/src/silx/image/marchingsquares/__init__.py
new file mode 100644
index 0000000..a47a7f6
--- /dev/null
+++ b/src/silx/image/marchingsquares/__init__.py
@@ -0,0 +1,117 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2018 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+"""
+This module provides implementations based on marching squares algorithms.
+
+The main implementation is done by :class:`MarchingSquaresMergeImpl`. It was
+designed to speed up the computation of iso surface using Cython and OpenMP.
+It also provides features like support of mask, and cache of min/max per tiles
+which is very efficient to find many iso contours from image gradient.
+
+Utilitary functions are provided as facade for simple use.
+:meth:`find_contours` to find iso contours from an image and using the same
+main signature as `find_contours` from `skimage`, but supporting mask.
+And :meth:`find_pixels` which returns a set of pixel coords containing the
+points of the iso contours.
+"""
+
+__authors__ = ["V. Valls"]
+__license__ = "MIT"
+__date__ = "02/07/2018"
+
+
+from ._mergeimpl import MarchingSquaresMergeImpl
+
+
+def _factory(engine, image, mask):
+    """Factory to create the marching square implementation from the engine
+    name"""
+    if engine == "merge":
+        return MarchingSquaresMergeImpl(image, mask)
+    elif engine == "skimage":
+        from _skimage import MarchingSquaresSciKitImage
+        return MarchingSquaresSciKitImage(image, mask)
+    else:
+        raise ValueError("Engine '%s' is not supported ('merge' or 'skimage' expected).")
+
+
+def find_pixels(image, level, mask=None):
+    """
+    Find the pixels following the iso contours at the given `level`.
+
+    These pixels are localized by the bound of the segment generated by the
+    iso contour algorithm.
+
+    The result is returned as a numpy array storing a list of coordinates y/x.
+
+    .. code-block:: python
+
+        # Example using a mask
+        shape = 100, 100
+        image = numpy.random.random(shape)
+        mask = numpy.random.random(shape) < 0.01
+        pixels = silx.image.marchingsquares.find_pixels(image, 0.5, mask=mask)
+
+    :param numpy.ndarray image: Image to process
+    :param float level: Level of the requested iso contours.
+    :param numpy.ndarray mask: An optional mask (a non-zero value invalidate
+        the pixels of the image)
+    :returns: An array of coordinates in y/x
+    :rtype: numpy.ndarray
+    """
+    assert(image is not None)
+    if mask is not None:
+        assert(image.shape == mask.shape)
+    engine = "merge"
+    impl = _factory(engine, image, mask)
+    return impl.find_pixels(level)
+
+
+def find_contours(image, level, mask=None):
+    """
+    Find the iso contours at the given `level`.
+
+    The result is returned as a list of polygons.
+
+    .. code-block:: python
+
+        # Example using a mask
+        shape = 100, 100
+        image = numpy.random.random(shape)
+        mask = numpy.random.random(shape) < 0.01
+        polygons = silx.image.marchingsquares.find_contours(image, 0.5, mask=mask)
+
+    :param numpy.ndarray image: Image to process
+    :param float level: Level of the requested iso contours.
+    :param numpy.ndarray mask: An optional mask (a non-zero value invalidate
+        the pixels of the image)
+    :returns: A list of array containing y-x coordinates of points
+    :rtype: List[numpy.ndarray]
+    """
+    assert(image is not None)
+    if mask is not None:
+        assert(image.shape == mask.shape)
+    engine = "merge"
+    impl = _factory(engine, image, mask)
+    return impl.find_contours(level)
diff --git a/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
diff --git a/src/silx/image/marchingsquares/_skimage.py b/src/silx/image/marchingsquares/_skimage.py
new file mode 100644
index 0000000..d49eeb0
--- /dev/null
+++ b/src/silx/image/marchingsquares/_skimage.py
@@ -0,0 +1,139 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2018 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+
+__authors__ = ["V. Valls"]
+__license__ = "MIT"
+__date__ = "05/04/2018"
+
+
+import numpy
+import skimage.measure
+
+
+class MarchingSquaresSciKitImage(object):
+    """Reference implementation of a marching squares using sci-kit image.
+
+    It uses `skimage.measure.find_contours` to find iso contours taking care of
+    an optional mask. As result the computation is not accurate but can be used
+    as reference for benchmark or for testing the API without compiling the
+    cython part of silx.
+
+    :param numpy.ndarray image: 2d-image containing the values
+    :param numpy.ndarray mask: Optional 2d-image containing mask to cancel
+        mask on part of the image. A `0` means the pixel at this location is
+        valid, else the pixel from the image will not be used.
+    """
+
+    def __init__(self, image, mask=None):
+        self._image = image
+        self._mask = mask
+
+    _deltas = [(0.0, 0.0), (0.99, 0.0), (0.0, 0.99), (0.99, 0.99)]
+
+    def _flag_coord_over_mask(self, coord):
+        """Flag coord over the mask as NaN"""
+        for dx, dy in self._deltas:
+            if self._mask[int(coord[0] + dx), int(coord[1] + dy)] != 0:
+                return float("nan"), float("nan")
+        return coord
+
+    def find_pixels(self, level):
+        """
+        Compute the pixels from the image over the requested iso contours
+        at this `level`.
+
+        This implementation have to use `skimage.measure.find_contours` then
+        it is not accurate nor efficient.
+
+        :param float level: Level of the requested iso contours.
+        :returns: An array of y-x coordinates.
+        :rtype: numpy.ndarray
+        """
+        polylines = skimage.measure.find_contours(self._image, level=level)
+        size = 0
+        for polyline in polylines:
+            size += len(polyline)
+        result = numpy.empty((size, 2), dtype=numpy.int32)
+        size = 0
+        delta = numpy.array([0.5, 0.5])
+        for polyline in polylines:
+            if len(polyline) == 0:
+                continue
+            integer_polyline = numpy.floor(polyline + delta)
+            result[size:size + len(polyline)] = integer_polyline
+            size += len(polyline)
+
+        if len(result) == 0:
+            return result
+
+        if self._mask is not None:
+            # filter out pixels over the mask
+            x_dim = self._image.shape[1]
+            indexes = result[:, 0] * x_dim + result[:, 1]
+            indexes = indexes.ravel()
+            mask = self._mask.ravel()
+            indexes = numpy.unique(indexes)
+            indexes = indexes[mask[indexes] == 0]
+            pixels = numpy.concatenate((indexes // x_dim, indexes % x_dim))
+            pixels.shape = 2, -1
+            pixels = pixels.T
+            result = pixels
+        else:
+            # Note: Cound be done using a single line numpy.unique(result, axis=0)
+            # But here it supports Debian 8
+            x_dim = self._image.shape[1]
+            indexes = result[:, 0] * x_dim + result[:, 1]
+            indexes = indexes.ravel()
+            indexes = numpy.unique(indexes)
+            pixels = numpy.concatenate((indexes // x_dim, indexes % x_dim))
+            pixels.shape = 2, -1
+            pixels = pixels.T
+            result = pixels
+        return result
+
+    def find_contours(self, level):
+        """
+        Compute the list of polygons of the iso contours at this `level`.
+
+        If no mask is involved, the result is the same as
+        `skimage.measure.find_contours`.
+
+        If the result have to be filtered with a mask, the result is not
+        accurate nor efficient. Polygons are not splited, but only points are
+        filtered out using NaN coordinates. This could create artefacts.
+
+        :param float level: Level of the requested iso contours.
+        :returns: A list of array containg y-x coordinates of points
+        :rtype: List[numpy.ndarray]
+        """
+        polylines = skimage.measure.find_contours(self._image, level=level)
+        if self._mask is None:
+            return polylines
+        result = []
+        for polyline in polylines:
+            polyline = map(self._flag_coord_over_mask, polyline)
+            polyline = list(polyline)
+            polyline = numpy.array(polyline)
+            result.append(polyline)
+        return result
diff --git a/src/silx/image/marchingsquares/include/patterns.h b/src/silx/image/marchingsquares/include/patterns.h
new file mode 100644
index 0000000..ff86cc1
--- /dev/null
+++ b/src/silx/image/marchingsquares/include/patterns.h
@@ -0,0 +1,89 @@
+# /*##########################################################################
+#
+# Copyright (c) 2018 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+
+/** This header provides static lookup table used by marching square
+ * algorithms.
+ */
+
+#ifndef __MARCHINGSQUARE_PATTERNS_H__
+#define __MARCHINGSQUARE_PATTERNS_H__
+
+/**
+ * Array containing pixel's coordinate linked together by an edge index.
+ */
+const unsigned char EDGE_TO_POINT[][2] = {
+	{0, 0},
+	{1, 0},
+	{1, 1},
+	{0, 1},
+	{0, 0}
+};
+
+/**
+ * Array of index containing 5 values: {nb, seg1b, seg2e, seg2b, seg2e}
+ * nb: number of segments (up to 2)
+ * seg1b: index of the start of the 1st edge
+ * seg1e: index of the end of the 1st edge
+ * seg2b: index of the start of the 2nd edge
+ * seg2e: index of the end of the 2nd edge
+ */
+const unsigned char CELL_TO_EDGE[][5] = {
+	{0, 0, 0, 0, 0},  // Case 0: 0000: nothing
+	{1, 0, 3, 0, 0},  // Case 1: 0001
+	{1, 0, 1, 0, 0},  // Case 2: 0010
+	{1, 1, 3, 0, 0},  // Case 3: 0011
+
+	{1, 1, 2, 0, 0},  // Case 4: 0100
+	{2, 0, 1, 2, 3},  // Case 5: 0101 > ambiguous
+	{1, 0, 2, 0, 0},  // Case 6: 0110
+	{1, 2, 3, 0, 0},  // Case 7: 0111
+
+	{1, 2, 3, 0, 0},  // Case 8: 1000
+	{1, 0, 2, 0, 0},  // Case 9: 1001
+	{2, 0, 3, 1, 2},  // Case 10: 1010 > ambiguous
+	{1, 1, 2, 0, 0},  // Case 11: 1011
+
+	{1, 1, 3, 0, 0},  // Case 12: 1100
+	{1, 0, 1, 0, 0},  // Case 13: 1101
+	{1, 0, 3, 0, 0},  // Case 14: 1110
+	{0, 0, 0, 0, 0},  // Case 15: 1111
+};
+
+
+typedef struct coord_t {
+	short x;
+	short y;
+
+	bool operator<(const coord_t& other) const {
+	    if (y < other.y) {
+	        return true;
+	    } else if (y == other.y) {
+	        return x < other.x;
+	    } else {
+	        return false;
+	    }
+	}
+} coord_t;
+
+#endif /*__MARCHINGSQUARE_PATTERNS_H__*/
diff --git a/src/silx/image/marchingsquares/setup.py b/src/silx/image/marchingsquares/setup.py
new file mode 100644
index 0000000..95998ab
--- /dev/null
+++ b/src/silx/image/marchingsquares/setup.py
@@ -0,0 +1,51 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2016-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+
+__authors__ = ["V. Valls"]
+__license__ = "MIT"
+__date__ = "23/04/2018"
+
+import os
+import numpy
+from numpy.distutils.misc_util import Configuration
+
+
+def configuration(parent_package='', top_path=None):
+    config = Configuration('marchingsquares', parent_package, top_path)
+    config.add_subpackage('test')
+
+    silx_include = os.path.join(top_path, "src", "silx", "utils", "include")
+    config.add_extension('_mergeimpl',
+                         sources=['_mergeimpl.pyx'],
+                         include_dirs=[numpy.get_include(), silx_include],
+                         language='c++',
+                         extra_link_args=['-fopenmp'],
+                         extra_compile_args=['-fopenmp'])
+
+    return config
+
+
+if __name__ == "__main__":
+    from numpy.distutils.core import setup
+    setup(configuration=configuration)
diff --git a/src/silx/image/marchingsquares/test/__init__.py b/src/silx/image/marchingsquares/test/__init__.py
new file mode 100644
index 0000000..776bb73
--- /dev/null
+++ b/src/silx/image/marchingsquares/test/__init__.py
@@ -0,0 +1,24 @@
+# -*- coding: utf-8 -*-
+#
+#    Project: silx
+#             https://github.com/silx-kit/silx
+#
+#    Copyright (C) 2012-2016  European Synchrotron Radiation Facility, Grenoble, France
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
diff --git a/src/silx/image/marchingsquares/test/test_funcapi.py b/src/silx/image/marchingsquares/test/test_funcapi.py
new file mode 100644
index 0000000..d1be584
--- /dev/null
+++ b/src/silx/image/marchingsquares/test/test_funcapi.py
@@ -0,0 +1,92 @@
+# -*- coding: utf-8 -*-
+#
+#    Project: silx
+#             https://github.com/silx-kit/silx
+#
+#    Copyright (C) 2012-2016  European Synchrotron Radiation Facility, Grenoble, France
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+
+__authors__ = ["V. Valls"]
+__license__ = "MIT"
+__date__ = "17/04/2018"
+
+import unittest
+import numpy
+import silx.image.marchingsquares
+
+
+class MockMarchingSquares(object):
+
+    last = None
+
+    def __init__(self, image, mask=None):
+        MockMarchingSquares.last = self
+        self.events = []
+        self.events.append(("image", image))
+        self.events.append(("mask", mask))
+
+    def find_pixels(self, level):
+        self.events.append(("find_pixels", level))
+        return None
+
+    def find_contours(self, level):
+        self.events.append(("find_contours", level))
+        return None
+
+
+class TestFunctionalApi(unittest.TestCase):
+    """Test that the default functional API is called using the right
+    parameters to the right location."""
+
+    def setUp(self):
+        self.old_impl = silx.image.marchingsquares.MarchingSquaresMergeImpl
+        silx.image.marchingsquares.MarchingSquaresMergeImpl = MockMarchingSquares
+
+    def tearDown(self):
+        silx.image.marchingsquares.MarchingSquaresMergeImpl = self.old_impl
+        del self.old_impl
+
+    def test_default_find_contours(self):
+        image = numpy.ones((2, 2), dtype=numpy.float32)
+        mask = numpy.zeros((2, 2), dtype=numpy.int32)
+        level = 2.5
+        silx.image.marchingsquares.find_contours(image=image, level=level, mask=mask)
+        events = MockMarchingSquares.last.events
+        self.assertEqual(len(events), 3)
+        self.assertEqual(events[0][0], "image")
+        self.assertEqual(events[0][1][0, 0], 1)
+        self.assertEqual(events[1][0], "mask")
+        self.assertEqual(events[1][1][0, 0], 0)
+        self.assertEqual(events[2][0], "find_contours")
+        self.assertEqual(events[2][1], level)
+
+    def test_default_find_pixels(self):
+        image = numpy.ones((2, 2), dtype=numpy.float32)
+        mask = numpy.zeros((2, 2), dtype=numpy.int32)
+        level = 3.5
+        silx.image.marchingsquares.find_pixels(image=image, level=level, mask=mask)
+        events = MockMarchingSquares.last.events
+        self.assertEqual(len(events), 3)
+        self.assertEqual(events[0][0], "image")
+        self.assertEqual(events[0][1][0, 0], 1)
+        self.assertEqual(events[1][0], "mask")
+        self.assertEqual(events[1][1][0, 0], 0)
+        self.assertEqual(events[2][0], "find_pixels")
+        self.assertEqual(events[2][1], level)
diff --git a/src/silx/image/marchingsquares/test/test_mergeimpl.py b/src/silx/image/marchingsquares/test/test_mergeimpl.py
new file mode 100644
index 0000000..07b94b5
--- /dev/null
+++ b/src/silx/image/marchingsquares/test/test_mergeimpl.py
@@ -0,0 +1,264 @@
+# -*- coding: utf-8 -*-
+#
+#    Project: silx
+#             https://github.com/silx-kit/silx
+#
+#    Copyright (C) 2012-2016  European Synchrotron Radiation Facility, Grenoble, France
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+
+__authors__ = ["V. Valls"]
+__license__ = "MIT"
+__date__ = "18/04/2018"
+
+import unittest
+import numpy
+from .._mergeimpl import MarchingSquaresMergeImpl
+
+
+class TestMergeImplApi(unittest.TestCase):
+
+    def test_image_not_an_array(self):
+        bad_image = 1
+        self.assertRaises(ValueError, MarchingSquaresMergeImpl, bad_image)
+
+    def test_image_bad_dim(self):
+        bad_image = numpy.array([[[1.0]]])
+        self.assertRaises(ValueError, MarchingSquaresMergeImpl, bad_image)
+
+    def test_image_not_big_enough(self):
+        bad_image = numpy.array([[1.0, 1.0, 1.0, 1.0]])
+        self.assertRaises(ValueError, MarchingSquaresMergeImpl, bad_image)
+
+    def test_mask_not_an_array(self):
+        image = numpy.array([[1.0, 1.0], [1.0, 1.0]])
+        bad_mask = 1
+        self.assertRaises(ValueError, MarchingSquaresMergeImpl, image, bad_mask)
+
+    def test_mask_not_match(self):
+        image = numpy.array([[1.0, 1.0], [1.0, 1.0]])
+        bad_mask = numpy.array([[1.0, 1.0]])
+        self.assertRaises(ValueError, MarchingSquaresMergeImpl, image, bad_mask)
+
+    def test_ok_anyway_bad_type(self):
+        image = numpy.array([[1.0, 1.0], [1.0, 1.0]], dtype=numpy.int32)
+        mask = numpy.array([[1.0, 1.0], [1.0, 1.0]], dtype=numpy.float32)
+        MarchingSquaresMergeImpl(image, mask)
+
+    def test_find_contours_result(self):
+        image = numpy.zeros((2, 2))
+        image[0, 0] = 1
+        ms = MarchingSquaresMergeImpl(image)
+        polygons = ms.find_contours(0.5)
+        self.assertIsInstance(polygons, list)
+        self.assertTrue(len(polygons), 1)
+        self.assertIsInstance(polygons[0], numpy.ndarray)
+        self.assertEqual(polygons[0].shape[1], 2)
+        self.assertEqual(polygons[0].dtype.kind, "f")
+
+    def test_find_pixels_result(self):
+        image = numpy.zeros((2, 2))
+        image[0, 0] = 1
+        ms = MarchingSquaresMergeImpl(image)
+        pixels = ms.find_pixels(0.5)
+        self.assertIsInstance(pixels, numpy.ndarray)
+        self.assertEqual(pixels.shape[1], 2)
+        self.assertEqual(pixels.dtype.kind, "i")
+
+    def test_find_contours_empty_result(self):
+        image = numpy.zeros((2, 2))
+        ms = MarchingSquaresMergeImpl(image)
+        polygons = ms.find_contours(0.5)
+        self.assertIsInstance(polygons, list)
+        self.assertEqual(len(polygons), 0)
+
+    def test_find_pixels_empty_result(self):
+        image = numpy.zeros((2, 2))
+        ms = MarchingSquaresMergeImpl(image)
+        pixels = ms.find_pixels(0.5)
+        self.assertIsInstance(pixels, numpy.ndarray)
+        self.assertEqual(pixels.shape[1], 2)
+        self.assertEqual(pixels.shape[0], 0)
+        self.assertEqual(pixels.dtype.kind, "i")
+
+    def test_find_contours_yx_result(self):
+        image = numpy.zeros((2, 2))
+        image[1, 0] = 1
+        ms = MarchingSquaresMergeImpl(image)
+        polygons = ms.find_contours(0.5)
+        polygon = polygons[0]
+        self.assertTrue((polygon == (0.5, 0)).any())
+        self.assertTrue((polygon == (1, 0.5)).any())
+
+    def test_find_pixels_yx_result(self):
+        image = numpy.zeros((2, 2))
+        image[1, 0] = 1
+        ms = MarchingSquaresMergeImpl(image)
+        pixels = ms.find_pixels(0.5)
+        self.assertTrue((pixels == (1, 0)).any())
+
+
+class TestMergeImplContours(unittest.TestCase):
+
+    def test_merge_segments(self):
+        image = numpy.zeros((4, 4))
+        image[(2, 3), :] = 1
+        ms = MarchingSquaresMergeImpl(image)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 1)
+
+    def test_merge_segments_2(self):
+        image = numpy.zeros((4, 4))
+        image[(2, 3), :] = 1
+        image[2, 2] = 0
+        ms = MarchingSquaresMergeImpl(image)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 1)
+
+    def test_merge_tiles(self):
+        image = numpy.zeros((4, 4))
+        image[(2, 3), :] = 1
+        ms = MarchingSquaresMergeImpl(image, group_size=2)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 1)
+
+    def test_fully_masked(self):
+        image = numpy.zeros((5, 5))
+        image[(2, 3), :] = 1
+        mask = numpy.ones((5, 5))
+        ms = MarchingSquaresMergeImpl(image, mask)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 0)
+
+    def test_fully_masked_minmax(self):
+        """This invalidates all the tiles. The route is not the same."""
+        image = numpy.zeros((5, 5))
+        image[(2, 3), :] = 1
+        mask = numpy.ones((5, 5))
+        ms = MarchingSquaresMergeImpl(image, mask, group_size=2, use_minmax_cache=True)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 0)
+
+    def test_masked_segments(self):
+        image = numpy.zeros((5, 5))
+        image[(2, 3, 4), :] = 1
+        mask = numpy.zeros((5, 5))
+        mask[:, 2] = 1
+        ms = MarchingSquaresMergeImpl(image, mask)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 2)
+
+    def test_closed_polygon(self):
+        image = numpy.zeros((5, 5))
+        image[2, 2] = 1
+        image[1, 2] = 1
+        image[3, 2] = 1
+        image[2, 1] = 1
+        image[2, 3] = 1
+        mask = None
+        ms = MarchingSquaresMergeImpl(image, mask)
+        polygons = ms.find_contours(0.9)
+        self.assertEqual(len(polygons), 1)
+        self.assertEqual(list(polygons[0][0]), list(polygons[0][-1]))
+
+    def test_closed_polygon_between_tiles(self):
+        image = numpy.zeros((5, 5))
+        image[2, 2] = 1
+        image[1, 2] = 1
+        image[3, 2] = 1
+        image[2, 1] = 1
+        image[2, 3] = 1
+        mask = None
+        ms = MarchingSquaresMergeImpl(image, mask, group_size=2)
+        polygons = ms.find_contours(0.9)
+        self.assertEqual(len(polygons), 1)
+        self.assertEqual(list(polygons[0][0]), list(polygons[0][-1]))
+
+    def test_open_polygon(self):
+        image = numpy.zeros((5, 5))
+        image[2, 2] = 1
+        image[1, 2] = 1
+        image[3, 2] = 1
+        image[2, 1] = 1
+        image[2, 3] = 1
+        mask = numpy.zeros((5, 5))
+        mask[1, 1] = 1
+        ms = MarchingSquaresMergeImpl(image, mask)
+        polygons = ms.find_contours(0.9)
+        self.assertEqual(len(polygons), 1)
+        self.assertNotEqual(list(polygons[0][0]), list(polygons[0][-1]))
+
+    def test_ambiguous_pattern(self):
+        image = numpy.zeros((6, 8))
+        image[(3, 4), :] = 1
+        image[:, (0, -1)] = 0
+        image[3, 3] = -0.001
+        image[4, 4] = 0.0
+        mask = None
+        ms = MarchingSquaresMergeImpl(image, mask)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 2)
+
+    def test_ambiguous_pattern_2(self):
+        image = numpy.zeros((6, 8))
+        image[(3, 4), :] = 1
+        image[:, (0, -1)] = 0
+        image[3, 3] = +0.001
+        image[4, 4] = 0.0
+        mask = None
+        ms = MarchingSquaresMergeImpl(image, mask)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 1)
+
+    def count_closed_polygons(self, polygons):
+        closed = 0
+        for polygon in polygons:
+            if list(polygon[0]) == list(polygon[-1]):
+                closed += 1
+        return closed
+
+    def test_image(self):
+        # example from skimage
+        x, y = numpy.ogrid[-numpy.pi:numpy.pi:100j, -numpy.pi:numpy.pi:100j]
+        image = numpy.sin(numpy.exp((numpy.sin(x)**3 + numpy.cos(y)**2)))
+        mask = None
+        ms = MarchingSquaresMergeImpl(image, mask)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 11)
+        self.assertEqual(self.count_closed_polygons(polygons), 3)
+
+    def test_image_tiled(self):
+        # example from skimage
+        x, y = numpy.ogrid[-numpy.pi:numpy.pi:100j, -numpy.pi:numpy.pi:100j]
+        image = numpy.sin(numpy.exp((numpy.sin(x)**3 + numpy.cos(y)**2)))
+        mask = None
+        ms = MarchingSquaresMergeImpl(image, mask, group_size=50)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 11)
+        self.assertEqual(self.count_closed_polygons(polygons), 3)
+
+    def test_image_tiled_minmax(self):
+        # example from skimage
+        x, y = numpy.ogrid[-numpy.pi:numpy.pi:100j, -numpy.pi:numpy.pi:100j]
+        image = numpy.sin(numpy.exp((numpy.sin(x)**3 + numpy.cos(y)**2)))
+        mask = None
+        ms = MarchingSquaresMergeImpl(image, mask, group_size=50, use_minmax_cache=True)
+        polygons = ms.find_contours(0.5)
+        self.assertEqual(len(polygons), 11)
+        self.assertEqual(self.count_closed_polygons(polygons), 3)
diff --git a/src/silx/image/medianfilter.py b/src/silx/image/medianfilter.py
new file mode 100644
index 0000000..857f73d
--- /dev/null
+++ b/src/silx/image/medianfilter.py
@@ -0,0 +1,114 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2017-2018 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+
+"""
+This module provides :func:`medfilt2d`, a 2D median filter function
+with the choice between 2 implementations: 'cpp' and 'opencl'.
+"""
+
+__authors__ = ["H. Payno"]
+__license__ = "MIT"
+__date__ = "04/05/2017"
+
+
+import logging
+
+from silx.math import medianfilter as medianfilter_cpp
+from silx.opencl import ocl as _ocl
+if _ocl is not None:
+    from silx.opencl import medfilt as medfilt_opencl
+else:  # No OpenCL device or pyopencl not installed
+    medfilt_opencl = None
+
+
+_logger = logging.getLogger(__name__)
+
+
+MEDFILT_ENGINES = ['cpp', 'opencl']
+
+
+def medfilt2d(image, kernel_size=3, engine='cpp'):
+    """Apply a median filter on an image.
+
+    This median filter is using a 'nearest' padding for values
+    past the array edges. If you want more padding options or
+    functionalities for the median filter (conditional filter 
+    for example) please have a look at
+    :mod:`silx.math.medianfilter`.
+
+    :param numpy.ndarray image: the 2D array for which we want to apply
+        the median filter.
+    :param kernel_size: the dimension of the kernel.
+        Kernel size must be odd.
+        If a scalar is given, then it is used as the size in both dimension.
+        Default: (3, 3)
+    :type kernel_size: A int or a list of 2 int (kernel_height, kernel_width)
+    :param engine: the type of implementation to use.
+        Valid values are: 'cpp' (default) and 'opencl'
+
+    :returns: the array with the median value for each pixel.
+
+    .. note::  if the opencl implementation is requested but
+        is not present or fails, the cpp implementation is called.
+
+    """
+    if engine not in MEDFILT_ENGINES:
+        err = 'silx doesn\'t have an implementation for the requested engine: '
+        err += '%s' % engine
+        raise ValueError(err)
+
+    if len(image.shape) != 2:
+        raise ValueError('medfilt2d deals with arrays of dimension 2 only')
+
+    if engine == 'cpp':
+        return medianfilter_cpp.medfilt(data=image,
+                                        kernel_size=kernel_size,
+                                        conditional=False)
+    elif engine == 'opencl':
+        if medfilt_opencl is None:
+            wrn = 'opencl median filter not available. '
+            wrn += 'Launching cpp implementation.'
+            _logger.warning(wrn)
+            # instead call the cpp implementation
+            return medianfilter_cpp.medfilt(data=image,
+                                            kernel_size=kernel_size,
+                                            conditional=False)
+        else:
+            try:
+                medianfilter = medfilt_opencl.MedianFilter2D(image.shape,
+                                                             devicetype="gpu")
+                res = medianfilter.medfilt2d(image, kernel_size)
+            except(RuntimeError, MemoryError, ImportError):
+                wrn = 'Exception occured in opencl median filter. '
+                wrn += 'To get more information see debug log.'
+                wrn += 'Launching cpp implementation.'
+                _logger.warning(wrn)
+                _logger.debug("median filter - openCL implementation issue.",
+                              exc_info=True)
+                # instead call the cpp implementation
+                res = medianfilter_cpp.medfilt(data=image,
+                                               kernel_size=kernel_size,
+                                               conditional=False)
+
+        return res
diff --git a/src/silx/image/phantomgenerator.py b/src/silx/image/phantomgenerator.py
new file mode 100644
index 0000000..10b249b
--- /dev/null
+++ b/src/silx/image/phantomgenerator.py
@@ -0,0 +1,160 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2016 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+
+__authors__ = ["N. Vigano", "H. Payno", "P. Paleo"]
+__license__ = "MIT"
+__date__ = "19/09/2017"
+
+import numpy
+
+
+class PhantomGenerator(object):
+    """
+    Class for generating Phantoms
+    """
+
+    class _Ellipsoid:
+        def __init__(self, a, b, c, x0, y0, z0, alpha, mu):
+            self.a = a
+            self.b = b
+            self.c = c
+            self.x0 = x0
+            self.y0 = y0
+            self.z0 = z0
+            self.alpha = alpha * numpy.pi / 180.0
+            self.mu = mu
+            self.cosAlpha = numpy.cos(self.alpha)
+            self.sinAlpha = numpy.sin(self.alpha)
+
+    SHEPP_LOGAN = [
+        #         a       b      c       x0     y0      z0     alpha  mu
+        _Ellipsoid(0.69, 0.92, 0.90, 0.0, 0.0, 0.0, 0.0, 0.10),
+        _Ellipsoid(0.6624, 0.874, 0.88, 0.0, -0.02, 0.0, 0.0, -0.08),
+        _Ellipsoid(0.11, 0.31, 0.21, 0.22, -0.0, 0.0, -18.0, -0.02),
+        _Ellipsoid(0.16, 0.41, 0.22, -0.22, 0.0, -0.25, 18.0, -0.02),
+        _Ellipsoid(0.21, 0.25, 0.35, 0.0, 0.35, -0.25, 0.0, 0.03),
+        _Ellipsoid(0.046, 0.046, 0.046, 0.0, 0.10, -0.25, 0.0, 0.01),
+        _Ellipsoid(0.046, 0.046, 0.02, 0.0, -0.10, -0.25, 0.0, 0.01),
+        _Ellipsoid(0.046, 0.023, 0.02, -0.08, -0.605, -0.25, 0.0, 0.01),
+        _Ellipsoid(0.023, 0.023, 0.10, 0.0, -0.605, -0.25, 0.0, 0.01),
+        _Ellipsoid(0.023, 0.046, 0.10, 0.06, -0.605, -0.25, 0.0, 0.01)
+    ]
+
+    @staticmethod
+    def get2DPhantomSheppLogan(n, ellipsoidID=None):
+        """
+        generate a classical 2D shepp logan phantom.
+
+        :param n: The width (and height) of the phantom to generate
+        :param ellipsoidID: The Id of the ellipsoid to pick. If None will
+                            produce every ellipsoid
+        :return numpy.ndarray: shepp logan phantom
+        """
+        assert(ellipsoidID is None or (ellipsoidID >= 0 and ellipsoidID < len(PhantomGenerator.SHEPP_LOGAN)))
+        if ellipsoidID is None:
+            area = PhantomGenerator._get2DPhantom(n,
+                                                  PhantomGenerator.SHEPP_LOGAN)
+        else:
+            area = PhantomGenerator._get2DPhantom(n,
+                                                  [PhantomGenerator.SHEPP_LOGAN[ellipsoidID]])
+
+        indices = numpy.abs(area) > 0
+        area[indices] = numpy.multiply(area[indices] + 0.1, 5)
+        return area / 100.0
+
+    @staticmethod
+    def _get2DPhantom(n, phantomSpec):
+        area = numpy.ndarray(shape=(n, n))
+        area.fill(0.)
+
+        count = 0
+        for ell in phantomSpec:
+            count = count+1
+            for x in range(n):
+                sumSquareXandY = PhantomGenerator._getSquareXandYsum(n, x, ell)
+                indices = sumSquareXandY <= 1
+                area[indices, x] = ell.mu
+        return area
+
+    @staticmethod
+    def _getSquareXandYsum(n, x, ell):
+        supportX1 = numpy.ndarray(shape=(n, ))
+        supportX2 = numpy.ndarray(shape=(n, ))
+        support_consts = numpy.ndarray(shape=(n, ))
+
+        xScaled = float(2*x-n)/float(n)
+        xCos = xScaled * ell.cosAlpha
+        xSin = -xScaled * ell.sinAlpha
+        supportX1.fill(xCos)
+        supportX2.fill(xSin)
+
+        supportY1 = numpy.arange(n)
+        support_consts.fill(2.)
+        supportY1 = numpy.multiply(support_consts, supportY1)
+        support_consts.fill(n)
+        supportY1 = numpy.subtract(supportY1, support_consts)
+        support_consts.fill(n)
+        supportY1 = numpy.divide(supportY1, support_consts)
+        supportY2 = numpy.array(supportY1)
+
+        support_consts.fill(ell.sinAlpha)
+        supportY1 = numpy.add(supportX1,
+                              numpy.multiply(supportY1, support_consts))
+        support_consts.fill(ell.cosAlpha)
+        supportY2 = numpy.add(supportX2,
+                              numpy.multiply(supportY2, support_consts))
+
+        support_consts.fill(ell.x0)
+        supportY1 = numpy.subtract(supportY1, support_consts)
+        support_consts.fill(ell.y0)
+        supportY2 = numpy.subtract(supportY2, support_consts)
+
+        support_consts.fill(ell.a)
+        supportY1 = numpy.power((numpy.divide(supportY1, support_consts)),
+                                2)
+        support_consts.fill(ell.b)
+        supportY2 = numpy.power(numpy.divide(supportY2, support_consts),
+                                2)
+
+        return numpy.add(supportY1, supportY2)
+
+    @staticmethod
+    def _getSquareZ(n, ell):
+        supportZ1 = numpy.arange(n)
+        support_consts = numpy.ndarray(shape=(n, ))
+        support_consts.fill(2.)
+        supportZ1 = numpy.multiply(support_consts, supportZ1)
+        support_consts.fill(n)
+        supportZ1 = numpy.subtract(supportZ1, support_consts)
+        support_consts.fill(n)
+        supportZ1 = numpy.divide(supportZ1, support_consts)
+
+        support_consts.fill(ell.z0)
+        supportZ1 = numpy.subtract(supportZ1, ell.z0)
+
+        support_consts.fill(ell.c)
+        return numpy.power(numpy.divide(supportZ1, support_consts),
+                           2)
+
diff --git a/src/silx/image/projection.py b/src/silx/image/projection.py
new file mode 100644
index 0000000..5c76c35
--- /dev/null
+++ b/src/silx/image/projection.py
@@ -0,0 +1,25 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+
+from silx.opencl.projection import *
diff --git a/src/silx/image/reconstruction.py b/src/silx/image/reconstruction.py
new file mode 100644
index 0000000..875b66b
--- /dev/null
+++ b/src/silx/image/reconstruction.py
@@ -0,0 +1,25 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+
+from silx.opencl.reconstruction import *
diff --git a/src/silx/image/setup.py b/src/silx/image/setup.py
new file mode 100644
index 0000000..69d5b1b
--- /dev/null
+++ b/src/silx/image/setup.py
@@ -0,0 +1,47 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2016-2018 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+
+__authors__ = ["J. Kieffer"]
+__license__ = "MIT"
+__date__ = "05/04/2018"
+
+from numpy.distutils.misc_util import Configuration
+
+
+def configuration(parent_package='', top_path=None):
+    config = Configuration('image', parent_package, top_path)
+    config.add_subpackage('test')
+    config.add_extension('bilinear',
+                         sources=["bilinear.pyx"],
+                         language='c')
+    config.add_extension('shapes',
+                         sources=["shapes.pyx"],
+                         language='c')
+    config.add_subpackage('marchingsquares')
+    return config
+
+
+if __name__ == "__main__":
+    from numpy.distutils.core import setup
+    setup(configuration=configuration)
diff --git a/src/silx/image/shapes.pyx b/src/silx/image/shapes.pyx
new file mode 100644
index 0000000..9284811
--- /dev/null
+++ b/src/silx/image/shapes.pyx
@@ -0,0 +1,321 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+"""This module provides functions making masks on an image.
+
+- :func:`circle_fill` function generates coordinates of a circle in an image.
+- :func:`draw_line` function generates coordinates of a line in an image.
+- :func:`polygon_fill_mask` function generates a mask from a set of points
+  defining a polygon.
+
+The :class:`Polygon` class provides checking if a point is inside a polygon.
+
+The whole module uses the (row, col) (i.e., (y, x))) convention
+for 2D coordinates.
+"""
+
+__authors__ = ["Jérôme Kieffer", "T. Vincent"]
+__license__ = "MIT"
+__date__ = "15/02/2019"
+__status__ = "dev"
+
+
+cimport cython
+import numpy
+from libc.math cimport ceil, fabs
+
+
+cdef class Polygon(object):
+    """Define a polygon that provides inside check and mask generation.
+
+    :param vertices: corners of the polygon
+    :type vertices: Nx2 array of floats of (row, col)
+    """
+
+    cdef float[:,:] vertices
+    cdef int nvert
+
+    def __init__(self, vertices):
+        self.vertices = numpy.ascontiguousarray(vertices, dtype=numpy.float32)
+        if self.vertices.ndim != 2 or self.vertices.shape[1] != 2:
+             raise ValueError("A list of 2d vertices is expected (n,2 dimensional array)")
+        self.nvert = self.vertices.shape[0]
+
+    def is_inside(self, row, col):
+        """Check if (row, col) is inside or outside the polygon
+
+        :param float row:
+        :param float col:
+        :return: True if position is inside polygon, False otherwise
+        """
+        return self.c_is_inside(row, col)
+
+    @cython.cdivision(True)
+    @cython.wraparound(False)
+    @cython.boundscheck(False)
+    cdef bint c_is_inside(self, float row, float col) nogil:
+        """Check if (row, col) is inside or outside the polygon
+
+        Pure C_Cython class implementation.
+
+        :param float row:
+        :param float col:
+        :return: True if position is inside polygon, False otherwise
+        """
+        cdef int index, is_inside
+        cdef float pt1x, pt1y, pt2x, pt2y, xinters
+        is_inside = 0
+
+        pt1x = self.vertices[self.nvert-1, 1]
+        pt1y = self.vertices[self.nvert-1, 0]
+        for index in range(self.nvert):
+            pt2x = self.vertices[index, 1]
+            pt2y = self.vertices[index, 0]
+
+            if (((pt1y <= row and row < pt2y) or
+                    (pt2y <= row and row < pt1y)) and
+                    # Extra (optional) condition to avoid some computation
+                    (col <= pt1x or col <= pt2x)):
+                xinters = (row - pt1y) * (pt2x - pt1x) / (pt2y - pt1y) + pt1x
+                is_inside ^= col < xinters
+            pt1x, pt1y = pt2x, pt2y
+        return is_inside
+
+    @cython.cdivision(True)
+    @cython.wraparound(False)
+    @cython.boundscheck(False)
+    def make_mask(self, int height, int width):
+        """Create a mask array representing the filled polygon
+
+        :param int height: Height of the mask array
+        :param int width: Width of the mask array
+        :return: 2D array (height, width)
+        """
+        cdef unsigned char[:, :] mask = numpy.zeros((height, width),
+                                                    dtype=numpy.uint8)
+        cdef int row_min, row_max, col_min, col_max  # mask subpart to update
+        cdef int row, col, index  # Loop indixes
+        cdef float pt1x, pt1y, pt2x, pt2y  # segment end points
+        cdef int xinters, is_inside, current
+
+        row_min = max(int(min(self.vertices[:, 0])), 0)
+        row_max = min(int(max(self.vertices[:, 0])) + 1, height)
+
+        # Can be replaced by prange(row_min, row_max, nogil=True)
+        with nogil:
+            for row in range(row_min, row_max):
+                # For each line of the image, mark intersection of all segments
+                # in the line and then run a xor scan to fill inner parts
+                # Adapted from http://alienryderflex.com/polygon_fill/
+                pt1x = self.vertices[self.nvert-1, 1]
+                pt1y = self.vertices[self.nvert-1, 0]
+                col_min = width - 1
+                col_max = 0
+                is_inside = 0  # Init with whether first col is inside or not
+
+                for index in range(self.nvert):
+                    pt2x = self.vertices[index, 1]
+                    pt2y = self.vertices[index, 0]
+
+                    if ((pt1y <= row and row < pt2y) or
+                            (pt2y <= row and row < pt1y)):
+                        # Intersection casted to int so that ]x, x+1] => x
+                        xinters = (<int>ceil(pt1x + (row - pt1y) *
+                                   (pt2x - pt1x) / (pt2y - pt1y))) - 1
+
+                        # Update column range to patch
+                        if xinters < col_min:
+                            col_min = xinters
+                        if xinters > col_max:
+                            col_max = xinters
+
+                        if xinters < 0:
+                            # Add an intersection to init value of xor scan
+                            is_inside ^= 1
+                        elif xinters < width:
+                            # Mark intersection in mask
+                            mask[row, xinters] ^= 1
+                        # else: do not consider intersection on the right
+
+                    pt1x, pt1y = pt2x, pt2y
+
+                if col_min < col_max:
+                    # Clip column range to mask
+                    if col_min < 0:
+                        col_min = 0
+                    if col_max > width - 1:
+                        col_max = width - 1
+
+                    # xor exclusive scan
+                    for col in range(col_min, col_max + 1):
+                        current = mask[row, col]
+                        mask[row, col] = is_inside
+                        is_inside = current ^ is_inside
+
+        # Ensures the result is exported as numpy array and not memory view.
+        return numpy.asarray(mask)
+
+
+def polygon_fill_mask(vertices, shape):
+    """Return a mask of boolean, True for pixels inside a polygon.
+
+    :param vertices: Strip of segments end points (row, column) or (y, x)
+    :type vertices: numpy.ndarray like container of dimension Nx2
+    :param shape: size of the mask as (height, width)
+    :type shape: 2-tuple of int
+    :return: Mask corresponding to the polygon
+    :rtype: numpy.ndarray of dimension shape
+    """
+    return Polygon(vertices).make_mask(shape[0], shape[1])
+
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+def draw_line(int row0, int col0, int row1, int col1, int width=1):
+    """Line includes both end points.
+    Width is handled by drawing parallel lines, so junctions of lines belonging
+    to different octant with width > 1 will not look nice.
+
+    Using Bresenham line algorithm:
+    Bresenham, J. E.
+    Algorithm for computer control of a digital plotter.
+    IBM Systems Journal. Vol 4 No 1. 1965. pp 25-30
+
+    :param int row0: Start point row
+    :param int col0: Start point col
+    :param int row1: End point row
+    :param int col1: End point col
+    :param int width: Thickness of the line in pixels (default 1)
+                      Width must be at least 1.
+    :return: Array coordinates of points inside the line (might be negative)
+    :rtype: 2-tuple of numpy.ndarray (rows, cols)
+    """
+    cdef int drow, dcol, invert_coords
+    cdef int db, da, delta, b, a, step_a, step_b
+    cdef int index, offset  # Loop indices
+
+    # Store coordinates of points of the line
+    cdef int[:, :] b_coords
+    cdef int[:, :] a_coords
+
+    dcol = abs(col1 - col0)
+    drow = abs(row1 - row0)
+    invert_coords = dcol < drow
+
+    if dcol == 0 and drow == 0:
+        return (numpy.array((row0,), dtype=numpy.int32),
+                numpy.array((col0,), dtype=numpy.int32))
+
+    if width < 1:
+        width = 1
+
+    # Set a and b according to segment octant
+    if not invert_coords:
+        da = dcol
+        db = drow
+        step_a = 1 if col1 > col0 else -1
+        step_b = 1 if row1 > row0 else -1
+        a = col0
+        b = row0
+
+    else:
+        da = drow
+        db = dcol
+        step_a = 1 if row1 > row0 else -1
+        step_b = 1 if col1 > col0 else -1
+        a = row0
+        b = col0
+
+    b_coords = numpy.empty((da + 1, width), dtype=numpy.int32)
+    a_coords = numpy.empty((da + 1, width), dtype=numpy.int32)
+
+    with nogil:
+        b -= (width - 1) // 2
+        delta = 2 * db - da
+        for index in range(da + 1):
+            for offset in range(width):
+                b_coords[index, offset] = b + offset
+                a_coords[index, offset] = a
+
+            if delta >= 0:  # M2: Move by step_a + step_b
+                b += step_b
+                delta -= 2 * da
+            # else M1: Move by step_a
+
+            a += step_a
+            delta += 2 * db
+
+    if not invert_coords:
+        return (numpy.asarray(b_coords).reshape(-1),
+                numpy.asarray(a_coords).reshape(-1))
+    else:
+        return (numpy.asarray(a_coords).reshape(-1),
+                numpy.asarray(b_coords).reshape(-1))
+
+
+def circle_fill(int crow, int ccol, float radius):
+    """Generates coordinate of image points lying in a disk.
+
+    :param int crow: Row of the center of the disk
+    :param int ccol: Column of the center of the disk
+    :param float radius: Radius of the disk
+    :return: Array coordinates of points inside the disk (might be negative)
+    :rtype: 2-tuple of numpy.ndarray (rows, cols)
+    """
+    cdef int i_radius, len_coords
+
+    radius = fabs(radius)
+    i_radius = <int>radius
+
+    coords = numpy.arange(-i_radius, ceil(radius) + 1,
+                          dtype=numpy.float32) ** 2
+    len_coords = len(coords)
+    # rows, cols = where(row**2 + col**2 < radius**2)
+    rows, cols = numpy.where(coords.reshape(1, len_coords) +
+                             coords.reshape(len_coords, 1) < radius ** 2)
+    return rows + crow - i_radius, cols + ccol - i_radius
+
+
+def ellipse_fill(int crow, int ccol, float radius_r, float radius_c):
+    """Generates coordinate of image points lying in a ellipse.
+
+    :param int crow: Row of the center of the ellipse
+    :param int ccol: Column of the center of the ellipse
+    :param float radius_r: Radius of the ellipse in the row
+    :param float radius_c: Radius of the ellipse in the column
+    :return: Array coordinates of points inside the ellipse (might be negative)
+    :rtype: 2-tuple of numpy.ndarray (rows, cols)
+    """
+    cdef int i_radius_r
+    cdef int i_radius_c
+
+    i_radius_r = <int>fabs(radius_r)
+    i_radius_c = <int>fabs(radius_c)
+
+    x_coords = numpy.arange(-i_radius_r, ceil(radius_r) + 1, dtype=numpy.float32).reshape(-1, 1)
+    y_coords = numpy.arange(-i_radius_c, ceil(radius_c) + 1, dtype=numpy.float32).reshape(1, -1)
+
+    # rows, cols = where(x**2 + col**2 < 1)
+    rows, cols = numpy.where(x_coords**2 / radius_r**2 + y_coords**2 / radius_c**2 < 1.0)
+    return rows + crow - i_radius_r, cols + ccol - i_radius_c
diff --git a/src/silx/image/sift.py b/src/silx/image/sift.py
new file mode 100644
index 0000000..cb1e6bd
--- /dev/null
+++ b/src/silx/image/sift.py
@@ -0,0 +1,25 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+
+from silx.opencl.sift import *
diff --git a/src/silx/image/test/__init__.py b/src/silx/image/test/__init__.py
new file mode 100644
index 0000000..40b11a1
--- /dev/null
+++ b/src/silx/image/test/__init__.py
@@ -0,0 +1,24 @@
+# -*- coding: utf-8 -*-
+#
+#    Project: silx
+#             https://github.com/silx-kit/silx
+#
+#    Copyright (C) 2012-2018  European Synchrotron Radiation Facility, Grenoble, France
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
diff --git a/src/silx/image/test/test_bb.py b/src/silx/image/test/test_bb.py
new file mode 100644
index 0000000..7427273
--- /dev/null
+++ b/src/silx/image/test/test_bb.py
@@ -0,0 +1,74 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2016 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""Basic tests for Bounding box"""
+
+__authors__ = ["H. Payno"]
+__license__ = "MIT"
+__date__ = "27/09/2019"
+
+
+import unittest
+import numpy
+from silx.image._boundingbox import _BoundingBox
+
+
+class TestBB(unittest.TestCase):
+    """Some simple test on the bounding box class"""
+    def test_creation(self):
+        """test some constructors"""
+        pts = numpy.array([(0, 0), (10, 20), (20, 0)])
+        bb = _BoundingBox.from_points(pts)
+        self.assertTrue(bb.bottom_left == (0, 0))
+        self.assertTrue(bb.top_right == (20, 20))
+        pts = numpy.array([(0, 10), (10, 20), (45, 30), (35, 0)])
+        bb = _BoundingBox.from_points(pts)
+        self.assertTrue(bb.bottom_left == (0, 0))
+        print(bb.top_right)
+        self.assertTrue(bb.top_right == (45, 30))
+
+    def test_isIn_pt(self):
+        """test the isIn function with points"""
+        bb = _BoundingBox(bottom_left=(6, 2), top_right=(12, 6))
+        self.assertTrue(bb.contains((10, 4)))
+        self.assertTrue(bb.contains((6, 2)))
+        self.assertTrue(bb.contains((12, 2)))
+        self.assertFalse(bb.contains((0, 0)))
+        self.assertFalse(bb.contains((20, 0)))
+        self.assertFalse(bb.contains((10, 0)))
+
+    def test_collide(self):
+        """test the collide function"""
+        bb1 = _BoundingBox(bottom_left=(6, 2), top_right=(12, 6))
+        self.assertTrue(bb1.collide(_BoundingBox(bottom_left=(6, 2), top_right=(12, 6))))
+        bb1 = _BoundingBox(bottom_left=(6, 2), top_right=(12, 6))
+        self.assertFalse(bb1.collide(_BoundingBox(bottom_left=(12, 2), top_right=(12, 2))))
+
+    def test_isIn_bb(self):
+        """test the isIn function with other bounding box"""
+        bb1 = _BoundingBox(bottom_left=(6, 2), top_right=(12, 6))
+        self.assertTrue(bb1.contains(_BoundingBox(bottom_left=(6, 2), top_right=(12, 6))))
+        bb1 = _BoundingBox(bottom_left=(6, 2), top_right=(12, 6))
+        self.assertTrue(bb1.contains(_BoundingBox(bottom_left=(12, 2), top_right=(12, 2))))
+        self.assertFalse(_BoundingBox(bottom_left=(12, 2), top_right=(12, 2)).contains(bb1))
diff --git a/src/silx/image/test/test_bilinear.py b/src/silx/image/test/test_bilinear.py
new file mode 100644
index 0000000..20ceb58
--- /dev/null
+++ b/src/silx/image/test/test_bilinear.py
@@ -0,0 +1,167 @@
+# -*- coding: utf-8 -*-
+#
+#    Project: silx (originally pyFAI)
+#             https://github.com/silx-kit/silx
+#
+#    Copyright (C) 2012-2017  European Synchrotron Radiation Facility, Grenoble, France
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+
+__authors__ = ["J. Kieffer"]
+__license__ = "MIT"
+__date__ = "25/11/2020"
+
+import unittest
+import numpy
+import logging
+logger = logging.getLogger(__name__)
+from ..bilinear import BilinearImage
+
+
+class TestBilinear(unittest.TestCase):
+    """basic maximum search test"""
+    N = 1000
+
+    def test_max_search_round(self):
+        """test maximum search using random points: maximum is at the pixel center"""
+        a = numpy.arange(100) - 40.
+        b = numpy.arange(100) - 60.
+        ga = numpy.exp(-a * a / 4000)
+        gb = numpy.exp(-b * b / 6000)
+        gg = numpy.outer(ga, gb)
+        b = BilinearImage(gg)
+
+        self.assertAlmostEqual(b.maxi, 1, 2, "maxi is almost 1")
+        self.assertLess(b.mini, 0.3, "mini should be around 0.23")
+
+        ok = 0
+        for s in range(self.N):
+            i, j = numpy.random.randint(100), numpy.random.randint(100)
+            k, l = b.local_maxi((i, j))
+            if abs(k - 40) > 1e-4 or abs(l - 60) > 1e-4:
+                logger.warning("Wrong guess maximum (%i,%i) -> (%.1f,%.1f)", i, j, k, l)
+            else:
+                logger.debug("Good guess maximum (%i,%i) -> (%.1f,%.1f)", i, j, k, l)
+                ok += 1
+        logger.debug("Success rate: %.1f", 100. * ok / self.N)
+        self.assertEqual(ok, self.N, "Maximum is always found")
+
+    def test_max_search_half(self):
+        """test maximum search using random points: maximum is at a pixel edge"""
+        a = numpy.arange(100) - 40.5
+        b = numpy.arange(100) - 60.5
+        ga = numpy.exp(-a * a / 4000)
+        gb = numpy.exp(-b * b / 6000)
+        gg = numpy.outer(ga, gb)
+        b = BilinearImage(gg)
+        ok = 0
+        for s in range(self.N):
+            i, j = numpy.random.randint(100), numpy.random.randint(100)
+            k, l = b.local_maxi((i, j))
+            if abs(k - 40.5) > 0.5 or abs(l - 60.5) > 0.5:
+                logger.warning("Wrong guess maximum (%i,%i) -> (%.1f,%.1f)", i, j, k, l)
+            else:
+                logger.debug("Good guess maximum (%i,%i) -> (%.1f,%.1f)", i, j, k, l)
+                ok += 1
+        logger.debug("Success rate: %.1f", 100. * ok / self.N)
+        self.assertEqual(ok, self.N, "Maximum is always found")
+
+    def test_map(self):
+        N = 6
+        y, x = numpy.ogrid[:N,:N + 10]
+        img = x + y
+        b = BilinearImage(img)
+        x2d = numpy.zeros_like(y) + x
+        y2d = numpy.zeros_like(x) + y
+        res1 = b.map_coordinates((y2d, x2d))
+        self.assertEqual(abs(res1 - img).max(), 0, "images are the same (corners)")
+
+        x2d = numpy.zeros_like(y) + (x[:,:-1] + 0.5)
+        y2d = numpy.zeros_like(x[:,:-1]) + y
+        res1 = b.map_coordinates((y2d, x2d))
+        self.assertEqual(abs(res1 - img[:,:-1] - 0.5).max(), 0, "images are the same (middle)")
+
+        x2d = numpy.zeros_like(y[:-1,:]) + (x[:,:-1] + 0.5)
+        y2d = numpy.zeros_like(x[:,:-1]) + (y[:-1,:] + 0.5)
+        res1 = b.map_coordinates((y2d, x2d))
+        self.assertEqual(abs(res1 - img[:-1, 1:]).max(), 0, "images are the same (center)")
+
+    def test_mask_grad(self):
+        N = 100
+        img = numpy.arange(N * N).reshape(N, N)
+        # No mask on the boundaries, makes the test complicated, pixel always separated
+        masked = 2 * numpy.random.randint(0, int((N - 1) / 2), size=(2, N)) + 1
+        mask = numpy.zeros((N, N), dtype=numpy.uint8)
+        mask[(masked[0], masked[1])] = 1
+        self.assertLessEqual(mask.sum(), N, "At most N pixels are masked")
+
+        b = BilinearImage(img, mask=mask)
+        self.assertEqual(b.has_mask, True, "interpolator has mask")
+        self.assertEqual(b.maxi, N * N - 1, "maxi is N²-1")
+        self.assertEqual(b.mini, 0, "mini is 0")
+
+        y, x = numpy.ogrid[:N,:N]
+        x2d = numpy.zeros_like(y) + x
+        y2d = numpy.zeros_like(x) + y
+        res1 = b.map_coordinates((y2d, x2d))
+        self.assertEqual(numpy.nanmax(abs(res1 - img)), 0, "images are the same (corners), or Nan ")
+
+        x2d = numpy.zeros_like(y) + (x[:,:-1] + 0.5)
+        y2d = numpy.zeros_like(x[:,:-1]) + y
+        res1 = b.map_coordinates((y2d, x2d))
+        self.assertLessEqual(numpy.max(abs(res1 - img[:, 1:] + 1 / 2.)), 0.5, "images are the same (middle) +/- 0.5")
+
+        x2d = numpy.zeros_like(y[:-1]) + (x[:,:-1] + 0.5)
+        y2d = numpy.zeros_like(x[:,:-1]) + (y[:-1] + 0.5)
+        res1 = b.map_coordinates((y2d, x2d))
+        exp = 0.25 * (img[:-1,:-1] + img[:-1, 1:] + img[1:,:-1] + img[1:, 1:])
+        self.assertLessEqual(abs(res1 - exp).max(), N / 4, "images are almost the same (center)")
+
+    def test_profile_grad(self):
+        N = 100
+        img = numpy.arange(N * N).reshape(N, N)
+        b = BilinearImage(img)
+        res1 = b.profile_line((0, 0), (N - 1, N - 1))
+        l = numpy.ceil(numpy.sqrt(2) * N)
+        self.assertEqual(len(res1), l, "Profile has correct length")
+        self.assertLess((res1[:-2] - res1[1:-1]).std(), 1e-3, "profile is linear (excluding last point)")
+
+    def test_profile_gaus(self):
+        N = 100
+        x = numpy.arange(N) - N // 2.0
+        g = numpy.exp(-x * x / (N * N))
+        img = numpy.outer(g, g)
+        b = BilinearImage(img)
+        res_hor = b.profile_line((N // 2, 0), (N // 2, N - 1))
+        res_ver = b.profile_line((0, N // 2), (N - 1, N // 2))
+        self.assertEqual(len(res_hor), N, "Profile has correct length")
+        self.assertEqual(len(res_ver), N, "Profile has correct length")
+        self.assertLess(abs(res_hor - g).max(), 1e-5, "correct horizontal profile")
+        self.assertLess(abs(res_ver - g).max(), 1e-5, "correct vertical profile")
+
+        # Profile with linewidth=3
+        expected_profile = img[:, N // 2 - 1:N // 2 + 2].mean(axis=1)
+        res_hor = b.profile_line((N // 2, 0), (N // 2, N - 1), linewidth=3)
+        res_ver = b.profile_line((0, N // 2), (N - 1, N // 2), linewidth=3)
+
+        self.assertEqual(len(res_hor), N, "Profile has correct length")
+        self.assertEqual(len(res_ver), N, "Profile has correct length")
+        self.assertLess(abs(res_hor - expected_profile).max(), 1e-5,
+                        "correct horizontal profile")
+        self.assertLess(abs(res_ver - expected_profile).max(), 1e-5,
+                        "correct vertical profile")
diff --git a/src/silx/image/test/test_medianfilter.py b/src/silx/image/test/test_medianfilter.py
new file mode 100644
index 0000000..d3386a4
--- /dev/null
+++ b/src/silx/image/test/test_medianfilter.py
@@ -0,0 +1,64 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2016-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+"""Tests that the different implementation of opencl (cpp, opencl) are
+    accessible
+"""
+
+__authors__ = ["H. Payno"]
+__license__ = "MIT"
+__date__ = "11/05/2017"
+
+import unittest
+from silx.image import medianfilter
+import numpy
+
+from silx.opencl.common import ocl
+
+
+class TestMedianFilterEngines(unittest.TestCase):
+    """Make sure we have access to all the different implementation of
+    median filter from image medfilt"""
+
+
+    IMG = numpy.arange(10000.).reshape(100, 100)
+
+    KERNEL = (1, 1)
+
+    def testCppMedFilt2d(self):
+        """test cpp engine for medfilt2d"""
+        res = medianfilter.medfilt2d(
+            image=TestMedianFilterEngines.IMG,
+            kernel_size=TestMedianFilterEngines.KERNEL,
+            engine='cpp')
+        self.assertTrue(numpy.array_equal(res, TestMedianFilterEngines.IMG))
+
+    @unittest.skipUnless(ocl, "PyOpenCl is missing")
+    def testOpenCLMedFilt2d(self):
+        """test cpp engine for medfilt2d"""
+        res = medianfilter.medfilt2d(
+            image=TestMedianFilterEngines.IMG,
+            kernel_size=TestMedianFilterEngines.KERNEL,
+            engine='opencl')
+        self.assertTrue(numpy.array_equal(res, TestMedianFilterEngines.IMG))
diff --git a/src/silx/image/test/test_shapes.py b/src/silx/image/test/test_shapes.py
new file mode 100644
index 0000000..63abc00
--- /dev/null
+++ b/src/silx/image/test/test_shapes.py
@@ -0,0 +1,354 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2016 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+"""Tests for polygon functions
+"""
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "15/02/2019"
+
+
+import logging
+import unittest
+import numpy
+
+from silx.utils.testutils import ParametricTestCase
+from silx.image import shapes
+
+_logger = logging.getLogger(__name__)
+
+
+class TestPolygonFill(ParametricTestCase):
+    """basic poylgon test"""
+
+    def test_squares(self):
+        """Test polygon fill for a square polygons"""
+        mask_shape = 4, 4
+        tests = {
+            # test name: [(row min, row max), (col min, col max)]
+            'square in': [(1, 3), (1, 3)],
+            'square out': [(1, 3), (1, 10)],
+            'square around': [(-1, 5), (-1, 5)],
+            }
+
+        for test_name, (rows, cols) in tests.items():
+            with self.subTest(msg=test_name, rows=rows, cols=cols,
+                              mask_shape=mask_shape):
+                ref_mask = numpy.zeros(mask_shape, dtype=numpy.uint8)
+                ref_mask[max(0, rows[0]):rows[1],
+                         max(0, cols[0]):cols[1]] = True
+
+                vertices = [(rows[0], cols[0]), (rows[1], cols[0]),
+                            (rows[1], cols[1]), (rows[0], cols[1])]
+                mask = shapes.polygon_fill_mask(vertices, ref_mask.shape)
+                is_equal = numpy.all(numpy.equal(ref_mask, mask))
+                if not is_equal:
+                    _logger.debug('%s failed with mask != ref_mask:',
+                                  test_name)
+                    _logger.debug('result:\n%s', str(mask))
+                    _logger.debug('ref:\n%s', str(ref_mask))
+                self.assertTrue(is_equal)
+
+    def test_eight(self):
+        """Tests with eight shape with different rotation and direction"""
+        ref_mask = numpy.array((
+            (1, 1, 1, 1, 1, 0),
+            (0, 1, 1, 1, 0, 0),
+            (0, 0, 1, 0, 0, 0),
+            (0, 0, 1, 0, 0, 0),
+            (0, 1, 1, 1, 0, 0),
+            (0, 0, 0, 0, 0, 0)), dtype=numpy.uint8)
+        ref_mask_rot = numpy.asarray(numpy.logical_not(ref_mask),
+                                     dtype=numpy.uint8)
+        ref_mask_rot[:, -1] = 0
+        ref_mask_rot[-1, :] = 0
+
+        tests = {
+            'dir 1': ([(0, 0), (5, 5), (5, 0), (0, 5)], ref_mask),
+            'dir 1, rot 90': ([(5, 0), (0, 5), (5, 5), (0, 0)], ref_mask_rot),
+            'dir 1, rot 180': ([(5, 5), (0, 0), (0, 5), (5, 0)], ref_mask),
+            'dir 1, rot -90': ([(0, 5), (5, 0), (0, 0), (5, 5)], ref_mask_rot),
+            'dir 2': ([(0, 0), (0, 5), (5, 0), (5, 5)], ref_mask),
+            'dir 2, rot 90': ([(5, 0), (0, 0), (5, 5), (0, 5)], ref_mask_rot),
+            'dir 2, rot 180': ([(5, 5), (5, 0), (0, 5), (0, 0)], ref_mask),
+            'dir 2, rot -90': ([(0, 5), (5, 5), (0, 0), (5, 0)], ref_mask_rot),
+        }
+
+        for test_name, (vertices, ref_mask) in tests.items():
+            with self.subTest(msg=test_name):
+                mask = shapes.polygon_fill_mask(vertices, ref_mask.shape)
+                is_equal = numpy.all(numpy.equal(ref_mask, mask))
+                if not is_equal:
+                    _logger.debug('%s failed with mask != ref_mask:',
+                                  test_name)
+                    _logger.debug('result:\n%s', str(mask))
+                    _logger.debug('ref:\n%s', str(ref_mask))
+                self.assertTrue(is_equal)
+
+    def test_shapes(self):
+        """Tests with shapes and reference mask"""
+        tests = {
+            # name: (
+            #     polygon corners as a list of (row, col),
+            #     ref_mask)
+            'concave polygon': (
+                [(1, 1), (4, 3), (1, 5), (2, 3)],
+                numpy.array((
+                    (0, 0, 0, 0, 0, 0, 0, 0),
+                    (0, 0, 0, 0, 0, 0, 0, 0),
+                    (0, 0, 1, 1, 1, 0, 0, 0),
+                    (0, 0, 0, 1, 0, 0, 0, 0),
+                    (0, 0, 0, 0, 0, 0, 0, 0),
+                    (0, 0, 0, 0, 0, 0, 0, 0)), dtype=numpy.uint8)),
+            'concave polygon partly outside mask': (
+                [(-1, -1), (4, 3), (1, 5), (2, 3)],
+                numpy.array((
+                    (1, 0, 0, 0, 0, 0),
+                    (0, 1, 0, 0, 0, 0),
+                    (0, 0, 1, 1, 1, 0),
+                    (0, 0, 0, 1, 0, 0),
+                    (0, 0, 0, 0, 0, 0),
+                    (0, 0, 0, 0, 0, 0),
+                    (0, 0, 0, 0, 0, 0),
+                    (0, 0, 0, 0, 0, 0)), dtype=numpy.uint8)),
+            'polygon surrounding mask': (
+                [(-1, -1), (-1, 7), (7, 7), (7, -1), (0, -1),
+                 (8, -2), (8, 8), (-2, 8)],
+                numpy.zeros((6, 6), dtype=numpy.uint8))
+            }
+
+        for test_name, (vertices, ref_mask) in tests.items():
+            with self.subTest(msg=test_name):
+                mask = shapes.polygon_fill_mask(vertices, ref_mask.shape)
+                is_equal = numpy.all(numpy.equal(ref_mask, mask))
+                if not is_equal:
+                    _logger.debug('%s failed with mask != ref_mask:',
+                                  test_name)
+                    _logger.debug('result:\n%s', str(mask))
+                    _logger.debug('ref:\n%s', str(ref_mask))
+                self.assertTrue(is_equal)
+
+
+class TestDrawLine(ParametricTestCase):
+    """basic draw line test"""
+
+    def test_aligned_lines(self):
+        """Test drawing horizontal, vertical and diagonal lines"""
+
+        lines = {  # test_name: (drow, dcol)
+            'Horizontal line, col0 < col1': (0, 10),
+            'Horizontal line, col0 > col1': (0, -10),
+            'Vertical line, row0 < row1': (10, 0),
+            'Vertical line, row0 > row1': (-10, 0),
+            'Diagonal col0 < col1 and row0 < row1': (10, 10),
+            'Diagonal col0 < col1 and row0 > row1': (-10, 10),
+            'Diagonal col0 > col1 and row0 < row1': (10, -10),
+            'Diagonal col0 > col1 and row0 > row1': (-10, -10),
+        }
+        row0, col0 = 1, 2  # Start point
+
+        for test_name, (drow, dcol) in lines.items():
+            row1 = row0 + drow
+            col1 = col0 + dcol
+            with self.subTest(msg=test_name, drow=drow, dcol=dcol):
+                # Build reference coordinates from drow and dcol
+                if drow == 0:
+                    rows = row0 * numpy.ones(abs(dcol) + 1)
+                else:
+                    step = 1 if drow > 0 else -1
+                    rows = row0 + numpy.arange(0, drow + step, step)
+
+                if dcol == 0:
+                    cols = col0 * numpy.ones(abs(drow) + 1)
+                else:
+                    step = 1 if dcol > 0 else -1
+                    cols = col0 + numpy.arange(0, dcol + step, step)
+                ref_coords = rows, cols
+
+                result = shapes.draw_line(row0, col0, row1, col1)
+                self.assertTrue(self.isEqual(test_name, result, ref_coords))
+
+    def test_noline(self):
+        """Test pt0 == pt1"""
+        for width in range(4):
+            with self.subTest(width=width):
+                result = shapes.draw_line(1, 2, 1, 2, width)
+                self.assertTrue(numpy.all(numpy.equal(result, [(1,), (2,)])))
+
+    def test_lines(self):
+        """Test lines not aligned with axes for 8 slopes and directions"""
+        row0, col0 = 1, 1
+
+        dy, dx = 3, 5
+        ref_coords = numpy.array(
+            [(0, 0), (1, 1), (1, 2), (2, 3), (2, 4), (3, 5)])
+
+        # Build lines for the 8 octants from this coordinantes
+        lines = {  # name: (drow, dcol, ref_coords)
+            '1st octant': (dy, dx, ref_coords),
+            '2nd octant': (dx, dy, ref_coords[:, (1, 0)]),  # invert x and y
+            '3rd octant': (dx, -dy, ref_coords[:, (1, 0)] * (1, -1)),
+            '4th octant': (dy, -dx, ref_coords * (1, -1)),
+            '5th octant': (-dy, -dx, ref_coords * (-1, -1)),
+            '6th octant': (-dx, -dy, ref_coords[:, (1, 0)] * (-1, -1)),
+            '7th octant': (-dx, dy, ref_coords[:, (1, 0)] * (-1, 1)),
+            '8th octant': (-dy, dx, ref_coords * (-1, 1))
+        }
+
+        # Test with different starting points with positive and negative coords
+        for row0, col0 in ((0, 0), (2, 3), (-4, 1), (-5, -6), (8, -7)):
+            for name, (drow, dcol, ref_coords) in lines.items():
+                row1 = row0 + drow
+                col1 = col0 + dcol
+                # Transpose from ((row0, col0), ...) to (rows, cols)
+                ref_coords = numpy.transpose(ref_coords + (row0, col0))
+
+                with self.subTest(msg=name,
+                                  pt0=(row0, col0), pt1=(row1, col1)):
+                    result = shapes.draw_line(row0, col0, row1, col1)
+                    self.assertTrue(self.isEqual(name, result, ref_coords))
+
+    def test_width(self):
+        """Test of line width"""
+
+        lines = {  # test_name: row0, col0, row1, col1, width, ref
+            'horizontal w=2':
+                (0, 0, 0, 1, 2, ((0, 1, 0, 1),
+                                 (0, 0, 1, 1))),
+            'horizontal w=3':
+                (0, 0, 0, 1, 3, ((-1, 0, 1, -1, 0, 1),
+                                 (0, 0, 0, 1, 1, 1))),
+            'vertical w=2':
+                (0, 0, 1, 0, 2, ((0, 0, 1, 1),
+                                 (0, 1, 0, 1))),
+            'vertical w=3':
+                (0, 0, 1, 0, 3, ((0, 0, 0, 1, 1, 1),
+                                 (-1, 0, 1, -1, 0, 1))),
+            'diagonal w=3':
+                (0, 0, 1, 1, 3, ((-1, 0, 1, 0, 1, 2),
+                                 (0, 0, 0, 1, 1, 1))),
+            '1st octant w=3':
+                (0, 0, 1, 2, 3,
+                 numpy.array(((-1, 0), (0, 0), (1, 0),
+                              (0, 1), (1, 1), (2, 1),
+                              (0, 2), (1, 2), (2, 2))).T),
+            '2nd octant w=3':
+                (0, 0, 2, 1, 3,
+                 numpy.array(((0, -1), (0, 0), (0, 1),
+                              (1, 0), (1, 1), (1, 2),
+                              (2, 0), (2, 1), (2, 2))).T),
+        }
+
+        for test_name, (row0, col0, row1, col1, width, ref) in lines.items():
+            with self.subTest(msg=test_name,
+                              pt0=(row0, col0), pt1=(row1, col1), width=width):
+                result = shapes.draw_line(row0, col0, row1, col1, width)
+                self.assertTrue(self.isEqual(test_name, result, ref))
+
+    def isEqual(self, test_name, result, ref):
+        """Test equality of two numpy arrays and log them if different"""
+        is_equal = numpy.all(numpy.equal(result, ref))
+        if not is_equal:
+            _logger.debug('%s failed with result != ref:',
+                          test_name)
+            _logger.debug('result:\n%s', str(result))
+            _logger.debug('ref:\n%s', str(ref))
+        return is_equal
+
+
+class TestCircleFill(ParametricTestCase):
+    """Tests for circle filling"""
+
+    def testCircle(self):
+        """Test circle_fill with different input parameters"""
+
+        square3x3 = numpy.array(((-1, -1, -1, 0, 0, 0, 1, 1, 1),
+                                 (-1, 0, 1, -1, 0, 1, -1, 0, 1)))
+
+        tests = [
+            # crow, ccol, radius, ref_coords = (ref_rows, ref_cols)
+            (0, 0, 1, ((0,), (0,))),
+            (10, 15, 1, ((10,), (15,))),
+            (0, 0, 1.5, square3x3),
+            (5, 10, 2, (5 + square3x3[0], 10 + square3x3[1])),
+            (10, 20, 3.5, (
+                10 + numpy.array((-3, -3, -3,
+                                  -2, -2, -2, -2, -2,
+                                  -1, -1, -1, -1, -1, -1, -1,
+                                  0, 0, 0, 0, 0, 0, 0,
+                                  1, 1, 1, 1, 1, 1, 1,
+                                  2, 2, 2, 2, 2,
+                                  3, 3, 3)),
+                20 + numpy.array((-1, 0, 1,
+                                  -2, -1, 0, 1, 2,
+                                  -3, -2, -1, 0, 1, 2, 3,
+                                  -3, -2, -1, 0, 1, 2, 3,
+                                  -3, -2, -1, 0, 1, 2, 3,
+                                  -2, -1, 0, 1, 2,
+                                  -1, 0, 1)))),
+        ]
+
+        for crow, ccol, radius, ref_coords in tests:
+            with self.subTest(crow=crow, ccol=ccol, radius=radius):
+                coords = shapes.circle_fill(crow, ccol, radius)
+                is_equal = numpy.all(numpy.equal(coords, ref_coords))
+                if not is_equal:
+                    _logger.debug('result:\n%s', str(coords))
+                    _logger.debug('ref:\n%s', str(ref_coords))
+                self.assertTrue(is_equal)
+
+
+class TestEllipseFill(unittest.TestCase):
+    """Tests for ellipse filling"""
+
+    def testPoint(self):
+        args = [1, 1, 1, 1]
+        result = shapes.ellipse_fill(*args)
+        expected = numpy.array(([1], [1]))
+        numpy.testing.assert_array_equal(result, expected)
+
+    def testTranslatedPoint(self):
+        args = [10, 10, 1, 1]
+        result = shapes.ellipse_fill(*args)
+        expected = numpy.array(([10], [10]))
+        numpy.testing.assert_array_equal(result, expected)
+
+    def testEllipse(self):
+        args = [0, 0, 20, 10]
+        rows, cols = shapes.ellipse_fill(*args)
+        self.assertEqual(len(rows), 617)
+        self.assertEqual(rows.mean(), 0)
+        self.assertAlmostEqual(rows.std(), 10.025575, places=3)
+        self.assertEqual(len(cols), 617)
+        self.assertEqual(cols.mean(), 0)
+        self.assertAlmostEqual(cols.std(), 4.897325, places=3)
+
+    def testTranslatedEllipse(self):
+        args = [0, 0, 20, 10]
+        expected_rows, expected_cols = shapes.ellipse_fill(*args)
+        args = [10, 50, 20, 10]
+        rows, cols = shapes.ellipse_fill(*args)
+        numpy.testing.assert_allclose(rows, expected_rows + 10)
+        numpy.testing.assert_allclose(cols, expected_cols + 50)
diff --git a/src/silx/image/test/test_tomography.py b/src/silx/image/test/test_tomography.py
new file mode 100644
index 0000000..f391a72
--- /dev/null
+++ b/src/silx/image/test/test_tomography.py
@@ -0,0 +1,54 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2016-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+"""
+Tests that the functions of tomography are valid
+"""
+
+__authors__ = ["H. Payno"]
+__license__ = "MIT"
+__date__ = "12/09/2017"
+
+import unittest
+import numpy
+from silx.test.utils import utilstest
+from silx.image import tomography
+
+class TestTomography(unittest.TestCase):
+    """
+
+    """
+
+    def setUp(self):
+        self.sinoTrueData = numpy.load(utilstest.getfile("sino500.npz"))["data"]
+
+    def testCalcCenterCentroid(self):
+        centerTD = tomography.calc_center_centroid(self.sinoTrueData)
+        self.assertTrue(numpy.isclose(centerTD, 256, rtol=0.01))
+
+    def testCalcCenterCorr(self):
+        centerTrueData = tomography.calc_center_corr(self.sinoTrueData,
+                                                     fullrot=False,
+                                                     props=1)
+        self.assertTrue(numpy.isclose(centerTrueData, 256, rtol=0.01))
diff --git a/src/silx/image/tomography.py b/src/silx/image/tomography.py
new file mode 100644
index 0000000..53855c1
--- /dev/null
+++ b/src/silx/image/tomography.py
@@ -0,0 +1,301 @@
+# coding: utf-8
+# /*##########################################################################
+# Copyright (C) 2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+"""
+This module contains utilitary functions for tomography
+"""
+
+__author__ = ["P. Paleo"]
+__license__ = "MIT"
+__date__ = "12/09/2017"
+
+
+import numpy as np
+from math import pi
+from functools import lru_cache
+from itertools import product
+from bisect import bisect
+from silx.math.fit import leastsq
+
+# ------------------------------------------------------------------------------
+# -------------------- Filtering-related functions -----------------------------
+# ------------------------------------------------------------------------------
+
+def compute_ramlak_filter(dwidth_padded, dtype=np.float32):
+    """
+    Compute the Ramachandran-Lakshminarayanan (Ram-Lak) filter, used in
+    filtered backprojection.
+
+    :param dwidth_padded: width of the 2D sinogram after padding
+    :param dtype: data type
+    """
+    L = dwidth_padded
+    h = np.zeros(L, dtype=dtype)
+    L2 = L//2+1
+    h[0] = 1/4.
+    j = np.linspace(1, L2, L2//2, False).astype(dtype) # np < 1.9.0
+    h[1:L2:2] = -1./(pi**2 * j**2)
+    h[L2:] = np.copy(h[1:L2-1][::-1])
+    return h
+
+
+def tukey(N, alpha=0.5):
+    """
+    Compute the Tukey apodization window.
+
+    :param int N: Number of points.
+    :param float alpha:
+    """
+    apod = np.zeros(N)
+    x = np.arange(N)/(N-1)
+    r = alpha
+    M1 = (0 <= x) * (x < r/2)
+    M2 = (r/2 <= x) * (x <= 1 - r/2)
+    M3 = (1 - r/2 < x) * (x <= 1)
+    apod[M1] = (1 + np.cos(2*pi/r * (x[M1] - r/2)))/2.
+    apod[M2] = 1.
+    apod[M3] = (1 + np.cos(2*pi/r * (x[M3] - 1 + r/2)))/2.
+    return apod
+
+
+def lanczos(N):
+    """
+    Compute the Lanczos window (truncated sinc) of width N.
+
+    :param int N: window width
+    """
+    x = np.arange(N)/(N-1)
+    return np.sin(pi*(2*x-1))/(pi*(2*x-1))
+
+
+def compute_fourier_filter(dwidth_padded, filter_name, cutoff=1.):
+    """
+    Compute the filter used for FBP.
+
+    :param dwidth_padded: padded detector width. As the filtering is done by the
+        Fourier convolution theorem, dwidth_padded should be at least 2*dwidth.
+    :param filter_name: Name of the filter. Available filters are:
+        Ram-Lak, Shepp-Logan, Cosine, Hamming, Hann, Tukey, Lanczos.
+    :param cutoff: Cut-off frequency, if relevant.
+    """
+    Nf = dwidth_padded
+    #~ filt_f = np.abs(np.fft.fftfreq(Nf))
+    rl = compute_ramlak_filter(Nf, dtype=np.float64)
+    filt_f = np.fft.fft(rl)
+
+    filter_name = filter_name.lower()
+    if filter_name in ["ram-lak", "ramlak"]:
+        return filt_f
+
+    w = 2 * pi * np.fft.fftfreq(dwidth_padded)
+    d = cutoff
+    apodization = {
+        # ~OK
+        "shepp-logan": np.sin(w[1:Nf]/(2*d))/(w[1:Nf]/(2*d)),
+        # ~OK
+        "cosine": np.cos(w[1:Nf]/(2*d)),
+        # OK
+        "hamming": 0.54*np.ones_like(filt_f)[1:Nf] + .46 * np.cos(w[1:Nf]/d),
+        # OK
+        "hann": (np.ones_like(filt_f)[1:Nf] + np.cos(w[1:Nf]/d))/2.,
+        # These one is not compatible with Astra - TODO investigate why
+        "tukey": np.fft.fftshift(tukey(dwidth_padded, alpha=d/2.))[1:Nf],
+        "lanczos": np.fft.fftshift(lanczos(dwidth_padded))[1:Nf],
+    }
+    if filter_name not in apodization:
+        raise ValueError("Unknown filter %s. Available filters are %s" %
+                         (filter_name, str(apodization.keys())))
+    filt_f[1:Nf] *= apodization[filter_name]
+    return filt_f
+
+
+@lru_cache(maxsize=1)
+def generate_powers():
+    """
+    Generate a list of powers of [2, 3, 5, 7],
+    up to (2**15)*(3**9)*(5**6)*(7**5).
+    """
+    primes = [2, 3, 5, 7]
+    maxpow = {2: 15, 3: 9, 5: 6, 7: 5}
+    valuations = []
+    for prime in primes:
+        # disallow any odd number (for R2C transform), and any number
+        # not multiple of 4 (Ram-Lak filter behaves strangely when
+        # dwidth_padded/2 is not even)
+        minval = 2 if prime == 2 else 0
+        valuations.append(range(minval, maxpow[prime]+1))
+    powers = product(*valuations)
+    res = []
+    for pw in powers:
+        res.append(np.prod(list(map(lambda x : x[0]**x[1], zip(primes, pw)))))
+    return np.unique(res)
+
+
+def get_next_power(n, powers=None):
+    """
+    Given a number, get the closest (upper) number p such that
+    p is a power of 2, 3, 5 and 7.
+    """
+    if powers is None:
+        powers = generate_powers()
+    idx = bisect(powers, n)
+    if powers[idx-1] == n:
+        return n
+    return powers[idx]
+
+
+# ------------------------------------------------------------------------------
+# ------------- Functions for determining the center of rotation  --------------
+# ------------------------------------------------------------------------------
+
+
+
+def calc_center_corr(sino, fullrot=False, props=1):
+    """
+    Compute a guess of the Center of Rotation (CoR) of a given sinogram.
+    The computation is based on the correlation between the line projections at
+    angle (theta = 0) and at angle (theta = 180).
+
+    Note that for most scans, the (theta=180) angle is not included,
+    so the CoR might be underestimated.
+    In a [0, 360[ scan, the projection angle at (theta=180) is exactly in the
+    middle for odd number of projections.
+
+    :param numpy.ndarray sino: Sinogram
+    :param bool fullrot: optional. If False (default), the scan is assumed to
+                         be [0, 180).
+                         If True, the scan is assumed to be [0, 380).
+    :param int props: optional. Number of propositions for the CoR
+    """
+
+    n_a, n_d = sino.shape
+    first = 0
+    last = -1 if not(fullrot) else n_a // 2
+    proj1 = sino[first, :]
+    proj2 = sino[last, :][::-1]
+
+    # Compute the correlation in the Fourier domain
+    proj1_f = np.fft.fft(proj1, 2 * n_d)
+    proj2_f = np.fft.fft(proj2, 2 * n_d)
+    corr = np.abs(np.fft.ifft(proj1_f * proj2_f.conj()))
+
+    if props == 1:
+        pos = np.argmax(corr)
+        if pos > n_d // 2:
+            pos -= n_d
+        return (n_d + pos) / 2.
+    else:
+        corr_argsorted = np.argsort(corr)[:props]
+        corr_argsorted[corr_argsorted > n_d // 2] -= n_d
+        return (n_d + corr_argsorted) / 2.
+
+
+def _sine_function(t, offset, amplitude, phase):
+    """
+    Helper function for calc_center_centroid
+    """
+    n_angles = t.shape[0]
+    res = amplitude * np.sin(2 * pi * (1. / (2 * n_angles)) * t + phase)
+    return offset + res
+
+
+def _sine_function_derivative(t, params, eval_idx):
+    """
+    Helper function for calc_center_centroid
+    """
+    offset, amplitude, phase = params
+    n_angles = t.shape[0]
+    w = 2.0 * pi * (1. / (2.0 * n_angles)) * t + phase
+    grad = (1.0, np.sin(w), amplitude*np.cos(w))
+    return grad[eval_idx]
+
+
+def calc_center_centroid(sino):
+    """
+    Compute a guess of the Center of Rotation (CoR) of a given sinogram.
+    The computation is based on the computation of the centroid of each
+    projection line, which should be a sine function according to the
+    Helgason-Ludwig condition.
+    This method is unlikely to work in local tomography.
+
+    :param numpy.ndarray sino: Sinogram
+    """
+
+    n_a, n_d = sino.shape
+    # Compute the vector of centroids of the sinogram
+    i = np.arange(n_d)
+    centroids = np.sum(sino*i, axis=1)/np.sum(sino, axis=1)
+
+    # Fit with a sine function : phase, amplitude, offset
+    # Using non-linear Levenberg–Marquardt  algorithm
+    angles = np.linspace(0, n_a, n_a, True)
+    # Initial parameter vector
+    cmax, cmin = centroids.max(), centroids.min()
+    offs = (cmax + cmin) / 2.
+    amp = (cmax - cmin) / 2.
+    phi = 1.1
+    p0 = (offs, amp, phi)
+
+    constraints = np.zeros((3, 3))
+
+    popt, _ = leastsq(model=_sine_function,
+                      xdata=angles,
+                      ydata=centroids,
+                      p0=p0,
+                      sigma=None,
+                      constraints=constraints,
+                      model_deriv=None,
+                      epsfcn=None,
+                      deltachi=None,
+                      full_output=0,
+                      check_finite=True,
+                      left_derivative=False,
+                      max_iter=100)
+    return popt[0]
+
+
+
+# ------------------------------------------------------------------------------
+# -------------------- Visualization-related functions -------------------------
+# ------------------------------------------------------------------------------
+
+
+def rescale_intensity(img, from_subimg=None, percentiles=None):
+    """
+    clamp intensity into the [2, 98] percentiles
+
+    :param img:
+    :param from_subimg:
+    :param percentiles:
+    :return: the rescale intensity
+    """
+    if percentiles is None:
+        percentiles = [2, 98]
+    else:
+        assert type(percentiles) in (tuple, list)
+        assert(len(percentiles) == 2)
+    data = from_subimg if from_subimg is not None else img
+    imin, imax = np.percentile(data, percentiles)
+    res = np.clip(img, imin, imax)
+    return res
+
diff --git a/src/silx/image/utils.py b/src/silx/image/utils.py
new file mode 100644
index 0000000..996d010
--- /dev/null
+++ b/src/silx/image/utils.py
@@ -0,0 +1,53 @@
+# -*- coding: utf-8 -*-
+# /*##########################################################################
+# Copyright (C) 2019 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ############################################################################*/
+
+import numpy as np
+from math import ceil
+
+def gaussian_kernel(sigma, cutoff=4, force_odd_size=False):
+    """
+    Generates a Gaussian convolution kernel.
+
+    :param sigma: Standard Deviation of the Gaussian curve.
+    :param cutoff: Parameter tuning the truncation of the Gaussian.
+        The higher cutoff, the biggest the array will be (and the closest to
+        a "true" Gaussian function).
+    :param force_odd_size: when set to True, the resulting array will always
+        have an odd size, regardless of the values of "sigma" and "cutoff".
+    :return: a numpy.ndarray containing the truncated Gaussian function.
+        The array size is 2*c*s+1 where c=cutoff, s=sigma.
+
+    Nota: due to the quick decay of the Gaussian function, small values of the
+        "cutoff" parameter are usually fine. The energy difference between a
+        Gaussian truncated to [-c, c] and a "true" one is
+            erfc(c/(sqrt(2)*s))
+        so choosing cutoff=4*sigma keeps the truncation error below 1e-4.
+    """
+    size = int(ceil(2 * cutoff * sigma + 1))
+    if force_odd_size and size % 2 == 0:
+        size += 1
+    x = np.arange(size) - (size - 1.0) / 2.0
+    g = np.exp(-(x / sigma) ** 2 / 2.0)
+    g /= g.sum()
+    return g