New upstream version 0.6.0+dfsg

5 files changed, 697 insertions, 4 deletions

diff --git a/silx/opencl/test/__init__.py b/silx/opencl/test/__init__.py
index 24aa06e..f6aadcd 100644
--- a/silx/opencl/test/__init__.py
+++ b/silx/opencl/test/__init__.py
@@ -24,18 +24,31 @@
 
 __authors__ = ["J. Kieffer"]
 __license__ = "MIT"
-__date__ = "15/03/2017"
+__date__ = "01/09/2017"
 
+import os
 import unittest
 from . import test_addition
 from . import test_medfilt
-from ..sift import test as test_sift
-
+from . import test_backprojection
+from . import test_projection
+from . import test_linalg
+from . import test_array_utils
 
 def suite():
     test_suite = unittest.TestSuite()
     test_suite.addTests(test_addition.suite())
     test_suite.addTests(test_medfilt.suite())
-    test_suite.addTests(test_sift.suite())
+    test_suite.addTests(test_backprojection.suite())
+    test_suite.addTests(test_projection.suite())
+    test_suite.addTests(test_linalg.suite())
+    test_suite.addTests(test_array_utils.suite())
+
+    # Allow to remove sift from the project
+    test_base_dir = os.path.dirname(__file__)
+    sift_dir = os.path.join(test_base_dir, "..", "sift")
+    if os.path.exists(sift_dir):
+        from ..sift import test as test_sift
+        test_suite.addTests(test_sift.suite())
 
     return test_suite
