[Tomo] FBP test fails under windows
investigation under way
diff --git a/silx/opencl/backprojection.py b/silx/opencl/backprojection.py
index 92ba89d47..15a03b936 100644
--- a/silx/opencl/backprojection.py
+++ b/silx/opencl/backprojection.py
@@ -29,10 +29,10 @@ from __future__ import absolute_import, print_function, with_statement, division
 
 __authors__ = ["A. Mirone, P. Paleo"]
 __license__ = "MIT"
-__date__ = "12/09/2017"
+__date__ = "05/10/2017"
 
 import logging
-import numpy as np
+import numpy
 
 from .common import pyopencl
 from .processing import EventDescription, OpenclProcessing, BufferDescription
@@ -51,13 +51,15 @@ try:
     _has_pyfft = True
 except ImportError:
     _has_pyfft = False
+# For silx v0.6 we disable the use FFT on GPU.
+_has_pyfft = False
 
 
 def _sizeof(Type):
     """
     return the size (in bytes) of a scalar type, like the C behavior
     """
-    return np.dtype(Type).itemsize
+    return numpy.dtype(Type).itemsize
 
 
 def _idivup(a, b):
@@ -67,6 +69,44 @@ def _idivup(a, b):
     return (a + (b - 1)) // b
 
 
+def fourier_filter(sino, filter_=None, fft_size=None):
+    """Simple numpy based implementation of fourier space filter
+    
+    :param sino: of shape shape = (num_projs, num_bins)
+    :param filter: filter function to apply in fourier space
+    :fft_size: size on which perform the fft. May be larger than the sino array 
+    :return: filtered sinogram
+    """
+    assert sino.ndim == 2
+    num_projs, num_bins = sino.shape
+    if fft_size is None:
+        fft_size = nextpow2(num_bins * 2 - 1)
+    else:
+        assert fft_size >= num_bins
+    if fft_size == num_bins:
+        sino_zeropadded = sino.astype(numpy.float32)
+    else:
+        sino_zeropadded = numpy.zeros((num_projs, fft_size),
+                                      dtype=numpy.complex64)
+        sino_zeropadded[:, :num_bins] = sino.astype(numpy.float32)
+
+    if filter_ is None:
+        h = numpy.zeros(fft_size, dtype=numpy.float32)
+        L2 = fft_size // 2 + 1
+        h[0] = 1 / 4.
+        j = numpy.linspace(1, L2, L2 // 2, False)
+        h[1:L2:2] = -1. / (numpy.pi ** 2 * j ** 2)
+        h[L2:] = numpy.copy(h[1:L2 - 1][::-1])
+        filter_ = numpy.fft.fft(h).astype(numpy.complex64)
+
+    # Linear convolution
+    sino_f = numpy.fft.fft(sino, fft_size)
+    sino_f = sino_f * filter_
+    sino_filtered = numpy.fft.ifft(sino_f)[:, :num_bins].real
+    # Send the filtered sinogram to device
+    return numpy.ascontiguousarray(sino_filtered.real, dtype=numpy.float32)
+
+
 class Backprojection(OpenclProcessing):
     """A class for performing the backprojection using OpenCL"""
     kernel_files = ["backproj.cl", "array_utils.cl"]
@@ -110,8 +150,8 @@ class Backprojection(OpenclProcessing):
                                   profile=profile)
         self.shape = sino_shape
 
-        self.num_bins = np.int32(sino_shape[1])
-        self.num_projs = np.int32(sino_shape[0])
+        self.num_bins = numpy.int32(sino_shape[1])
+        self.num_projs = numpy.int32(sino_shape[0])
         self.angles = angles
         if slice_shape is None:
             self.slice_shape = (self.num_bins, self.num_bins)
@@ -121,24 +161,25 @@ class Backprojection(OpenclProcessing):
             _idivup(self.slice_shape[0], 32) * 32,
             _idivup(self.slice_shape[1], 32) * 32
         )
-        self.slice = np.zeros(self.dimrec_shape, dtype=np.float32)
+        self.slice = numpy.zeros(self.dimrec_shape, dtype=numpy.float32)
         self.filter_name = filter_name if filter_name else "Ram-Lak"
         if axis_position:
-            self.axis_pos = np.float32(axis_position)
+            self.axis_pos = numpy.float32(axis_position)
         else:
-            self.axis_pos = np.float32((sino_shape[1] - 1.) / 2)
+            self.axis_pos = numpy.float32((sino_shape[1] - 1.) / 2)
         self.axis_array = None  # TODO: add axis correction front-end
