diff options
Diffstat (limited to 'silx/opencl/test/test_backprojection.py')
| -rw-r--r-- | silx/opencl/test/test_backprojection.py | 78 |
1 files changed, 62 insertions, 16 deletions
diff --git a/silx/opencl/test/test_backprojection.py b/silx/opencl/test/test_backprojection.py index 70ce2ae..7ab416d 100644 --- a/silx/opencl/test/test_backprojection.py +++ b/silx/opencl/test/test_backprojection.py @@ -35,8 +35,9 @@ __date__ = "19/01/2018" import time import logging -import numpy +import numpy as np import unittest +from math import pi try: import mako except ImportError: @@ -44,6 +45,7 @@ except ImportError: from ..common import ocl if ocl: from .. import backprojection + from ..sinofilter import compute_fourier_filter from silx.test.utils import utilstest logger = logging.getLogger(__name__) @@ -55,7 +57,7 @@ def generate_coords(img_shp, center=None): 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] + R, C = np.mgrid[:l_r, :l_c] if center is None: center0, center1 = l_r / 2., l_c / 2. else: @@ -71,7 +73,7 @@ def clip_circle(img, center=None, radius=None): """ R, C = generate_coords(img.shape, center) M = R * R + C * C - res = numpy.zeros_like(img) + res = np.zeros_like(img) if radius is None: radius = img.shape[0] / 2. - 1 mask = M < radius * radius @@ -85,23 +87,29 @@ 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 * 16: - self.skipTest("Current implementation of OpenCL backprojection is not supported on this platform yet") + self.skipTest("Current implementation of OpenCL backprojection is " + "not supported on this platform yet") + # Astra does not use the same backprojector implementation. + # Therefore, we cannot expect results to be the "same" (up to float32 + # numerical error) + self.tol = 5e-2 + if self.fbp.is_cpu: + # Precision is less when using CPU + self.tol *= 2 def tearDown(self): self.sino = None -# self.fbp.log_profile() + # 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"] + self.sino = np.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"] + self.reference_rec = np.load(utilstest.getfile("rec_astra_500.npz"))["data"] def measure(self): "Common measurement of timings" @@ -124,8 +132,8 @@ class TestFBP(unittest.TestCase): 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)) + logger.debug("Absolute difference: %s with %s outlier pixels out of %s" + "", delta.max(), bad.sum(), np.prod(bad.shape)) return delta.max() @unittest.skipUnless(ocl and mako, "pyopencl is missing") @@ -143,21 +151,59 @@ class TestFBP(unittest.TestCase): 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") + self.assertTrue(err < self.tol, "Max error is too high") + # Test multiple reconstructions # ----------------------------- - res0 = numpy.copy(res) + res0 = np.copy(res) for i in range(10): res = self.fbp.filtered_backprojection(self.sino) - errmax = numpy.max(numpy.abs(res - res0)) + errmax = np.max(np.abs(res - res0)) self.assertTrue(errmax < 1.e-6, "Max error is too high") + @unittest.skipUnless(ocl and mako, "pyopencl is missing") + def test_fbp_filters(self): + """ + Test the different available filters of silx FBP. + """ + avail_filters = [ + "ramlak", "shepp-logan", "cosine", "hamming", + "hann" + ] + # Create a Dirac delta function at a single angle view. + # As the filters are radially invarant: + # - backprojection yields an image where each line is a Dirac. + # - FBP yields an image where each line is the spatial filter + # One can simply filter "dirac" without backprojecting it, but this + # test will also ensure that backprojection behaves well. + dirac = np.zeros_like(self.sino) + na, dw = dirac.shape + dirac[0, dw//2] = na / pi * 2 + + for filter_name in avail_filters: + B = backprojection.Backprojection(dirac.shape, filter_name=filter_name) + r = B(dirac) + # Check that radial invariance is kept + std0 = np.max(np.abs(np.std(r, axis=0))) + self.assertTrue( + std0 < 5.e-6, + "Something wrong with FBP(filter=%s)" % filter_name + ) + # Check that the filter is retrieved + r_f = np.fft.fft(np.fft.fftshift(r[0])).real / 2. # filter factor + ref_filter_f = compute_fourier_filter(dw, filter_name) + errmax = np.max(np.abs(r_f - ref_filter_f)) + logger.info("FBP filter %s: max error=%e" % (filter_name, errmax)) + self.assertTrue( + errmax < 1.e-3, + "Something wrong with FBP(filter=%s)" % filter_name + ) + def suite(): testSuite = unittest.TestSuite() testSuite.addTest(TestFBP("test_fbp")) + testSuite.addTest(TestFBP("test_fbp_filters")) return testSuite |