Python gputools.OCLProgram类代码示例

def bilateral3(data, size_filter, sigma_p, sigma_x = 10.):
    """bilateral filter """
    
    dtype = data.dtype.type
    dtypes_kernels = {np.float32:"bilat3_float",}

    if not dtype in dtypes_kernels.keys():
        logger.info("data type %s not supported yet (%s), casting to float:"%(dtype,dtypes_kernels.keys()))
        data = data.astype(np.float32)
        dtype = data.dtype.type


    img = OCLImage.from_array(data)
    res = OCLArray.empty_like(data)

    
    prog = OCLProgram(abspath("kernels/bilateral3.cl"))

    print img.shape

    prog.run_kernel(dtypes_kernels[dtype],
                    img.shape,None,
                    img,res.data,
                    np.int32(img.shape[0]),np.int32(img.shape[1]),
                    np.int32(size_filter),np.float32(sigma_x),np.float32(sigma_p))


    return res.get()

def gpu_kuwahara(data, N=5):
    """Function to convolve an imgage with the Kuwahara filter on GPU."""
    # create numpy arrays


    if (N%2==0):       
        raise ValueError("Data has to be a (2n+1)x(2n+1) array.")

    
    data_g = OCLArray.from_array(data.astype(float32)) 
       
    res_g = OCLArray.empty((data.shape[0],data.shape[1]),float32) 
    
    prog = OCLProgram("./OpenCL/gpu_kernels/gpu_kuwahara.cl")
    
    # start kernel on gput
    prog.run_kernel("kuwahara",   # the name of the kernel in the cl file
                   data_g.shape[::-1], # global size, the number of threads e.g. (128,128,) 
                    None,   # local size, just leave it to None
                    data_g.data,res_g.data,
                    int32(N)) 
                    
    
#                    
    
    return res_g.get()

def create_dn_buffer(size, units,points,
                     dn_inner = .0, rad_inner = 0,
                     dn_outer = .1, rad_outer = .4):

    Nx, Ny, Nz = size
    dx, dy, dz = units

    program = OCLProgram(absPath("kernels/bpm_3d_spheres.cl"))


    dn_g = OCLArray.empty((Nz,Ny,Nx),dtype=np.float32)

    # sort by z
    ps = np.array(points)
    ps = ps[np.argsort(ps[:,2]),:]

    Np = ps.shape[0]

    pointsBuf = OCLArray.from_array(ps.flatten().astype(np.float32))

    program.run_kernel("fill_dn",(Nx,Ny,Nz),None,dn_g.data,
                       pointsBuf.data,np.int32(Np),
                       np.float32(dx),np.float32(dy),np.float32(dz),
                       np.float32(dn_inner),np.float32(rad_inner),
                       np.float32(dn_outer),np.float32(rad_outer))


    return dn_g

def test_bessel(n,x):
    x_g = OCLArray.from_array(x.astype(float32))
    res_g = OCLArray.empty_like(x.astype(float32))
    
    p = OCLProgram(absPath("kernels/bessel.cl"))
    p.run_kernel("bessel_fill",x_g.shape,None,
                 x_g.data,res_g.data,int32(n))

    return res_g.get()

def perlin2(size, units, repeat = (10.,)*2):
    wx, wy = repeat
    dx, dy = units

    prog = OCLProgram(abspath("perlin.cl"))

    d = OCLArray.empty(size[::-1],np.float32)
    prog.run_kernel("perlin2d",d.shape[::-1],None,
                    d.data,
                    np.float32(dx),np.float32(dy),
                    np.float32(wx),np.float32(wy))

    return d.get()

def focus_field_cylindrical_plane(shape = (128,128),
                            units = (.1,.1),
                            z = 0.,
                            lam = .5, NA = .6, n0 = 1.,
                            ex_g = None,
                            n_integration_steps = 200):
    """
    calculates the x component of the electric field  at a given z position z for a perfect, aberration free optical system
    via the vectorial debye diffraction integral for a cylindrical lens

    see
    Colin J. R. Sheppard,
    Cylindrical lenses—focusing and imaging: a review

    Appl. Opt. 52, 538-545 (2013)


    if ex_g is a valid OCLArray it fills it and returns None
    otherwise returns ex as a numpy array


    """

    p = OCLProgram(absPath("kernels/psf_cylindrical.cl"),build_options = str("-I %s -D INT_STEPS=%s"%(absPath("."),n_integration_steps)))


    Nx, Ny = shape
    dx, dy = units

    alpha = np.arcsin(NA/n0)

    if ex_g is None:
        use_buffer = False
        ex_g = OCLArray.empty((Ny,Nx),np.complex64)
    else:
        use_buffer = True

    assert ex_g.shape[::-1] == shape


    p.run_kernel("psf_cylindrical_plane",(Nx,Ny),None,
                 ex_g.data,
                 np.float32(-dy*(Ny-1)/2.),np.float32(dy*(Ny-1)/2.),
                 np.float32(z),
                 np.float32(lam/n0),
                 np.float32(alpha))

    if not use_buffer:
        return ex_g.get()