+
         self.is_cpu = False
         if self.device.type == "CPU":
             self.is_cpu = True
 
         self.compute_fft_plans()
         self.buffers = [
-                       BufferDescription("_d_slice", np.prod(self.dimrec_shape), np.float32, mf.READ_WRITE),
-                       BufferDescription("d_sino", self.num_projs * self.num_bins, np.float32, mf.READ_WRITE),  # before transferring to texture (if available)
-                       BufferDescription("d_cos", self.num_projs, np.float32, mf.READ_ONLY),
-                       BufferDescription("d_sin", self.num_projs, np.float32, mf.READ_ONLY),
-                       BufferDescription("d_axes", self.num_projs, np.float32, mf.READ_ONLY),
+                       BufferDescription("_d_slice", numpy.prod(self.dimrec_shape), numpy.float32, mf.READ_WRITE),
+                       BufferDescription("d_sino", self.num_projs * self.num_bins, numpy.float32, mf.READ_WRITE),  # before transferring to texture (if available)
+                       BufferDescription("d_cos", self.num_projs, numpy.float32, mf.READ_ONLY),
+                       BufferDescription("d_sin", self.num_projs, numpy.float32, mf.READ_ONLY),
+                       BufferDescription("d_axes", self.num_projs, numpy.float32, mf.READ_ONLY),
                       ]
         self.allocate_buffers()
         if not(self.is_cpu):
@@ -149,10 +190,10 @@ class Backprojection(OpenclProcessing):
                 "d_filter": self.d_filter,
                 "d_sino_z": self.d_sino_z
             })
-        self.d_sino = self.cl_mem["d_sino"] # shorthand
+        self.d_sino = self.cl_mem["d_sino"]  # shorthand
         self.compute_angles()
 
-        self.local_mem = 256 * 3 * _sizeof(np.float32)  # constant for all image sizes
+        self.local_mem = 256 * 3 * _sizeof(numpy.float32)  # constant for all image sizes
         OpenclProcessing.compile_kernels(self, self.kernel_files)
         # check that workgroup can actually be (16, 16)
         self.check_workgroup_size("backproj_cpu_kernel")
@@ -165,19 +206,19 @@ class Backprojection(OpenclProcessing):
 
     def compute_angles(self):
         if self.angles is None:
-            self.angles = np.linspace(0, np.pi, self.num_projs, False)
-        h_cos = np.cos(self.angles).astype(np.float32)
-        h_sin = np.sin(self.angles).astype(np.float32)
+            self.angles = numpy.linspace(0, numpy.pi, self.num_projs, False)
+        h_cos = numpy.cos(self.angles).astype(numpy.float32)
+        h_sin = numpy.sin(self.angles).astype(numpy.float32)
         pyopencl.enqueue_copy(self.queue, self.cl_mem["d_cos"], h_cos)
         pyopencl.enqueue_copy(self.queue, self.cl_mem["d_sin"], h_sin)
         if self.axis_array:
             pyopencl.enqueue_copy(self.queue,
                                   self.cl_mem["d_axes"],
-                                  self.axis_array.astype(np.float32))
+                                  self.axis_array.astype(numpy.float32))
         else:
             pyopencl.enqueue_copy(self.queue,
                                   self.cl_mem["d_axes"],
-                                  np.ones(self.num_projs, dtype=np.float32) * self.axis_pos)
+                                  numpy.ones(self.num_projs, dtype=numpy.float32) * self.axis_pos)
 
     def allocate_textures(self):
         """
@@ -190,7 +231,7 @@ class Backprojection(OpenclProcessing):
                                                              pyopencl.channel_order.INTENSITY,
                                                              pyopencl.channel_type.FLOAT
                                                             ),
-                                        hostbuf=np.zeros(self.shape[::-1], dtype=np.float32)
+                                        hostbuf=numpy.zeros(self.shape[::-1], dtype=numpy.float32)
                                         )
 
     def compute_fft_plans(self):
@@ -203,10 +244,11 @@ class Backprojection(OpenclProcessing):
         if _has_pyfft:
             logger.debug("pyfft is available. Computing FFT plans...")
             # batched 1D transform
-            self.pyfft_plan = pyfft_Plan(self.fft_size, queue=self.queue)
+            self.pyfft_plan = pyfft_Plan(self.fft_size, queue=self.queue,
+                                         wait_for_finish=True)
             self.d_sino_z = parray.zeros(self.queue,
                                          (self.num_projs, self.fft_size),
-                                         dtype=np.complex64)
+                                         dtype=numpy.complex64)
             logger.debug("... done")
         else:
             logger.debug("pyfft not available, using numpy.fft")
@@ -219,80 +261,84 @@ class Backprojection(OpenclProcessing):
         """
         if self.filter_name == "Ram-Lak":
             L = self.fft_size
-            h = np.zeros(L, dtype=np.float32)
+            h = numpy.zeros(L, dtype=numpy.float32)
             L2 = L // 2 + 1
             h[0] = 1 / 4.
