summaryrefslogtreecommitdiff
path: root/silx/opencl/test/test_linalg.py
diff options
context:
space:
mode:
Diffstat (limited to 'silx/opencl/test/test_linalg.py')
-rw-r--r--silx/opencl/test/test_linalg.py215
1 files changed, 215 insertions, 0 deletions
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")