def focus_field_cylindrical(shape,units,lam = .5,NA = .3, n0=1.,
                            n_integration_steps = 100):
    """computes focus field of cylindrical lerns with given NA

    see:
    Colin J. R. Sheppard,
    Cylindrical lenses—focusing and imaging: a review

    Appl. Opt. 52, 538-545 (2013)

    return u,ex,ey,ez   with u being the intensity
    """

    p = OCLProgram(absPath("kernels/psf_cylindrical.cl"),build_options = str("-I %s -D INT_STEPS=%s"%(absPath("."),n_integration_steps)))

    
    Nx, Ny, Nz = shape
    dx, dy, dz = units

    alpha = np.arcsin(NA/n0)
    
    u_g = OCLArray.empty((Nz,Ny),np.float32)
    ex_g = OCLArray.empty((Nz,Ny),np.complex64)
    ey_g = OCLArray.empty((Nz,Ny),np.complex64)
    ez_g = OCLArray.empty((Nz,Ny),np.complex64)

    t = time.time()
    
    p.run_kernel("psf_cylindrical",u_g.shape[::-1],None,
                 ex_g.data,
                 ey_g.data,
                 ez_g.data,
                 u_g.data,
                 np.float32(-dy*(Ny-1)/2.),np.float32(dy*(Ny-1)/2.),
                 np.float32(-dz*(Nz-1)/2.),np.float32(dz*(Nz-1)/2.),
                 np.float32(lam/n0),
                 np.float32(alpha))

    u = np.array(np.repeat(u_g.get()[...,np.newaxis],Nx,axis=-1))
    ex = np.array(np.repeat(ex_g.get()[...,np.newaxis],Nx,axis=-1))
    ey = np.array(np.repeat(ey_g.get()[...,np.newaxis],Nx,axis=-1))
    ez = np.array(np.repeat(ez_g.get()[...,np.newaxis],Nx,axis=-1))

    
    print "time in secs:" , time.time()-t
    

    return u, ex, ey, ez

def _convolve_sep2_gpu(data_g, hx_g, hy_g, res_g = None):

    assert_bufs_type(np.float32,data_g,hx_g,hy_g)

    prog = OCLProgram(abspath("kernels/convolve_sep.cl"))

    Ny,Nx = hy_g.shape[0],hx_g.shape[0]

    tmp_g = OCLArray.empty_like(data_g)

    if res_g is None:
        res_g = OCLArray.empty_like(data_g)
    
    prog.run_kernel("conv_sep2_x",data_g.shape[::-1],None,data_g.data,hx_g.data,tmp_g.data,np.int32(Nx))
    prog.run_kernel("conv_sep2_y",data_g.shape[::-1],None,tmp_g.data,hy_g.data,res_g.data,np.int32(Ny))

    return res_g

def _perlin3_single(size,units = (1.,)*3,repeat = (10.,)*3,offz = 0,Nz0 = None):
    if Nz0 is None:
        Nz0 = size[-1]

    dx, dy, dz = units
    wx, wy, wz = repeat

    prog = OCLProgram(abspath("perlin.cl"))

    d = OCLArray.empty(size[::-1],np.float32)
    prog.run_kernel("perlin3d",d.shape[::-1],None,
                    d.data,
                    np.int32(offz),
                    np.float32(dx),np.float32(dy),np.float32(dz),
                    np.float32(wx),np.float32(wy),np.float32(wz) )

    return d.get()

def focus_field_debye_at(x,y,z,lam, NA, n0 = 1., n_integration_steps = 200):
    """ the same as focus_field_debye but for the coordinates given in x, y, z (arrays of same shape)

        slower than focus_field_debye as it doesnt assume the coordinates to be on a grid
    """

    print absPath("kernels/psf_debye.cl")
    p = OCLProgram(absPath("kernels/psf_debye.cl"),
                   build_options = str("-I %s -D INT_STEPS=%s"%(absPath("."),n_integration_steps)))

    if np.isscalar(NA):
        NA = [0.,NA]

    alphas = np.arcsin(np.array(NA)/n0)
    assert len(alphas)%2 ==0

    assert x.shape == y.shape == z.shape
    dshape =x.shape
    N = np.prod(dshape)

    x_g = OCLArray.from_array(x.flatten().astype(np.float32))
    y_g = OCLArray.from_array(y.flatten().astype(np.float32))
    z_g = OCLArray.from_array(z.flatten().astype(np.float32))

    u_g = OCLArray.empty(N,np.float32)
    ex_g = OCLArray.empty(N,np.complex64)
    ey_g = OCLArray.empty(N,np.complex64)
    ez_g = OCLArray.empty(N,np.complex64)

    alpha_g = OCLArray.from_array(alphas.astype(np.float32))

    p.run_kernel("debye_wolf_at",(N,),None,
                 x_g.data,y_g.data,z_g.data,
                 ex_g.data,ey_g.data,ez_g.data, u_g.data,
                 np.float32(1.),np.float32(0.),
                 np.float32(lam/n0),
                 alpha_g.data, np.int32(len(alphas)))

    u = u_g.get().reshape(dshape)
    ex = ex_g.get().reshape(dshape)
    ey = ey_g.get().reshape(dshape)
    ez = ez_g.get().reshape(dshape)

    return u, ex, ey, ez