diff --git a/silx/opencl/test/test_array_utils.py b/silx/opencl/test/test_array_utils.py
new file mode 100644
index 0000000..833d828
--- /dev/null
+++ b/silx/opencl/test/test_array_utils.py
@@ -0,0 +1,161 @@
+#!/usr/bin/env python
+# 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.
+#
+# ###########################################################################*/
+"""Test of the OpenCL array_utils"""
+
+from __future__ import division, print_function
+
+__authors__ = ["Pierre paleo"]
+__license__ = "MIT"
+__copyright__ = "2013-2017 European Synchrotron Radiation Facility, Grenoble, France"
+__date__ = "14/06/2017"
+
+
+import time
+import logging
+import numpy as np
+import unittest
+try:
+    import mako
+except ImportError:
+    mako = None
+from ..common import ocl
+if ocl:
+    import pyopencl as cl
+    import pyopencl.array as parray
+    from .. import linalg
+from ..utils import get_opencl_code
+from silx.test.utils import utilstest
+
+logger = logging.getLogger(__name__)
+try:
+    from scipy.ndimage.filters import laplace
+    _has_scipy = True
+except ImportError:
+    _has_scipy = False
+
+
+
+@unittest.skipUnless(ocl and mako, "PyOpenCl is missing")
+class TestCpy2d(unittest.TestCase):
+
+    def setUp(self):
+        if ocl is None:
+            return
+        self.ctx = ocl.create_context()
+        if logger.getEffectiveLevel() <= logging.INFO:
+            self.PROFILE = True
+            self.queue = cl.CommandQueue(
+                            self.ctx,
+                            properties=cl.command_queue_properties.PROFILING_ENABLE)
+        else:
+            self.PROFILE = False
+            self.queue = cl.CommandQueue(self.ctx)
+        self.allocate_arrays()
+        self.program = cl.Program(self.ctx, get_opencl_code("array_utils")).build()
+
+    def allocate_arrays(self):
+        """
+        Allocate various types of arrays for the tests
+        """
+        self.prng_state = np.random.get_state()
+        # Generate arrays of random shape
+        self.shape1 = np.random.randint(20, high=512, size=(2,))
+        self.shape2 = np.random.randint(20, high=512, size=(2,))
+        self.array1 = np.random.rand(*self.shape1).astype(np.float32)
+        self.array2 = np.random.rand(*self.shape2).astype(np.float32)
+        self.d_array1 = parray.to_device(self.queue, self.array1)
+        self.d_array2 = parray.to_device(self.queue, self.array2)
+        # Generate random offsets
+        offset1_y = np.random.randint(2, high=min(self.shape1[0], self.shape2[0]) - 10)
+        offset1_x = np.random.randint(2, high=min(self.shape1[1], self.shape2[1]) - 10)
+        offset2_y = np.random.randint(2, high=min(self.shape1[0], self.shape2[0]) - 10)
+        offset2_x = np.random.randint(2, high=min(self.shape1[1], self.shape2[1]) - 10)
+        self.offset1 = (offset1_y, offset1_x)
+        self.offset2 = (offset2_y, offset2_x)
+        # Compute the size of the rectangle to transfer
+        size_y = np.random.randint(2, high=min(self.shape1[0], self.shape2[0]) - max(offset1_y, offset2_y) + 1)
+        size_x = np.random.randint(2, high=min(self.shape1[1], self.shape2[1]) - max(offset1_x, offset2_x) + 1)
+        self.transfer_shape = (size_y, size_x)
+
+    def tearDown(self):
+        self.array1 = None
+        self.array2 = None
+        self.d_array1.data.release()
+        self.d_array2.data.release()
+        self.d_array1 = None
+        self.d_array2 = None
+        self.ctx = None
+        self.queue = None
+
+    def compare(self, result, reference):
+        errmax = np.max(np.abs(result - reference))
+        logger.info("Max error = %e" % (errmax))
+        self.assertTrue(errmax == 0, str("Max error is too high"))#. PRNG state was %s" % str(self.prng_state)))
+
+    @unittest.skipUnless(ocl and mako, "pyopencl is missing")
+    def test_cpy2d(self):
+        """
+        Test rectangular transfer of self.d_array1 to self.d_array2
+        """
+        # Reference
+        o1 = self.offset1
+        o2 = self.offset2
+        T = self.transfer_shape
+        logger.info("""Testing D->D rectangular copy with (N1_y, N1_x) = %s,
+                    (N2_y, N2_x) = %s:
+                    array2[%d:%d, %d:%d] = array1[%d:%d, %d:%d]""" %
+                        (
+                            str(self.shape1), str(self.shape2),
+                            o2[0], o2[0] + T[0],
+                            o2[1], o2[1] + T[1],
+                            o1[0], o1[0] + T[0],
+                            o1[1], o1[1] + T[1]
+                        )
+                    )
+        self.array2[o2[0]:o2[0] + T[0], o2[1]:o2[1] + T[1]] = self.array1[o1[0]:o1[0] + T[0], o1[1]:o1[1] + T[1]]
+        kernel_args = (
+            self.d_array2.data,
+            self.d_array1.data,
+            np.int32(self.shape2[1]),
+            np.int32(self.shape1[1]),
+            np.int32(self.offset2[::-1]),
+            np.int32(self.offset1[::-1]),
+            np.int32(self.transfer_shape[::-1])
+        )
+        wg = None
+        ndrange = self.transfer_shape[::-1]
+        self.program.cpy2d(self.queue, ndrange, wg, *kernel_args)
+        res = self.d_array2.get()
+        self.compare(res, self.array2)
+
+
+def suite():
+    testSuite = unittest.TestSuite()
+    testSuite.addTest(TestCpy2d("test_cpy2d"))
+    return testSuite
+
+if __name__ == '__main__':
+    unittest.main(defaultTest="suite")
diff --git a/silx/opencl/test/test_backprojection.py b/silx/opencl/test/test_backprojection.py
new file mode 100644
index 0000000..342bd2f
--- /dev/null
+++ b/silx/opencl/test/test_backprojection.py
@@ -0,0 +1,165 @@
+#!/usr/bin/env python
+# 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.
+#
+# ###########################################################################*/
+"""Test of the filtered backprojection module"""
+
+from __future__ import division, print_function
+
+__authors__ = ["Pierre paleo"]
+__license__ = "MIT"
+__copyright__ = "2013-2017 European Synchrotron Radiation Facility, Grenoble, France"
+__date__ = "05/10/2017"
+
+
+import time
+import logging
+import numpy
+import unittest
+try:
+    import mako
+except ImportError:
+    mako = None
+from ..common import ocl
+if ocl:
+    from .. import backprojection
+from silx.test.utils import utilstest
+
+logger = logging.getLogger(__name__)
+
+
+def generate_coords(img_shp, center=None):
+    """
+    Return two 2D arrays containing the indexes of an image.
+    The zero is at the center of the image.
+    """
+    l_r, l_c = float(img_shp[0]), float(img_shp[1])
+    R, C = numpy.mgrid[:l_r, :l_c]
+    if center is None:
+        center0, center1 = l_r / 2., l_c / 2.
+    else:
+        center0, center1 = center
+    R = R + 0.5 - center0
+    C = C + 0.5 - center1
+    return R, C
+
+
+def clip_circle(img, center=None, radius=None):
+    """
+    Puts zeros outside the inscribed circle of the image support.
+    """
+    R, C = generate_coords(img.shape, center)
+    M = R * R + C * C
+    res = numpy.zeros_like(img)
+    if radius is None:
+        radius = img.shape[0] / 2. - 1
+    mask = M < radius * radius
+    res[mask] = img[mask]
+    return res
+
+
+@unittest.skipUnless(ocl and mako, "PyOpenCl is missing")
+class TestFBP(unittest.TestCase):
+
+    def setUp(self):
+        if ocl is None:
+            return
+        # ~ if sys.platform.startswith('darwin'):
+            # ~ self.skipTest("Backprojection is not implemented on CPU for OS X yet")
+        self.getfiles()
+        self.fbp = backprojection.Backprojection(self.sino.shape, profile=True)
+        if self.fbp.compiletime_workgroup_size < 16:
+            self.skipTest("Current implementation of OpenCL backprojection is not supported on this platform yet")
+
+    def tearDown(self):
+        self.sino = None
+#         self.fbp.log_profile()
+        self.fbp = None
+
+    def getfiles(self):
+        # load sinogram of 512x512 MRI phantom
+        self.sino = numpy.load(utilstest.getfile("sino500.npz"))["data"]
+        # load reconstruction made with ASTRA FBP (with filter designed in spatial domain)
+        self.reference_rec = numpy.load(utilstest.getfile("rec_astra_500.npz"))["data"]
+
+    def measure(self):
+        "Common measurement of timings"
+        t1 = time.time()
+        try:
+            result = self.fbp.filtered_backprojection(self.sino)
+        except RuntimeError as msg:
+            logger.error(msg)
+            return
+        t2 = time.time()
+        return t2 - t1, result
+
+    def compare(self, res):
+        """
+        Compare a result with the reference reconstruction.
+        Only the valid reconstruction zone (inscribed circle) is taken into
+        account
+        """
+        res_clipped = clip_circle(res)
+        ref_clipped = clip_circle(self.reference_rec)
+        delta = abs(res_clipped - ref_clipped)
+        bad = delta > 1
+#         numpy.save("/tmp/bad.npy", bad.astype(int))
+        logger.debug("Absolute difference: %s with %s outlier pixels out of %s", delta.max(), bad.sum(), numpy.prod(bad.shape))
+        return delta.max()
+
+    @unittest.skipUnless(ocl and mako, "pyopencl is missing")
+    def test_fbp(self):
+        """
+        tests FBP
+        """
+        # Test single reconstruction
+        # --------------------------
+        t, res = self.measure()
+        if t is None:
+            logger.info("test_fp: skipped")
+        else:
+            logger.info("test_backproj: time = %.3fs" % t)
+            err = self.compare(res)
+            msg = str("Max error = %e" % err)
+            logger.info(msg)
+            # TODO: cannot do better than 1e0 ?
+            # The plain backprojection was much better, so it must be an issue in the filtering process
+            self.assertTrue(err < 1., "Max error is too high")
+        # Test multiple reconstructions
+        # -----------------------------
+        res0 = numpy.copy(res)
+        for i in range(10):
+            res = self.fbp.filtered_backprojection(self.sino)
+            errmax = numpy.max(numpy.abs(res - res0))
+            self.assertTrue(errmax < 1.e-6, "Max error is too high")
+
+
+def suite():
+    testSuite = unittest.TestSuite()
+    testSuite.addTest(TestFBP("test_fbp"))
+    return testSuite
+
+
+if __name__ == '__main__':
+    unittest.main(defaultTest="suite")
diff --git a/silx/opencl/test/test_linalg.py b/silx/opencl/test/test_linalg.py
new file mode 100644
index 0000000..7d03983
--- /dev/null
+++ b/silx/opencl/test/test_linalg.py
@@ -0,0 +1,215 @@
+#!/usr/bin/env python
+# 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.
+#
+# ###########################################################################*/
+"""Test of the linalg module"""
+
+from __future__ import division, print_function
+
+__authors__ = ["Pierre paleo"]
+__license__ = "MIT"
+__copyright__ = "2013-2017 European Synchrotron Radiation Facility, Grenoble, France"
+__date__ = "14/06/2017"
+
+
+import time
+import logging
+import numpy as np
+import unittest
+try:
+    import mako
+except ImportError:
+    mako = None
+from ..common import ocl
+if ocl:
+    import pyopencl as cl
+    import pyopencl.array as parray
+    from .. import linalg
+from silx.test.utils import utilstest
+
+logger = logging.getLogger(__name__)
+try:
+    from scipy.ndimage.filters import laplace
+    _has_scipy = True
+except ImportError:
+    _has_scipy = False
+
+
+# TODO move this function in math or image ?
+def gradient(img):
+    '''
+    Compute the gradient of an image as a numpy array
+    Code from https://github.com/emmanuelle/tomo-tv/
+    '''
+    shape = [img.ndim, ] + list(img.shape)
+    gradient = np.zeros(shape, dtype=img.dtype)
+    slice_all = [0, slice(None, -1),]
+    for d in range(img.ndim):
+        gradient[slice_all] = np.diff(img, axis=d)
+        slice_all[0] = d + 1
+        slice_all.insert(1, slice(None))
+    return gradient
+
+
+# TODO move this function in math or image ?
+def divergence(grad):
+    '''
+    Compute the divergence of a gradient
+    Code from https://github.com/emmanuelle/tomo-tv/
+    '''
+    res = np.zeros(grad.shape[1:])
+    for d in range(grad.shape[0]):
+        this_grad = np.rollaxis(grad[d], d)
+        this_res = np.rollaxis(res, d)
+        this_res[:-1] += this_grad[:-1]
+        this_res[1:-1] -= this_grad[:-2]
+        this_res[-1] -= this_grad[-2]
+    return res
+
+
+@unittest.skipUnless(ocl and mako, "PyOpenCl is missing")
+class TestLinAlg(unittest.TestCase):
+
+    def setUp(self):
+        if ocl is None:
+            return
+        self.getfiles()
+        self.la = linalg.LinAlg(self.image.shape)
+        self.allocate_arrays()
+
+    def allocate_arrays(self):
+        """
+        Allocate various types of arrays for the tests
+        """
+        # numpy images
+        self.grad = np.zeros(self.image.shape, dtype=np.complex64)
+        self.grad2 = np.zeros((2,) + self.image.shape, dtype=np.float32)
+        self.grad_ref = gradient(self.image)
+        self.div_ref = divergence(self.grad_ref)
+        self.image2 = np.zeros_like(self.image)
+        # Device images
+        self.gradient_parray = parray.zeros(self.la.queue, self.image.shape, np.complex64)
+        # we should be using cl.Buffer(self.la.ctx, cl.mem_flags.READ_WRITE, size=self.image.nbytes*2),
+        # but platforms not suporting openCL 1.2 have a problem with enqueue_fill_buffer,
+        # so we use the parray "fill" utility
+        self.gradient_buffer = self.gradient_parray.data
+        # Do the same for image
+        self.image_parray = parray.to_device(self.la.queue, self.image)
+        self.image_buffer = self.image_parray.data
+        # Refs
+        tmp = np.zeros(self.image.shape, dtype=np.complex64)
+        tmp.real = np.copy(self.grad_ref[0])
+        tmp.imag = np.copy(self.grad_ref[1])
+        self.grad_ref_parray = parray.to_device(self.la.queue, tmp)
+        self.grad_ref_buffer = self.grad_ref_parray.data
+
+    def tearDown(self):
+        self.image = None
+        self.image2 = None
+        self.grad = None
+        self.grad2 = None
+        self.grad_ref = None
+        self.div_ref = None
+        self.gradient_parray.data.release()
+        self.gradient_parray = None
+        self.gradient_buffer = None
+        self.image_parray.data.release()
+        self.image_parray = None
+        self.image_buffer = None
+        self.grad_ref_parray.data.release()
+        self.grad_ref_parray = None
+        self.grad_ref_buffer = None
+
+    def getfiles(self):
+        # load 512x512 MRI phantom - TODO include Lena or ascent once a .npz is available
+        self.image = np.load(utilstest.getfile("Brain512.npz"))["data"]
+
+    def compare(self, result, reference, abstol, name):
+        errmax = np.max(np.abs(result - reference))
+        logger.info("%s: Max error = %e" % (name, errmax))
+        self.assertTrue(errmax < abstol, str("%s: Max error is too high" % name))
+
+    @unittest.skipUnless(ocl and mako, "pyopencl is missing")
+    def test_gradient(self):
+        arrays = {
+            "numpy.ndarray": self.image,
+            "buffer": self.image_buffer,
+            "parray": self.image_parray
+        }
+        for desc, image in arrays.items():
+            # Test with dst on host (numpy.ndarray)
+            res = self.la.gradient(image, return_to_host=True)
+            self.compare(res, self.grad_ref, 1e-6, str("gradient[src=%s, dst=numpy.ndarray]" % desc))
+            # Test with dst on device (pyopencl.Buffer)
+            self.la.gradient(image, dst=self.gradient_buffer)
+            cl.enqueue_copy(self.la.queue, self.grad, self.gradient_buffer)
+            self.grad2[0] = self.grad.real
+            self.grad2[1] = self.grad.imag
+            self.compare(self.grad2, self.grad_ref, 1e-6, str("gradient[src=%s, dst=buffer]" % desc))
+            # Test with dst on device (pyopencl.Array)
+            self.la.gradient(image, dst=self.gradient_parray)
+            self.grad = self.gradient_parray.get()
+            self.grad2[0] = self.grad.real
+            self.grad2[1] = self.grad.imag
+            self.compare(self.grad2, self.grad_ref, 1e-6, str("gradient[src=%s, dst=parray]" % desc))
+
+    @unittest.skipUnless(ocl and mako, "pyopencl is missing")
+    def test_divergence(self):
+        arrays = {
+            "numpy.ndarray": self.grad_ref,
+            "buffer": self.grad_ref_buffer,
+            "parray": self.grad_ref_parray
+        }
+        for desc, grad in arrays.items():
+            # Test with dst on host (numpy.ndarray)
+            res = self.la.divergence(grad, return_to_host=True)
+            self.compare(res, self.div_ref, 1e-6, str("divergence[src=%s, dst=numpy.ndarray]" % desc))
+            # Test with dst on device (pyopencl.Buffer)
+            self.la.divergence(grad, dst=self.image_buffer)
+            cl.enqueue_copy(self.la.queue, self.image2, self.image_buffer)
+            self.compare(self.image2, self.div_ref, 1e-6, str("divergence[src=%s, dst=buffer]" % desc))
+            # Test with dst on device (pyopencl.Array)
+            self.la.divergence(grad, dst=self.image_parray)
+            self.image2 = self.image_parray.get()
+            self.compare(self.image2, self.div_ref, 1e-6, str("divergence[src=%s, dst=parray]" % desc))
+
+    @unittest.skipUnless(ocl and mako and _has_scipy, "pyopencl and/or scipy is missing")
+    def test_laplacian(self):
+        laplacian_ref = laplace(self.image)
+        # Laplacian = div(grad)
+        self.la.gradient(self.image)
+        laplacian_ocl = self.la.divergence(self.la.d_gradient, return_to_host=True)
+        self.compare(laplacian_ocl, laplacian_ref, 1e-6, "laplacian")
+
+
+def suite():
+    testSuite = unittest.TestSuite()
+    testSuite.addTest(TestLinAlg("test_gradient"))
+    testSuite.addTest(TestLinAlg("test_divergence"))
+    testSuite.addTest(TestLinAlg("test_laplacian"))
+    return testSuite
+
+
+if __name__ == '__main__':
+    unittest.main(defaultTest="suite")
diff --git a/silx/opencl/test/test_projection.py b/silx/opencl/test/test_projection.py
new file mode 100644
index 0000000..c9a3a1c
--- /dev/null
+++ b/silx/opencl/test/test_projection.py
@@ -0,0 +1,139 @@
+#!/usr/bin/env python
+# 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.
+#
+# ###########################################################################*/
+"""Test of the forward projection module"""
+
+from __future__ import division, print_function
+
+__authors__ = ["Pierre paleo"]
+__license__ = "MIT"
+__copyright__ = "2013-2017 European Synchrotron Radiation Facility, Grenoble, France"
+__date__ = "14/06/2017"
+
+
+import time
+import logging
+import numpy as np
+import unittest
+try:
+    import mako
+except ImportError:
+    mako = None
+from ..common import ocl
+if ocl:
+    from .. import projection
+from silx.test.utils import utilstest
+
+logger = logging.getLogger(__name__)
+
+
+
+@unittest.skipUnless(ocl and mako, "PyOpenCl is missing")
+class TestProj(unittest.TestCase):
+
+    def setUp(self):
+        if ocl is None:
+            return
+        #~ if sys.platform.startswith('darwin'):
+            #~ self.skipTest("Projection is not implemented on CPU for OS X yet")
+        self.getfiles()
+        n_angles = self.sino.shape[0]
+        self.proj = projection.Projection(self.phantom.shape, n_angles)
+        if self.proj.compiletime_workgroup_size < 16:
+            self.skipTest("Current implementation of OpenCL projection is not supported on this platform yet")
+
+
+    def tearDown(self):
+        self.phantom = None
+        self.sino = None
+        self.proj = None
+
+
+    def getfiles(self):
+        # load 512x512 MRI phantom
+        self.phantom = np.load(utilstest.getfile("Brain512.npz"))["data"]
+        # load sinogram computed with PyHST
+        self.sino = np.load(utilstest.getfile("sino500_pyhst.npz"))["data"]
+
+
+    def measure(self):
+        "Common measurement of timings"
+        t1 = time.time()
+        try:
+            result = self.proj.projection(self.phantom)
+        except RuntimeError as msg:
+            logger.error(msg)
+            return
+        t2 = time.time()
+        return t2 - t1, result
+
+
+    def compare(self, res):
+        """
+        Compare a result with the reference reconstruction.
+        Only the valid reconstruction zone (inscribed circle) is taken into account
+        """
+        # Compare with the original phantom.
+        # TODO: compare a standard projection
+        ref = self.sino
+        return np.max(np.abs(res - ref))
+
+
+    @unittest.skipUnless(ocl and mako, "pyopencl is missing")
+    def test_proj(self):
+        """
+        tests Projection
+        """
+        # Test single reconstruction
+        # --------------------------
+        t, res = self.measure()
+        if t is None:
+            logger.info("test_proj: skipped")
+        else:
+            logger.info("test_proj: time = %.3fs" % t)
+            err = self.compare(res)
+            msg = str("Max error = %e" % err)
+            logger.info(msg)
+            # Interpolation differs at some lines, giving relative error of 10/50000
+            self.assertTrue(err < 20., "Max error is too high")
+        # Test multiple reconstructions
+        # -----------------------------
+        res0 = np.copy(res)
+        for i in range(10):
+            res = self.proj.projection(self.phantom)
+            errmax = np.max(np.abs(res - res0))
+            self.assertTrue(errmax < 1.e-6, "Max error is too high")
+
+
+
+def suite():
+    testSuite = unittest.TestSuite()
+    testSuite.addTest(TestProj("test_proj"))
+    return testSuite
+
+
+
+if __name__ == '__main__':
+    unittest.main(defaultTest="suite")