-            j = np.linspace(1, L2, L2 // 2, False)
-            h[1:L2:2] = -1. / (np.pi ** 2 * j ** 2)
-            h[L2:] = np.copy(h[1:L2 - 1][::-1])
+            j = numpy.linspace(1, L2, L2 // 2, False)
+            h[1:L2:2] = -1. / (numpy.pi ** 2 * j ** 2)
+            h[L2:] = numpy.copy(h[1:L2 - 1][::-1])
         else:
             # TODO: other filters
             raise ValueError("Filter %s is not available" % self.filter_name)
         self.filter = h
         if self.pyfft_plan:
-            self.d_filter = parray.to_device(self.queue, h.astype(np.complex64))
+            self.d_filter = parray.to_device(self.queue, h.astype(numpy.complex64))
             self.pyfft_plan.execute(self.d_filter.data)
         else:
-            self.filter = np.fft.fft(h).astype(np.complex64)
+            self.filter = numpy.fft.fft(h).astype(numpy.complex64)
             self.d_filter = None
 
     def _get_local_mem(self):
         return pyopencl.LocalMemory(self.local_mem)  # constant for all image sizes
 
     def cpy2d_to_slice(self, dst):
-        ndrange = (int(self.slice_shape[1]), int(self.slice_shape[0])) # pyopencl < 2015.2
-        slice_shape_ocl = np.int32(ndrange)
+        ndrange = (int(self.slice_shape[1]), int(self.slice_shape[0]))  # pyopencl < 2015.2
+        slice_shape_ocl = numpy.int32(ndrange)
         wg = None
         kernel_args = (
             dst.data,
             self.cl_mem["_d_slice"],
-            np.int32(self.slice_shape[1]),
-            np.int32(self.dimrec_shape[1]),
-            np.int32((0, 0)),
-            np.int32((0, 0)),
+            numpy.int32(self.slice_shape[1]),
+            numpy.int32(self.dimrec_shape[1]),
+            numpy.int32((0, 0)),
+            numpy.int32((0, 0)),
             slice_shape_ocl
         )
         return self.kernels.cpy2d(self.queue, ndrange, wg, *kernel_args)
 
-
     def transfer_to_texture(self, sino):
         sino2 = sino
-        if not(sino.flags["C_CONTIGUOUS"] and sino.dtype == np.float32):
-            sino2 = np.ascontiguousarray(sino, dtype=np.float32)
+        if not(sino.flags["C_CONTIGUOUS"] and sino.dtype == numpy.float32):
+            sino2 = numpy.ascontiguousarray(sino, dtype=numpy.float32)
         if self.is_cpu:
-            return pyopencl.enqueue_copy(
-                self.queue,
-                self.d_sino,
-                sino2
-            )
+            ev = pyopencl.enqueue_copy(
+                                        self.queue,
+                                        self.d_sino,
+                                        sino2
+                                        )
+            what = "transfer filtered sino H->D buffer"
         else:
-            return pyopencl.enqueue_copy(
-                self.queue,
-                self.d_sino_tex,
-                sino2,
-                origin=(0, 0),
-                region=self.shape[::-1]
-            )
+            ev = pyopencl.enqueue_copy(
+                                       self.queue,
+                                       self.d_sino_tex,
+                                       sino2,
+                                       origin=(0, 0),
+                                       region=self.shape[::-1]
+                                       )
+            what = "transfer filtered sino H->D texture"
+        return EventDescription(what, ev)
 
     def transfer_device_to_texture(self, d_sino):
         if self.is_cpu:
             if id(self.d_sino) == id(d_sino):
                 return
-            return pyopencl.enqueue_copy(
-                self.queue,
-                self.d_sino,
-                d_sino
-            )
+            ev = pyopencl.enqueue_copy(
+                                       self.queue,
+                                       self.d_sino,
+                                       d_sino
+                                       )
+            what = "transfer filtered sino D->D buffer"
         else:
-            return pyopencl.enqueue_copy(
-                self.queue,
-                self.d_sino_tex,
-                d_sino,
-                offset=0,
-                origin=(0, 0),
-                region=self.shape[::-1]
-            )
-
+            ev = pyopencl.enqueue_copy(
+                                       self.queue,
+                                       self.d_sino_tex,
+                                       d_sino,
+                                       offset=0,
+                                       origin=(0, 0),
+                                       region=self.shape[::-1]
+                                       )
+            what = "transfer filtered sino D->D texture"
+        return EventDescription(what, ev)
 
     def backprojection(self, sino=None, dst=None):
         """Perform the backprojection on an input sinogram
@@ -304,8 +350,8 @@ class Backprojection(OpenclProcessing):
         events = []
         with self.sem:
 
-            if sino is not None: # assuming numpy.ndarray
-                self.transfer_to_texture(sino)
+            if sino is not None:  # assuming numpy.ndarray
+                events.append(self.transfer_to_texture(sino))
             # Prepare arguments for the kernel call
             if self.is_cpu:
                 d_sino_ref = self.d_sino
@@ -317,8 +363,8 @@ class Backprojection(OpenclProcessing):
                 self.axis_pos,  # axis position (float32)
                 self.cl_mem["_d_slice"],  # d_slice (__global float32*)
                 d_sino_ref,  # d_sino (__read_only image2d_t or float*)
-                np.float32(0),  # gpu_offset_x (float32)
-                np.float32(0),  # gpu_offset_y (float32)
+                numpy.float32(0),  # gpu_offset_x (float32)
+                numpy.float32(0),  # gpu_offset_y (float32)
                 self.cl_mem["d_cos"],  # d_cos (__global float32*)
                 self.cl_mem["d_sin"],  # d_sin (__global float32*)
                 self.cl_mem["d_axes"],  # d_axis  (__global float32*)
@@ -342,7 +388,7 @@ class Backprojection(OpenclProcessing):
                                            self.cl_mem["_d_slice"])
                 events.append(EventDescription("copy D->H result", ev))
                 ev.wait()
-                res = np.copy(self.slice)
+                res = numpy.copy(self.slice)
                 if self.dimrec_shape[0] > self.slice_shape[0] or self.dimrec_shape[1] > self.slice_shape[1]:
                     res = res[:self.slice_shape[0], :self.slice_shape[1]]
                 # if the slice is backprojected onto a bigger grid
@@ -371,30 +417,33 @@ class Backprojection(OpenclProcessing):
         if sino.shape[0] != self.num_projs or sino.shape[1] != self.num_bins:
             raise ValueError("Expected sinogram with (projs, bins) = (%d, %d)" % (self.num_projs, self.num_bins))
         if rescale:
-            sino = sino * np.pi / self.num_projs
+            sino = sino * numpy.pi / self.num_projs
+        events = []
         # if pyfft is available, all can be done on the device
         if self.d_filter is not None:
-            events = []
+
             # Zero-pad the sinogram.
             # TODO: this can be done on GPU with a "Memcpy2D":
             #  cl.enqueue_copy(queue, dst, src, host_origin=(0,0), buffer_origin=(0,0), region=shape, host_pitches=(sino.shape[1],), buffer_pitches=(self.fft_size,))
             # However it does not work properly, and raises an error for pyopencl < 2017.1
-            sino_zeropadded = np.zeros((sino.shape[0], self.fft_size), dtype=np.complex64)
-            sino_zeropadded[:, :self.num_bins] = sino.astype(np.float32)
-            sino_zeropadded = np.ascontiguousarray(sino_zeropadded, dtype=np.complex64)
+            sino_zeropadded = numpy.zeros((sino.shape[0], self.fft_size), dtype=numpy.complex64)
+            sino_zeropadded[:, :self.num_bins] = sino.astype(numpy.float32)
+            sino_zeropadded = numpy.ascontiguousarray(sino_zeropadded, dtype=numpy.complex64)
             with self.sem:
                 # send to GPU
                 ev = pyopencl.enqueue_copy(self.queue, self.d_sino_z.data, sino_zeropadded)
                 events.append(EventDescription("Send sino H->D", ev))
+
                 # FFT (in-place)
                 self.pyfft_plan.execute(self.d_sino_z.data, batch=self.num_projs)
+
                 # Multiply (complex-wise) with the the filter
                 ev = self.kernels.mult(self.queue,
                                        tuple(int(i) for i in self.d_sino_z.shape[::-1]),
                                        None,
                                        self.d_sino_z.data,
                                        self.d_filter.data,
-                                       np.int32(self.fft_size),
+                                       numpy.int32(self.fft_size),
                                        self.num_projs
                                        )
                 events.append(EventDescription("complex 2D-1D multiplication", ev))
@@ -406,28 +455,23 @@ class Backprojection(OpenclProcessing):
                                             self.d_sino_z.data,
                                             self.num_bins,
                                             self.num_projs,
-                                            np.int32(self.fft_size)
+                                            numpy.int32(self.fft_size)
                                             )
                 events.append(EventDescription("conversion from complex padded sinogram to sinogram", ev))
+                # debug
+#                 ev.wait()
+#                 h_sino = numpy.zeros(sino.shape, dtype=numpy.float32)
+#                 ev = pyopencl.enqueue_copy(self.queue, h_sino, self.d_sino)
+#                 ev.wait()
+#                 numpy.save("/tmp/filtered_sinogram_%s.npy" % self.ctx.devices[0].platform.name.split()[0], h_sino)
                 events.append(self.transfer_device_to_texture(self.d_sino))
-            if self.profile:
-                self.events += events
             # ------
         else:  # no pyfft
-            # Zero-padding of the sinogram
-            sino_zeropadded = np.zeros((sino.shape[0],
-                                        self.fft_size),
-                                       dtype=np.complex64)
-            sino_zeropadded[:, :self.num_bins] = sino.astype(np.float32)
-            # Linear convolution
-            sino_f = np.fft.fft(sino, self.fft_size)
-            sino_f = sino_f * self.filter
-            sino_filtered = np.fft.ifft(sino_f)[:, :self.num_bins].real
-            # Send the filtered sinogram to device
-            sino_filtered = np.ascontiguousarray(sino_filtered,
-                                                 dtype=np.float32)
+            sino_filtered = fourier_filter(sino, filter_=self.filter, fft_size=self.fft_size)
             with self.sem:
-                pyopencl.enqueue_copy(self.queue, self.d_sino, sino_filtered)
+                events.append(self.transfer_to_texture(sino_filtered))
+        if self.profile:
+            self.events += events
 
     def filtered_backprojection(self, sino):
         """
diff --git a/silx/opencl/common.py b/silx/opencl/common.py
index 0e94f1af5..ebf50c77b 100644
--- a/silx/opencl/common.py
+++ b/silx/opencl/common.py
@@ -34,7 +34,7 @@ __author__ = "Jerome Kieffer"
 __contact__ = "Jerome.Kieffer@ESRF.eu"
 __license__ = "MIT"
 __copyright__ = "2012-2017 European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "02/08/2017"
+__date__ = "05/10/2017"
 __status__ = "stable"
 __all__ = ["ocl", "pyopencl", "mf", "release_cl_buffers", "allocate_cl_buffers",
            "measure_workgroup_size", "kernel_workgroup_size"]
@@ -419,6 +419,10 @@ class OpenCL(object):
         if (platformid is not None) and (deviceid is not None):
             platformid = int(platformid)
             deviceid = int(deviceid)
+        elif "PYOPENCL_CTX" in os.environ:
+            pyopencl_ctx = [int(i) if i.isdigit() else 0 for i in os.environ["PYOPENCL_CTX"].split(":")]
+            pyopencl_ctx += [0] * (2 - len(pyopencl_ctx))  # pad with 0
+            platformid, deviceid = pyopencl_ctx
         else:
             if useFp64:
                 ids = ocl.select_device(type=devicetype, extensions=["cl_khr_int64_base_atomics"])
diff --git a/silx/opencl/processing.py b/silx/opencl/processing.py
index b7eda265a..1997a5518 100644
--- a/silx/opencl/processing.py
+++ b/silx/opencl/processing.py
@@ -41,7 +41,7 @@ __author__ = "Jerome Kieffer"
 __contact__ = "Jerome.Kieffer@ESRF.eu"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "12/09/2017"
+__date__ = "03/10/2017"
 __status__ = "stable"
 
 
@@ -167,8 +167,8 @@ class OpenclProcessing(object):
             ualloc = 0
             for buf in buffers:
                 ualloc += numpy.dtype(buf.dtype).itemsize * buf.size
-            logger.info("%.3fMB are needed on device which has %.3fMB",
-                        ualloc / 1.0e6, self.device.memory / 1.0e6)
+            logger.info("%.3fMB are needed on device: %s,  which has %.3fMB",
+                        ualloc / 1.0e6, self.device, self.device.memory / 1.0e6)
 
             if ualloc >= self.device.memory:
                 raise MemoryError("Fatal error in allocate_buffers. Not enough "
diff --git a/silx/opencl/test/test_backprojection.py b/silx/opencl/test/test_backprojection.py
index cf0ef69b3..342bd2f8e 100644
--- a/silx/opencl/test/test_backprojection.py
+++ b/silx/opencl/test/test_backprojection.py
@@ -30,12 +30,12 @@ from __future__ import division, print_function
 __authors__ = ["Pierre paleo"]
 __license__ = "MIT"
 __copyright__ = "2013-2017 European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "14/06/2017"
+__date__ = "05/10/2017"
 
 
 import time
 import logging
-import numpy as np
+import numpy
 import unittest
 try:
     import mako
@@ -55,7 +55,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 = np.mgrid[:l_r, :l_c]
+    R, C = numpy.mgrid[:l_r, :l_c]
     if center is None:
         center0, center1 = l_r / 2., l_c / 2.
     else:
@@ -70,11 +70,12 @@ 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**2+C**2
-    res = np.zeros_like(img)
+    M = R * R + C * C
+    res = numpy.zeros_like(img)
     if radius is None:
-        radius = img.shape[0]/2.-1
-    res[M < radius**2] = img[M < radius**2]
+        radius = img.shape[0] / 2. - 1
+    mask = M < radius * radius
+    res[mask] = img[mask]
     return res
 
 
@@ -84,22 +85,23 @@ 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")
+        # ~ 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)
+        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 = np.load(utilstest.getfile("sino500.npz"))["data"]
+        self.sino = numpy.load(utilstest.getfile("sino500.npz"))["data"]
         # load reconstruction made with ASTRA FBP (with filter designed in spatial domain)
-        self.reference_rec = np.load(utilstest.getfile("rec_astra_500.npz"))["data"]
+        self.reference_rec = numpy.load(utilstest.getfile("rec_astra_500.npz"))["data"]
 
     def measure(self):
         "Common measurement of timings"
@@ -120,7 +122,11 @@ class TestFBP(unittest.TestCase):
         """
         res_clipped = clip_circle(res)
         ref_clipped = clip_circle(self.reference_rec)
-        return np.max(np.abs(res_clipped - ref_clipped))
+        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):
@@ -142,10 +148,10 @@ class TestFBP(unittest.TestCase):
             self.assertTrue(err < 1., "Max error is too high")
         # Test multiple reconstructions
         # -----------------------------
-        res0 = np.copy(res)
+        res0 = numpy.copy(res)
         for i in range(10):
             res = self.fbp.filtered_backprojection(self.sino)
-            errmax = np.max(np.abs(res - res0))
+            errmax = numpy.max(numpy.abs(res - res0))
             self.assertTrue(errmax < 1.e-6, "Max error is too high")
 
 
diff --git a/silx/resources/opencl/addition.cl b/silx/resources/opencl/addition.cl
index 8ecfd4e8f..35d799649 100644
--- a/silx/resources/opencl/addition.cl
+++ b/silx/resources/opencl/addition.cl
@@ -27,10 +27,16 @@
  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  * OTHER DEALINGS IN THE SOFTWARE.
  */
-__kernel void addition(__global float* a, __global float* b, __global float* res, int N)
+
+
+// "Hello_world" kernel to test if OpenCL is actually working
+kernel void addition(global float* a,
+                     global float* b,
+                     global float* res,
+                     int N)
 {
-    unsigned int i = get_global_id(0);
+    int i = get_global_id(0);
     if( i<N ){
         res[i] = a[i] + b[i];
     }
-}
\ No newline at end of file
+}
diff --git a/silx/resources/opencl/array_utils.cl b/silx/resources/opencl/array_utils.cl
index 7e279d7b1..60677dc5e 100644
--- a/silx/resources/opencl/array_utils.cl
+++ b/silx/resources/opencl/array_utils.cl
@@ -14,18 +14,19 @@
  * :param transfer_shape: shape of the transfer array in the form (x, y)
  *
  */
-__kernel void cpy2d(
-                    __global float* dst,
-                    __global float* src,
+kernel void cpy2d(
+                    global float* dst,
+                    global float* src,
                     int dst_width,
                     int src_width,
                     int2 dst_offset,
                     int2 src_offset,
                     int2 transfer_shape)
 {
-    uint gidx = get_global_id(0);
-    uint gidy = get_global_id(1);
-    if (gidx < transfer_shape.x && gidy < transfer_shape.y) {
+    int gidx = get_global_id(0),
+        gidy = get_global_id(1);
+    if (gidx < transfer_shape.x && gidy < transfer_shape.y)
+    {
         dst[(dst_offset.y + gidy)*dst_width + (dst_offset.x + gidx)] = src[(src_offset.y + gidy)*src_width + (src_offset.x + gidx)];
     }
 }
diff --git a/silx/resources/opencl/backproj.cl b/silx/resources/opencl/backproj.cl
index 162d94e59..6fadc2c3d 100644
--- a/silx/resources/opencl/backproj.cl
+++ b/silx/resources/opencl/backproj.cl
@@ -36,37 +36,33 @@
 /*******************************************************************************/
 
 
-
-
-
-const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;
-
-__kernel void backproj_kernel(
+kernel void backproj_kernel(
     int num_proj,
     int num_bins,
     float axis_position,
-    __global float *d_SLICE,
-    __read_only image2d_t d_sino,
+    global float *d_SLICE,
+    read_only image2d_t d_sino,
     float gpu_offset_x,
     float gpu_offset_y,
-    __global float * d_cos_s, // precalculated cos(theta[i])
-    __global float * d_sin_s, // precalculated sin(theta[i])
-    __global float * d_axis_s, // array of axis positions (n_projs)
-    __local float* shared2)     // 768B of local mem
+    global float * d_cos_s,  // precalculated cos(theta[i])
+    global float* d_sin_s,   // precalculated sin(theta[i])
+    global float* d_axis_s,  // array of axis positions (n_projs)
+    local float* shared2)    // 768B of local mem
 {
+    const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;
     const int tidx = get_local_id(0); //threadIdx.x;
     const int bidx = get_group_id(0); //blockIdx.x;
     const int tidy = get_local_id(1); //threadIdx.y;
     const int bidy = get_group_id(1); //blockIdx.y;
 
-    //~ __local float shared[768];
+    //~ local float shared[768];
     //~ float  * sh_sin  = shared;
     //~ float  * sh_cos  = shared+256;
     //~ float  * sh_axis = sh_cos+256;
 
-    __local float sh_cos[256];
-    __local float sh_sin[256];
-    __local float sh_axis[256];
+    local float sh_cos[256];
+    local float sh_sin[256];
+    local float sh_axis[256];
 
     float pcos, psin;
     float h0, h1, h2, h3;
@@ -102,9 +98,9 @@ __kernel void backproj_kernel(
         h3 = (acorr05 + (bx00+1)*pcos - (by00+1)*psin);
 
         if(h0>=0 && h0<num_bins) res0 += read_imagef(d_sino, sampler, (float2) (h0 +0.5f,proj +0.5f)).x; // tex2D(texprojs,h0 +0.5f,proj +0.5f);
-        if(h1>=0 && h1<num_bins) res1 += read_imagef(d_sino, sampler, (float2) (h1 +0.5f,proj +0.5f)).x; //tex2D(texprojs,h1 +0.5f,proj +0.5f);
-        if(h2>=0 && h2<num_bins) res2 += read_imagef(d_sino, sampler, (float2) (h2 +0.5f,proj +0.5f)).x; //tex2D(texprojs,h2 +0.5f,proj +0.5f);
-        if(h3>=0 && h3<num_bins) res3 += read_imagef(d_sino, sampler, (float2) (h3 +0.5f,proj +0.5f)).x; //tex2D(texprojs,h3 +0.5f,proj +0.5f);
+        if(h1>=0 && h1<num_bins) res1 += read_imagef(d_sino, sampler, (float2) (h1 +0.5f,proj +0.5f)).x; // tex2D(texprojs,h1 +0.5f,proj +0.5f);
+        if(h2>=0 && h2<num_bins) res2 += read_imagef(d_sino, sampler, (float2) (h2 +0.5f,proj +0.5f)).x; // tex2D(texprojs,h2 +0.5f,proj +0.5f);
+        if(h3>=0 && h3<num_bins) res3 += read_imagef(d_sino, sampler, (float2) (h3 +0.5f,proj +0.5f)).x; // tex2D(texprojs,h3 +0.5f,proj +0.5f);
     }
     d_SLICE[ 32*get_num_groups(0)*(bidy*32+tidy*2+0) + bidx*32 + tidx*2 + 0] = res0;
     d_SLICE[ 32*get_num_groups(0)*(bidy*32+tidy*2+1) + bidx*32 + tidx*2 + 0] = res1;
@@ -130,11 +126,16 @@ __kernel void backproj_kernel(
 
 #define ADJACENT_PIXELS_VALS(arr, Nx, y, xm, xp) ((float2) (arr[y*Nx+xm], arr[y*Nx+xp]))
 
-float linear_interpolation(float2 vals, float x, int xm, int xp) {
+//Simple linear interpolator for working on the GPU
+static float linear_interpolation(float2 vals,
+                                  float x,
+                                  int xm,
+                                  int xp)
+{
     if (xm == xp)
-	return vals.s0; 
+        return vals.s0;
     else 
-	return (vals.s0 * (xp - x)) + (vals.s1 * (x - xm)); 
+        return (vals.s0 * (xp - x)) + (vals.s1 * (x - xm));
 }
 
 /**
@@ -142,27 +143,27 @@ float linear_interpolation(float2 vals, float x, int xm, int xp) {
  *  Same kernel as backproj_kernel, but targets the CPU (no texture)
  *
 **/
-__kernel void backproj_cpu_kernel(
+kernel void backproj_cpu_kernel(
     int num_proj,
     int num_bins,
     float axis_position,
-    __global float *d_SLICE,
-    __global float* d_sino,
+    global float *d_SLICE,
+    global float* d_sino,
     float gpu_offset_x,
     float gpu_offset_y,
-    __global float * d_cos_s, // precalculated cos(theta[i])
-    __global float * d_sin_s, // precalculated sin(theta[i])
-    __global float * d_axis_s, // array of axis positions (n_projs)
-    __local float* shared2)     // 768B of local mem
+    global float * d_cos_s, // precalculated cos(theta[i])
+    global float * d_sin_s, // precalculated sin(theta[i])
+    global float * d_axis_s, // array of axis positions (n_projs)
+    local float* shared2)     // 768B of local mem
 {
     const int tidx = get_local_id(0); //threadIdx.x;
     const int bidx = get_group_id(0); //blockIdx.x;
     const int tidy = get_local_id(1); //threadIdx.y;
     const int bidy = get_group_id(1); //blockIdx.y;
 
-    __local float sh_cos[256];
-    __local float sh_sin[256];
-    __local float sh_axis[256];
+    local float sh_cos[256];
+    local float sh_sin[256];
+    local float sh_axis[256];
 
     float pcos, psin;
     float h0, h1, h2, h3;
diff --git a/silx/resources/opencl/backproj_helper.cl b/silx/resources/opencl/backproj_helper.cl
index f0232c379..b1590f83d 100644
--- a/silx/resources/opencl/backproj_helper.cl
+++ b/silx/resources/opencl/backproj_helper.cl
@@ -30,27 +30,28 @@
  */
 
 
-__kernel void mult(
-                    __global float2* d_sino,
-                    __global float2* d_filter,
+kernel void mult(   global float2* d_sino,
+                    global float2* d_filter,
                     int num_bins,
                     int num_projs)
+{
+    int gid0 = get_global_id(0);
+    int gid1 = get_global_id(1);
+    if (gid0 < num_bins && gid1 < num_projs)
     {
-      int gid0 = get_global_id(0);
-      int gid1 = get_global_id(1);
-      if (gid0 < num_bins && gid1 < num_projs) {
         // d_sino[gid1*num_bins+gid0] *= d_filter[gid0];
         d_sino[gid1*num_bins+gid0].x *= d_filter[gid0].x;
         d_sino[gid1*num_bins+gid0].y *= d_filter[gid0].x;
       }
-    }
+}
 
-__kernel void cpy2d_c2r(
-                    __global float* d_sino,
-                    __global float2* d_sino_complex,
-                    int num_bins,
-                    int num_projs,
-                    int fft_size)
+// copy only the real part of the valid data to the real array
+kernel void cpy2d_c2r(
+                      global float* d_sino,
+                      global float2* d_sino_complex,
+                      int num_bins,
+                      int num_projs,
+                      int fft_size)
 {
     int gid0 = get_global_id(0);
     int gid1 = get_global_id(1);
diff --git a/silx/resources/opencl/proj.cl b/silx/resources/opencl/proj.cl
index e64f63b2a..afc58ffc5 100644
--- a/silx/resources/opencl/proj.cl
+++ b/silx/resources/opencl/proj.cl
@@ -28,19 +28,18 @@
 /************************ GPU VERSION (with textures) **************************/
 /*******************************************************************************/
 
-const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;
 
-__kernel void  forward_kernel(
-        __global float *d_Sino,
-        __read_only image2d_t d_slice,
+kernel void  forward_kernel(
+        global float *d_Sino,
+        read_only image2d_t d_slice,
         int dimslice,
         int num_bins,
-        __global float* angles_per_project ,
+        global float* angles_per_project ,
         float axis_position,
-        __global float *d_axis_corrections,
-        __global int* d_beginPos    ,
-        __global int* d_strideJoseph,
-        __global int* d_strideLine  ,
+        global float *d_axis_corrections,
+        global int* d_beginPos    ,
+        global int* d_strideJoseph,
+        global int* d_strideLine  ,
         int num_projections,
         int  dimrecx,
         int  dimrecy,
@@ -49,7 +48,7 @@ __kernel void  forward_kernel(
         int josephnoclip,
         int normalize)
 {
-
+    const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;
     const int tidx = get_local_id(0);
     const int bidx = get_group_id(0);
     const int tidy = get_local_id(1);
@@ -57,10 +56,10 @@ __kernel void  forward_kernel(
     float angle;
     float cos_angle,sin_angle ;
 
-    __local float corrections[16];
-    __local int beginPos[16*2];
-    __local int strideJoseph[16*2];
-    __local int strideLine[16*2];
+    local float corrections[16];
+    local int beginPos[16*2];
+    local int strideJoseph[16*2];
+    local int strideLine[16*2];
 
     // thread will use corrections[tidy]
     // All are read by first warp
@@ -159,7 +158,8 @@ __kernel void  forward_kernel(
 
     if((bidy*16 + tidy) < num_projections && (bidx*16 + tidx) < num_bins) {
         res *= Area;
-        if (normalize) res *= M_PI*0.5/num_projections;
+        if (normalize)
+            res *= M_PI_F * 0.5f / num_projections;
         d_Sino[dimrecx*(bidy*16 + tidy) + (bidx*16 + tidx)] = res;
     }
 }
@@ -171,17 +171,17 @@ __kernel void  forward_kernel(
 /*******************************************************************************/
 
 
-__kernel void  forward_kernel_cpu(
-        __global float *d_Sino,
-        __global float* d_slice,
+kernel void  forward_kernel_cpu(
+        global float *d_Sino,
+        global float* d_slice,
         int dimslice,
         int num_bins,
-        __global float* angles_per_project ,
+        global float* angles_per_project ,
         float axis_position,
-        __global float *d_axis_corrections,
-        __global int* d_beginPos    ,
-        __global int* d_strideJoseph,
-        __global int* d_strideLine  ,
+        global float *d_axis_corrections,
+        global int* d_beginPos    ,
+        global int* d_strideJoseph,
+        global int* d_strideLine  ,
         int num_projections,
         int  dimrecx,
         int  dimrecy,
@@ -198,10 +198,10 @@ __kernel void  forward_kernel_cpu(
     float angle;
     float cos_angle,sin_angle ;
 
-    __local float corrections[16];
-    __local int beginPos[16*2];
-    __local int strideJoseph[16*2];
-    __local int strideLine[16*2];
+    local float corrections[16];
+    local int beginPos[16*2];
+    local int strideJoseph[16*2];
+    local int strideLine[16*2];
 
     // thread will use corrections[tidy]
     // All are read by first warp
@@ -338,7 +338,8 @@ __kernel void  forward_kernel_cpu(
 
     if((bidy*16 + tidy) < num_projections && (bidx*16 + tidx) < num_bins) {
         res *= Area;
-        if (normalize) res *= M_PI*0.5/num_projections;
+        if (normalize)
+            res *= M_PI_F * 0.5f / num_projections;
         d_Sino[dimrecx*(bidy*16 + tidy) + (bidx*16 + tidx)] = res;
     }
 }
['silx/resources/opencl/addition.cl', 'silx/opencl/common.py', 'silx/opencl/test/test_backprojection.py', 'silx/resources/opencl/backproj.cl', 'silx/resources/opencl/backproj_helper.cl', 'silx/resources/opencl/array_utils.cl', 'silx/opencl/backprojection.py', 'silx/resources/opencl/proj.cl', 'silx/opencl/processing.py']