def gpu_structure(data):
    """Function to convolve an imgage with a structure filter on GPU."""
    # create numpy arrays
    
    
    data_g = OCLArray.from_array(data.astype(float32)) 
       
    res_g = OCLArray.empty((data.shape[0],data.shape[1],2),float32) 
    
    prog = OCLProgram("./OpenCL/gpu_kernels/gpu_structure.cl")
    
    # start kernel on gput
    prog.run_kernel("structure",   # the name of the kernel in the cl file
                    data_g.shape[::-1], # global size, the number of threads e.g. (128,128,) 
                    None,   # local size, just leave it to None
                    data_g.data,res_g.data) 
                    
                    
    return res_g.get()

def _convolve_buf(data_g, h_g , res_g = None):
    """
    buffer variant
    """
    assert_bufs_type(np.float32,data_g,h_g)

    prog = OCLProgram(abspath("kernels/convolve.cl"))

    if res_g is None:
        res_g = OCLArray.empty(data_g.shape,dtype=np.float32)

    Nhs = [np.int32(n) for n in h_g.shape]
    
    kernel_name = "convolve%sd_buf"%(len(data_g.shape)) 
    prog.run_kernel(kernel_name,data_g.shape[::-1],None,
                    data_g.data,h_g.data,res_g.data,
                    *Nhs)

    return res_g

def gpu_mean(data, Nx=10,Ny=10):
    """Function to convolve an imgage with a mean filter on GPU."""
    # create numpy arrays
    
    
    data_g = OCLArray.from_array(data.astype(float32)) 
       
    res_g = OCLArray.empty(data.shape,float32) 
    
    prog = OCLProgram("./OpenCL/gpu_kernels/gpu_mean.cl")
    
    # start kernel on gput
    prog.run_kernel("mean",   # the name of the kernel in the cl file
                    data_g.shape[::-1], # global size, the number of threads e.g. (128,128,) 
                    None,   # local size, just leave it to None
                    data_g.data,res_g.data,
                    int32(Nx),int32(Ny)) 
                    
                    
    return res_g.get()

def perlin2(size, units=None, repeat=(10.,)*2, scale=None, shift=(0, 0)):
    """
        2d perlin noise
        either scale =(10.,10.) or units (5.,5.) have to be given....

        scale is the characteristic length in pixels
    Parameters
    ----------
    size:

    units
    repeat
    scale
    shift

    Returns
    -------

    """

    if scale:
        if np.isscalar(scale):
            scale = (scale,)*2
        repeat = scale
        units = (1.,)*2

    wx, wy = repeat
    dx, dy = units
    offset_x, offset_y = shift

    prog = OCLProgram(abspath("kernels/perlin.cl"))

    d = OCLArray.empty(size[::-1], np.float32)
    prog.run_kernel("perlin2d", d.shape[::-1], None,
                    d.data,
                    np.float32(dx), np.float32(dy),
                    np.float32(wx), np.float32(wy),
                    np.float32(offset_x), np.float32(offset_y),
                    )

    return d.get()

def affine(data, mat = np.identity(4), mode ="linear"):
    """affine transform data with matrix mat

    """ 

    bop = {"linear":"","nearest":"-D USENEAREST"}

    if not mode in bop.keys():
        raise KeyError("mode = '%s' not defined ,valid: %s"%(mode, bop.keys()))
    
    d_im = OCLImage.from_array(data)
    res_g = OCLArray.empty(data.shape,np.float32)
    mat_g = OCLArray.from_array(np.linalg.inv(mat).astype(np.float32,copy=False))

    prog = OCLProgram(abspath("kernels/transformations.cl")
                      , build_options=[bop[mode]])

    prog.run_kernel("affine",
                    data.shape[::-1],None,
                    d_im,res_g.data,mat_g.data)

    return res_g.get()