def __init__(self, queue, discr, dtype, allocator):
        context = queue.context
        self.discr = discr
        import pyopencl as cl

        self.allocator = allocator

        dtype4 = cl.array.vec.types[np.dtype(dtype), 4]

        l = discr.ldis
        drdsdt_unvec = np.zeros((l.Np, l.Np, 4), dtype)
        for i, mat in enumerate([l.Dr, l.Ds, l.Dt]):
            drdsdt_unvec[:, :, i] = mat

        self.drdsdt = cl.array.to_device(
            queue,
            drdsdt_unvec
            .view(dtype=dtype4)[:, :, 0].copy(order="F"))
        self.drdsdt_img = cl.image_from_array(context, drdsdt_unvec.view(dtype=dtype4)[:, :, 0])

        drst_dx_unvec = np.zeros((discr.K, 4), dtype)
        drst_dy_unvec = np.zeros((discr.K, 4), dtype)
        drst_dz_unvec = np.zeros((discr.K, 4), dtype)

        for i in range(3):
            drst_dx_unvec[:, i] = discr.drst_dxyz[i, 0][:,0]
            drst_dy_unvec[:, i] = discr.drst_dxyz[i, 1][:,0]
            drst_dz_unvec[:, i] = discr.drst_dxyz[i, 2][:,0]

        self.drst_dx = cl.array.to_device(queue, drst_dx_unvec.view(dtype=dtype4)[:, 0])
        self.drst_dy = cl.array.to_device(queue, drst_dy_unvec.view(dtype=dtype4)[:, 0])
        self.drst_dz = cl.array.to_device(queue, drst_dz_unvec.view(dtype=dtype4)[:, 0])

        self.vmapP = cl.array.to_device(queue,
                discr.vmapP.astype(np.int32).copy().reshape(discr.K, -1))
        self.vmapM = cl.array.to_device(queue,
                discr.vmapM.astype(np.int32).copy().reshape(discr.K, -1))

        self.nx = cl.array.to_device(queue, discr.nx.astype(dtype))
        self.ny = cl.array.to_device(queue, discr.ny.astype(dtype))
        self.nz = cl.array.to_device(queue, discr.nz.astype(dtype))

        self.Fscale = cl.array.to_device(queue, discr.Fscale.astype(dtype))
        self.bc = cl.array.to_device(queue, discr.bc.astype(dtype))

        self.LIFT = cl.array.to_device(queue, l.LIFT.copy(order="F").astype(dtype))
        self.LIFT_img = cl.image_from_array(context, l.LIFT.astype(dtype))

        self.volume_events = []
        self.surface_events = []

    def __call__(self, ctx, src, kernel):
        self.build(ctx)
        kernel = np.array(kernel, copy=False, dtype=np.float32)
        halflen = kernel.shape[0] / 2
        kernelf = kernel.flatten()
        kernelf_buf = cl.Buffer(self.ctx, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=kernelf)

        src_padded = np.zeros((src.shape[0]+2*halflen, src.shape[1]+2*halflen, 4), dtype=src.dtype)
        src_padded[halflen:-halflen,halflen:-halflen,:src.shape[2]] = src[:,:,:src.shape[2]]

        src_padded[halflen:-halflen,:halflen,:src.shape[2]] = src_padded[halflen:-halflen,halflen:halflen*2,:src.shape[2]][:,::-1]
        src_padded[halflen:-halflen,-halflen:,:src.shape[2]] = src_padded[halflen:-halflen,-halflen*2:-halflen,:src.shape[2]][:,::-1]

        src_padded[:halflen,:,:src.shape[2]] = src_padded[halflen:halflen*2,:,:src.shape[2]][::-1,...]
        src_padded[-halflen:,:,:src.shape[2]] = src_padded[-halflen*2:-halflen,:,:src.shape[2]][::-1,...]

        norm = np.issubdtype(src.dtype, np.integer)
        src_buf = cl.image_from_array(self.ctx, src_padded, 4, norm_int=norm)
        dest = np.zeros((src.shape[0], src.shape[1], 4), dtype=src.dtype)
        dest_buf = init_image(self.ctx, dest, 4, mode="w", norm_int=norm)

        queue = cl.CommandQueue(self.ctx)
        self.prg.convolve2d_naive(queue, (dest.shape[1], dest.shape[0]), None, src_buf, dest_buf, kernelf_buf, np.int32(kernel.shape[0]))
        cl.enqueue_copy(queue, dest, dest_buf, origin=(0, 0), region=(src.shape[1], src.shape[0])).wait()

        # src_buf.release()
        # dest_buf.release()
        # kernelf_buf.release()
        return dest[:,:,:src.shape[2]]

    def test_rotate_image3d_1(self):
        
        shape = (8, 6, 5)
        np_image = np.zeros(shape, dtype=np.float32)
        np_image[0, 0, 0] = 1
        np_image[0, 0, 1] = 1
        np_image[0, 0, 2] = 1

        # 90 degree rotation around z-axis
        rotmat = [[0, 1, 0],
                  [1, 0, 0],
                  [0, 0, 1]]

        np_out = np.zeros_like(np_image)

        expected = np.zeros_like(np_image)
        expected[0, 0, 0] = 1
        expected[0, 1, 0] = 1
        expected[0, 2, 0] = 1

        cl_image = cl.image_from_array(self.queue.context, np_image)
        cl_out = cl_array.to_device(self.queue, np_out)
        cl_sampler = cl.Sampler(self.queue.context, False, cl.addressing_mode.CLAMP, cl.filter_mode.LINEAR)

        self.kernels.rotate_image3d(self.queue, cl_sampler, cl_image, rotmat, cl_out)

        self.assertTrue(np.allclose(expected, cl_out.get()))

 def mask(self, mask):
     BaseCorrelator.mask.fset(self, mask)
     self._norm_factor = np.float32(self._norm_factor)
     self._rmax = np.int32(self._rmax)
     self._gtemplate = cl.image_from_array(
             self._ctx, self._template.astype(np.float32)
             )
     self._gmask = cl.image_from_array(
             self._ctx, self._mask.astype(np.float32)
             )
     max_items = self._queue.device.max_compute_units * 32 * 16
     gws = [0] * 3
     gws[0] = min(2 * self._rmax, max_items)
     gws[1] = min(max_items // gws[0], 2 * self._rmax)
     gws[2] = min(max(max_items // (gws[0] * gws[0]), 1), 2 * self._rmax)
     self._gws = tuple(gws)

    def test_clashvol(self):

        NROT = np.random.randint(self.rotations.shape[0] + 1)
        rotmat = self.rotations[NROT]
        cpu_lsurf = np.zeros_like(self.im_lsurf.array)
        disvis.libdisvis.rotate_image3d(self.im_lsurf.array, self.vlength, np.linalg.inv(rotmat), self.im_center, cpu_lsurf)

        cpu_clashvol = numpy.fft.irfftn(numpy.fft.rfftn(cpu_lsurf).conj() * numpy.fft.rfftn(self.rcore.array), s=self.shape)

        gpu_rcore = cl_array.to_device(self.queue, np.asarray(self.rcore.array, dtype=np.float32))
        gpu_im_lsurf = cl.image_from_array(self.queue.context, np.asarray(self.im_lsurf.array, dtype=np.float32))
        gpu_lsurf = cl_array.zeros(self.queue, self.shape, dtype=np.float32)

        self.kernels.rotate_image3d(self.queue, self.sampler, gpu_im_lsurf, rotmat, gpu_lsurf, self.im_center)

        gpu_ft_lsurf = cl_array.zeros(self.queue, self.ft_shape, dtype=np.complex64)
        gpu_ft_rcore = cl_array.zeros(self.queue, self.ft_shape, dtype=np.complex64)
        gpu_ft_clashvol = cl_array.zeros(self.queue, self.ft_shape, dtype=np.complex64)
        gpu_clashvol = cl_array.zeros(self.queue, self.shape, dtype=np.float32)

        self.kernels.rfftn(self.queue, gpu_rcore, gpu_ft_rcore)
        self.kernels.rfftn(self.queue, gpu_lsurf, gpu_ft_lsurf)
        self.kernels.c_conj_multiply(self.queue, gpu_ft_lsurf, gpu_ft_rcore, gpu_ft_clashvol)
        self.kernels.irfftn(self.queue, gpu_ft_clashvol, gpu_clashvol)

        self.assertTrue(np.allclose(cpu_clashvol, gpu_clashvol.get(), atol=0.8))

    def _upload_image(self, image):
        assert image.max() <= 1.0

        # Check the number of channels in the image
        if image.ndim == 2:
            num_channels = 1
        else:
            if sys.platform.startswith('win') and 'geforce' in self.ctx.devices[0].name.lower() and image.shape[2] == 3:
                # This is a hack for Windows/nVidia, as we believe and found so
                # far for various GeFoce cards that the nvidia OpenCL
                # implementation sucks. Reporting an out-of-resources error when
                # trying to upload an RGB three channel image to the GPU
                # Quite counterintuitively adding an unneeded fourth channel
                # makes the out-of-resources error go away. FIXME if you can.
                tmp = image
                image = np.ones((tmp.shape[0], tmp.shape[1], 4))
                num_channels = 4
                image[:, :, :3] = tmp[:]
            else:
                num_channels = image.shape[2]

        # Tell OpenCL to copy the image into device memory
        image_gpu = cl.image_from_array(self.ctx, image.astype(np.float32),
                                        num_channels=num_channels, mode="r")
        return image_gpu

 def __call__(self, ctx, src):
     self.build(ctx)
     src = np.asarray(src)
     src2 = np.zeros((src.shape[0], src.shape[1], 4),dtype=src.dtype)
     src2[:,:,0:src.shape[2]] = src[:,:,0:src.shape[2]]
     norm = np.issubdtype(src2.dtype, np.integer)
     src2_buf = cl.image_from_array(self.ctx, src2, 4, norm_int=norm)
     dest_buf = cl.image_from_array(self.ctx, src2, 4, mode="w", norm_int=norm)
     dest = np.empty_like(src2)
     queue = cl.CommandQueue(self.ctx)
     self.prg.YCrCb2RGB(queue, (src2.shape[1], src2.shape[0]), None, src2_buf, dest_buf)
     cl.enqueue_copy(queue, dest, dest_buf, origin=(0, 0), region=(src2.shape[1], src2.shape[0])).wait()
     dest = dest[:,:,0:src.shape[2]].copy()
     src2_buf.release()
     dest_buf.release()
     return dest

	def get_color(self, img):
		# OpenCL only supports RGBA images, not RGB, so add an alpha channel
		src = np.array(img.convert('RGBA'))
		src.shape = w, h, _ = img.width, img.height, 4

		w = int(w * self.SCALE_FACTOR)
		h = int(h * self.SCALE_FACTOR)

		local_size = self.max_work_item_sizes
		global_size = (math.ceil(h / local_size[0]), math.ceil(w / local_size[1]))

		total_work_groups = global_size[0] * global_size[1]

		mf = cl.mem_flags
		src_buf = cl.image_from_array(self.ctx, src, 4, norm_int=True)

		out = np.zeros(4 * total_work_groups, dtype=np.int32)
		out_buf = cl.Buffer(self.ctx, mf.WRITE_ONLY, size=out.itemsize * 4 * total_work_groups)

		kernel = self.prg.get_color
		kernel.set_scalar_arg_dtypes([None, None, np.uint32, np.uint32])
		kernel(self.queue, global_size, local_size, src_buf, out_buf, w, h, g_times_l=True)

		cl.enqueue_copy(self.queue, dest=out, src=out_buf, is_blocking=True)

		# this sum takes .1 ms at 3440x1440, don't even bother OpenCL-ifying it
		resized_out = np.reshape(out, (out.shape[0] / 4, 4))
		summed_out = np.sum(resized_out, axis=0)

		avg_out = (summed_out / summed_out[3])[:3].astype(int)

		return avg_out

    def find_starburst_ray_boundaries(self, im, seed_point, cutoff_index,
                                      threshold, n_rays, n_samples, ray_step):

        if self.cached_shape != im.shape:
            self.setup_device(im.shape)

        #(im_, _, _) = self.sobel3x3_separable(im.astype(np.float32))
        im_ = im.astype(np.float32)
        self.clIm2D = cl.image_from_array(self.ctx, im_, num_channels=1)
        # # load im to memory
        # cl.enqueue_copy(self.q, self.clIm2D, clIm.data, offset=0,
        #                 origin=(0, 0), region=clIm.shape)

        seed_point_ = (seed_point[1], seed_point[0])

        # sample the rays, computing the "ray-wise" gradient and mean + std along the
        # way.  We'll pull back the mean and running stds and compute the thresholds
        # on the CPU
        sampled = self.cl_find_ray_boundaries(self.clIm2D, n_rays,
                                     n_samples, ray_step, seed_point_, cutoff_index, threshold)

        # run through the resampled values to find cutoffs

        # pull back the cutoff and return them
        return sampled

    def _gpu_init(self):
        """Method to initialize all the data for GPU-accelerate search"""

        self.gpu_data = {}
        g = self.gpu_data
        d = self.data
        q = self.queue

        # move data to the GPU. All should be float32, as these is the native
        # lenght for GPUs
        g['rcore'] = cl_array.to_device(q, float32array(d['rcore'].array))
        g['rsurf'] = cl_array.to_device(q, float32array(d['rsurf'].array))
        # Make the scanning chain object an Image, as this is faster to rotate
        g['im_lsurf'] = cl.image_from_array(q.context, float32array(d['lsurf'].array))
        g['sampler'] = cl.Sampler(q.context, False, cl.addressing_mode.CLAMP,
                                  cl.filter_mode.LINEAR)

        if self.distance_restraints:
            g['restraints'] = cl_array.to_device(q, float32array(d['restraints']))

        # Allocate arrays on the GPU
        g['lsurf'] = cl_array.zeros_like(g['rcore'])
        g['clashvol'] = cl_array.zeros_like(g['rcore'])
        g['intervol'] = cl_array.zeros_like(g['rcore'])
        g['interspace'] = cl_array.zeros(q, d['shape'], dtype=np.int32)
        g['restspace'] = cl_array.zeros_like(g['interspace'])
        g['access_interspace'] = cl_array.zeros_like(g['interspace'])
        g['best_access_interspace'] = cl_array.zeros_like(g['interspace'])

        # arrays for counting
        # Reductions are typically tedious on GPU, and we need to define the
        # workgroupsize to allocate the correct amount of data
        WORKGROUPSIZE = 32
        nsubhists = int(np.ceil(g['rcore'].size/WORKGROUPSIZE))
        g['subhists'] = cl_array.zeros(q, (nsubhists, d['nrestraints'] + 1), dtype=np.float32)
        g['viol_counter'] = cl_array.zeros(q, (nsubhists, d['nrestraints'], d['nrestraints']), dtype=np.float32)

        # complex arrays
        g['ft_shape'] = list(d['shape'])
        g['ft_shape'][0] = d['shape'][0]//2 + 1
        g['ft_rcore'] = cl_array.zeros(q, g['ft_shape'], dtype=np.complex64)
        g['ft_rsurf'] = cl_array.zeros_like(g['ft_rcore'])
        g['ft_lsurf'] = cl_array.zeros_like(g['ft_rcore'])
        g['ft_clashvol'] = cl_array.zeros_like(g['ft_rcore'])
        g['ft_intervol'] = cl_array.zeros_like(g['ft_rcore'])

        # other miscellanious data
        g['nrot'] = d['nrot']
        g['max_clash'] = d['max_clash']
        g['min_interaction'] = d['min_interaction']

        # kernels
        g['k'] = Kernels(q.context)
        g['k'].rfftn = pyclfft.RFFTn(q.context, d['shape'])
        g['k'].irfftn = pyclfft.iRFFTn(q.context, d['shape'])

        # initial calculations
        g['k'].rfftn(q, g['rcore'], g['ft_rcore'])
        g['k'].rfftn(q, g['rsurf'], g['ft_rsurf'])

    def test_touch(self):

        MAX_CLASH = 100 + 0.9
        MIN_INTER = 300 + 0.9

        NROT = np.random.randint(self.rotations.shape[0] + 1)
        rotmat = self.rotations[0]
        cpu_lsurf = np.zeros_like(self.im_lsurf.array)
        disvis.libdisvis.rotate_image3d(self.im_lsurf.array, self.vlength, np.linalg.inv(rotmat), self.im_center, cpu_lsurf)

        cpu_clashvol = numpy.fft.irfftn(numpy.fft.rfftn(cpu_lsurf).conj() * numpy.fft.rfftn(self.rcore.array))

        gpu_rcore = cl_array.to_device(self.queue, np.asarray(self.rcore.array, dtype=np.float32))
        gpu_im_lsurf = cl.image_from_array(self.queue.context, np.asarray(self.im_lsurf.array, dtype=np.float32))
        gpu_lsurf = cl_array.zeros(self.queue, self.shape, dtype=np.float32)

        self.kernels.rotate_image3d(self.queue, self.sampler, gpu_im_lsurf, rotmat, gpu_lsurf, self.im_center)

        gpu_ft_lsurf = cl_array.zeros(self.queue, self.ft_shape, dtype=np.complex64)
        gpu_ft_rcore = cl_array.zeros(self.queue, self.ft_shape, dtype=np.complex64)
        gpu_ft_clashvol = cl_array.zeros(self.queue, self.ft_shape, dtype=np.complex64)
        gpu_clashvol = cl_array.zeros(self.queue, self.shape, dtype=np.float32)

        self.kernels.rfftn(self.queue, gpu_rcore, gpu_ft_rcore)
        self.kernels.rfftn(self.queue, gpu_lsurf, gpu_ft_lsurf)
        self.kernels.c_conj_multiply(self.queue, gpu_ft_lsurf, gpu_ft_rcore, gpu_ft_clashvol)
        self.kernels.irfftn(self.queue, gpu_ft_clashvol, gpu_clashvol)
        
        cpu_intervol = numpy.fft.irfftn(numpy.fft.rfftn(cpu_lsurf).conj() * numpy.fft.rfftn(self.rsurf.array))

        gpu_rsurf = cl_array.to_device(self.queue, np.asarray(self.rsurf.array, dtype=np.float32))

        gpu_ft_rsurf = cl_array.zeros(self.queue, self.ft_shape, dtype=np.complex64)
        gpu_ft_intervol = cl_array.zeros(self.queue, self.ft_shape, dtype=np.complex64)
        gpu_intervol = cl_array.zeros(self.queue, self.shape, dtype=np.float32)

        cpu_interspace = np.zeros(self.shape, dtype=np.int32)
        gpu_interspace = cl_array.zeros(self.queue, self.shape, dtype=np.int32)

        self.kernels.rfftn(self.queue, gpu_rsurf, gpu_ft_rsurf)
        self.kernels.rfftn(self.queue, gpu_lsurf, gpu_ft_lsurf)
        self.kernels.c_conj_multiply(self.queue, gpu_ft_lsurf, gpu_ft_rsurf, gpu_ft_intervol)
        self.kernels.irfftn(self.queue, gpu_ft_intervol, gpu_intervol)

        self.kernels.touch(self.queue, gpu_clashvol, MAX_CLASH, gpu_intervol, MIN_INTER, gpu_interspace)

        np.logical_and(cpu_clashvol < MAX_CLASH, cpu_intervol > MIN_INTER, cpu_interspace)

        disvis.volume.Volume(cpu_interspace, self.im_lsurf.voxelspacing, self.im_lsurf.origin).tofile('cpu_interspace.mrc')
        disvis.volume.Volume(gpu_interspace.get(), self.im_lsurf.voxelspacing, self.im_lsurf.origin).tofile('gpu_interspace.mrc')
        disvis.volume.Volume(cpu_interspace - gpu_interspace.get(), self.im_lsurf.voxelspacing, self.im_lsurf.origin).tofile('diff.mrc')
        print()
        print(cpu_interspace.sum(), gpu_interspace.get().sum())
        print(np.abs(cpu_interspace - gpu_interspace.get()).sum())
                           

        self.assertTrue(np.allclose(gpu_interspace.get(), cpu_interspace))

    def cl_load_data(self, population, world):
        mf = cl.mem_flags


        out = cl.Buffer(self.ctx, mf.WRITE_ONLY, population.nbytes)

        population_cl = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=population)

#        world_cl = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=world.flatten())
        world_cl = cl.image_from_array(self.ctx, world, mode="r")
        
        return population_cl, world_cl, out

  def __init__(self,volumeNode,contextPreference='GPU',renderSize=(512,512)):
    self.volumeNode = volumeNode
    self.volumeArray = slicer.util.array(self.volumeNode.GetID())
    self.renderSize = renderSize

    try:
      import pyopencl
      import numpy
    except ImportError:
      raise "No OpenCL for you!\nInstall pyopencl in slicer's python installation."

    import os
    os.environ['PYOPENCL_COMPILER_OUTPUT'] = '1'

    self.ctx = None
    for platform in pyopencl.get_platforms():
        for device in platform.get_devices():
            if pyopencl.device_type.to_string(device.type) == contextPreference:
               self.ctx = pyopencl.Context([device])
               break;

    if not self.ctx:
      self.ctx = pyopencl.create_some_context()
    self.queue = pyopencl.CommandQueue(self.ctx)

    inPath = os.path.dirname(slicer.modules.rendercl.path) + "/Render.cl.in"
    fp = open(inPath)
    sourceIn = fp.read()
    fp.close()

    source = sourceIn % { 
        'rayStepSize' : '0.01f',
        'rayMaxSteps' : '500',
        }
    self.prg = pyopencl.Program(self.ctx, source).build()

    # create a 3d image from the volume
    num_channels = 1
    self.volumeImage_dev = pyopencl.image_from_array(self.ctx, self.volumeArray, num_channels)

    # create a 2d array for the render buffer
    self.renderArray = numpy.zeros(self.renderSize,dtype=numpy.dtype('uint32'))
    self.renderArray_dev = pyopencl.array.to_device(self.queue, self.renderArray)

    self.volumeSampler = pyopencl.Sampler(self.ctx,False,
                              pyopencl.addressing_mode.REPEAT,
                              pyopencl.filter_mode.LINEAR)

    # TODO make 2D image of transfer function

    self.transferFunctionSampler = pyopencl.Sampler(self.ctx,False,
                              pyopencl.addressing_mode.REPEAT,
                              pyopencl.filter_mode.LINEAR)

 def from_array(cls,arr, *args, **kwargs):
     ctx = get_device().context
     if not arr.ndim in [1,2,3,4]:
         raise ValueError("dimension of array wrong, should be 1...4 but is %s"%arr.ndim)        
     elif arr.ndim == 4:
         num_channels = arr.shape[-1]
     else:
         num_channels = None
     
     res =  pyopencl.image_from_array(ctx, arr,num_channels = num_channels,
                                          *args, **kwargs)
     res.dtype = arr.dtype
     return res

 def __call__(self, ctx, src2, kernel):
     if self.ctx != ctx:
         self.ctx = ctx
         self.prg = cl.Program(self.ctx, pkg_resources.resource_string(__name__, "convolve2d.cl")).build()
     src2 = np.asarray(src2)
     src = np.zeros((src2.shape[0], src2.shape[1], 4),dtype=src2.dtype)
     src[:,:,0:src2.shape[2]] = src2[:,:,0:src2.shape[2]]
     norm = np.issubdtype(src.dtype, np.integer)
     src_buf = cl.image_from_array(self.ctx, src, 4, norm_int=norm)
     dest_buf = cl.image_from_array(self.ctx, src, 4, mode="w", norm_int=norm)
     dest = np.empty_like(src)
     kernel = np.array(kernel, dtype=np.float32)
     kernelf = kernel.flatten()
     kernelf_buf = cl.Buffer(self.ctx, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=kernelf)
     halflen = (kernelf.shape[0]>>1)
     queue = cl.CommandQueue(self.ctx)
     self.prg.convolve2d_local(queue, (src.shape[1]-halflen, src.shape[0]-halflen), None, src_buf, dest_buf, kernelf_buf, np.int_(kernelf.shape[0]))
     cl.enqueue_copy(queue, dest, dest_buf, origin=(0, 0), region=(src.shape[1], src.shape[0])).wait()
     dest = dest[:,:,0:src2.shape[2]].copy()
     src_buf.release()
     dest_buf.release()
     kernelf_buf.release()
     return dest

    def __call__(self, ctx, pyramid, blendmap, x, y):
        self.build(ctx)

        norm = np.issubdtype(pyramid.dtype, np.integer)

        pyramid_buf = cl.image_from_array(self.ctx, pyramid, 4, mode="r", norm_int=norm)

        dest = np.zeros_like(pyramid[:,:,0,:],dtype=pyramid.dtype)
        dest_buf = init_image(self.ctx, dest, 4, mode="w", norm_int=norm)

        xoff = dest.shape[1] - x
        yoff = dest.shape[0] - y
        blendmap_buf = cl.image_from_array(self.ctx, blendmap[yoff:yoff+dest.shape[0],xoff:xoff+dest.shape[1],:].copy(), 4, mode="r")

        queue = cl.CommandQueue(self.ctx)
        self.prg.blend(queue, (dest.shape[1], dest.shape[0]), None,
                       pyramid_buf, blendmap_buf, dest_buf)
        cl.enqueue_copy(queue, dest, dest_buf, origin=(0, 0), region=(dest.shape[1], dest.shape[0])).wait()

        dest_buf.release()
        pyramid_buf.release()
        blendmap_buf.release()
        return dest

	def popCorn(self):
		mf = cl.mem_flags

		#initialize client side (CPU) arrays
		#self.a = numpy.array(range(10), dtype=numpy.float32)
		#self.b = numpy.array(range(10), dtype=numpy.float32)
		self.data = pyfits.getdata(os.path.join(self.path,self.file1))

		#create OpenCL buffers
		#self.imagebuf = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.ra1)
		
		self.imagebuf = cl.image_from_array(self.ctx,self.data)
		
		self.dest_buf = cl.Buffer(self.ctx, mf.WRITE_ONLY, self.dec2.nbytes)

    def test_rotate_grids_and_multiply(self):
        shape = (5, 5, 5)
        template = np.zeros(shape, dtype=np.float32)
        template[2, 2, 1:4] = 1
        template[2, 1:4, 2] = 1
        mask = template * 2
        np_out_template = np.zeros(shape, dtype=np.float32)
        np_out_template[0, 0, :2] = 1
        np_out_template[0, 0, -1] = 1
        np_out_template[0, :2, 0] = 1
        np_out_template[0, -1, 0] = 1
        np_out_mask = np_out_template * 2 
        np_out_mask2 = np_out_mask ** 2


        cl_template = cl.image_from_array(self.ctx, template)
        cl_mask = cl.image_from_array(self.ctx, mask)
        cl_rotmat = np.asarray([1, 0, 0, 0, 1, 0, 0, 0, 1] + [0] * 7, dtype=np.float32)
        cl_center = np.asarray([2, 2, 2, 0], dtype=np.float32)
        cl_shape = np.asarray([5, 5, 5, 125], dtype=np.int32)
        cl_radius = np.int32(2)

        cl_out_template = cl_array.to_device(self.queue, np.zeros(shape, dtype=np.float32))
        cl_out_mask = cl_array.to_device(self.queue, np.zeros(shape, dtype=np.float32))
        cl_out_mask2 = cl_array.to_device(self.queue, np.zeros(shape, dtype=np.float32))

        gws = tuple([int(2 * cl_radius + 1)] * 3)
        args = (cl_template, cl_mask, cl_rotmat, self.s_linear, self.s_nearest,
                cl_center, cl_shape, cl_radius, cl_out_template.data, cl_out_mask.data,
                cl_out_mask2.data)
        self.k.rotate_grids_and_multiply(self.queue, gws, None, *args)
        self.queue.finish()

        self.assertTrue(np.allclose(np_out_template, cl_out_template.get()))
        self.assertTrue(np.allclose(np_out_mask, cl_out_mask.get()))
        self.assertTrue(np.allclose(np_out_mask2, cl_out_mask2.get()))

    def test_rotate_image3d(self):
        # CPU
        NROT = np.random.randint(self.rotations.shape[0] + 1)
        rotmat = self.rotations[NROT]

        #print(rotmat)

        cpu_lsurf = np.zeros_like(self.im_lsurf.array)
        disvis.libdisvis.rotate_image3d(self.im_lsurf.array, self.vlength, np.linalg.inv(rotmat), self.im_center, cpu_lsurf)

        gpu_im_lsurf = cl.image_from_array(self.queue.context, np.asarray(self.im_lsurf.array, dtype=np.float32))
        gpu_lsurf = cl_array.zeros(self.queue, self.shape, dtype=np.float32)
        self.kernels.rotate_image3d(self.queue, self.sampler, gpu_im_lsurf, rotmat, gpu_lsurf, self.im_center)

        self.assertTrue(np.allclose(cpu_lsurf, gpu_lsurf.get(), atol=0.01))

def LoadImage(context, fileName):
    im = Image.open(fileName)
    img = np.array(im)
    IMG1 = scale_img(img, 8)
    im = Image.fromarray(IMG1)
    # Make sure the image is RGBA formatted
    if im.mode != "RGBA":
        im = im.convert("RGBA")
    IMG1 = np.array(im)
    if len(IMG1.shape) > 2:
        nchannels = IMG1.shape[-1]
    else:
        nchannels = None
    t0 = time.clock()
    clImage = cl.image_from_array(context, IMG1, num_channels=nchannels,
                                  mode="r", norm_int=False)
    t1 = time.clock()
    print t1 - t0, " Load to GPU..."
    return clImage, im.size, IMG1