From bfa4dba15485b4192f8bbe13345e9658c97ecf76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Picca=20Fr=C3=A9d=C3=A9ric-Emmanuel?= Date: Sat, 7 Oct 2017 07:59:01 +0200 Subject: New upstream version 0.6.0+dfsg --- silx/opencl/test/__init__.py | 21 +++- silx/opencl/test/test_array_utils.py | 161 ++++++++++++++++++++++++ silx/opencl/test/test_backprojection.py | 165 ++++++++++++++++++++++++ silx/opencl/test/test_linalg.py | 215 ++++++++++++++++++++++++++++++++ silx/opencl/test/test_projection.py | 139 +++++++++++++++++++++ 5 files changed, 697 insertions(+), 4 deletions(-) create mode 100644 silx/opencl/test/test_array_utils.py create mode 100644 silx/opencl/test/test_backprojection.py create mode 100644 silx/opencl/test/test_linalg.py create mode 100644 silx/opencl/test/test_projection.py (limited to 'silx/opencl/test') 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") -- cgit v1.2.3