lamino-filters/src/lamino_bp_generic.cl

// all kernels must process volme voxelwise
// please be careful with local and global workers
// we need a 3D processing (not 2D)

// using tests we show max volume grosse: 1024 x 1024 x 256 or 2048 x 2048 x 64
// theory: 223 MB (CL_DEVICE_MAX_MEM_ALLOC_SIZE) x 5 cards / 4 Bytes (Float) = 1024 x 1024 x 278 voxels

#include <lamino-filter-def.h>

kernel void
lamino_bp_generic (global float *proj,
				   global float *volume,
				   global CLParameters *param)
{
    const ushort vX = get_global_id(0);
    const ushort vY = get_global_id(1);
    const ushort vZ = get_global_id(2);

    /* const ushort vSX = get_global_size(0); */
    /* const ushort vSY = get_global_size(1); */
    /* const ushort vSZ = get_global_size(2); */

    /* const long int idx = (vZ * vSY * vSX) + (vY * vSX) + vX; */

    const long int idx = (vZ * get_global_size(1) * get_global_size(0)) +
                         (vY * get_global_size(0)) + vX;
    //const int idx = (vY * vSX) + vX;

    const float newz = (float)vZ - param->vol_oz;

    /* float newz_matr02 = newz * param->mat_2 + param->proj_ox; */
    /* float newz_matr12 = newz * param->mat_5 + param->proj_oy; */

    /// prepare y info
    const float newy = (float)vY - param->vol_oy;

    /* float newy_matr01 = newy * param->mat_1 + newz_matr02; */
    /* float newy_matr11 = newy * param->mat_4 + newz_matr12; */
    /* float newy_matr01 = newy * param->mat_1 + newz * param->mat_2 + param->proj_ox; */
    /* float newy_matr11 = newy * param->mat_4 + newz * param->mat_5 + param->proj_oy; */

    /// prepare x info
    const float newx = (float)vX - param->vol_ox;

    /* float oldx = newx*param->mat_0 + newy_matr01; */
    const float oldx = newx * param->mat_0 + newy * param->mat_1 + newz * param->mat_2 + param->proj_ox;

    /* float oldy = newx * param->mat_3 + newy_matr11; */
    const float oldy = newx * param->mat_3 + newy * param->mat_4 + newz * param->mat_5 + param->proj_oy;

    /* const float yo = floor(oldy); */
    /* const float xo = floor(oldx); */

    // bilinear interpolation
    /* const float yf_1 = oldy - yo; */
    const float yf_1 = oldy - floor(oldy);
    const float yf_0 = 1.0f - yf_1;

    /* const float xf_1 = oldx - xo; */
    const float xf_1 = oldx - floor(oldx);
    const float xf_0 = 1.0f  - xf_1;

    const int base = ((int) floor(oldx)) + ((int) floor(oldy)) * param->proj_sx;
    /* const int base = (int)xo + (int)yo * param->proj_sx; */
    float result;
// TODO: check that out of proj plain
    result  = proj[base    ] * xf_0 * yf_0;
    result += proj[base + 1] * xf_1 * yf_0;
    result += proj[base + param->proj_sx    ] * xf_0 * yf_1;
    result += proj[base + param->proj_sx + 1] * xf_1 * yf_1;
    volume[idx] += result;
}

__kernel void lamino_clean_vol(__global float *volume)
{
    const ushort vX = get_global_id(0);
    const ushort vY = get_global_id(1);
    const ushort vZ = get_global_id(2);

    const ushort vSX = get_global_size(0);
    const ushort vSY = get_global_size(1);
    const ushort vSZ = get_global_size(2);

    const int idx = (vZ * vSY * vSX) + (vY * vSX) + vX;
    //const int idx = (vY * vSX) + vX;

    volume[idx] = 0;
}

__kernel void lamino_norm_vol(__global float *volume,
				const float factor)
{
    const ushort vX = get_global_id(0);
    const ushort vY = get_global_id(1);
    const ushort vZ = get_global_id(2);

    const ushort vSX = get_global_size(0);
    const ushort vSY = get_global_size(1);
    const ushort vSZ = get_global_size(2);

    const int idx = (vZ * vSY * vSX) + (vY * vSX) + vX;
    //const int idx = (vY * vSX) + vX;

    float val = volume[idx] * factor;
    volume[idx] = val;
}