123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139 |
- // 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>
- #define PI 3.1415926535897932384626433832795028841971693993751
- #define PI_2 1.5707963267948966192313216916397514420985846996876
- float
- interpolate (global float *proj,
- global CLParameters *param,
- float oldx,
- float oldy)
- {
- /* Check if vertically out of projection plane boundaries */
- if (oldy < 0)
- oldy = 0.5;;
- if (oldy > param->proj_sy)
- oldy = param->proj_sy - 0.5;
- /* bilinear interpolation */
- const float yf_1 = oldy - floor(oldy);
- const float yf_0 = 1.0f - yf_1;
- 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;
- float result;
- 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;
- return result;
- }
- 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 long int idx = (vZ * get_global_size(1) * get_global_size(0)) +
- (vY * get_global_size(0)) + vX;
- const float newz = (float) vZ * param->z_spacing - param->vol_oz;
- const float newy = (float) vY - param->vol_oy;
- const float newx = (float) vX - param->vol_ox;
- const float oldx = newx * param->mat_0 + newy * param->mat_1 + newz * param->mat_2 + param->proj_ox;
- const float oldy = newx * param->mat_3 + newy * param->mat_4 + newz * param->mat_5 + param->proj_oy;
- volume[idx] += interpolate (proj, param, oldx, oldy);
- }
- kernel void
- lamino_bp_varied (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 long int idx = (vZ * get_global_size(1) * get_global_size(0)) +
- (vY * get_global_size(0)) + vX;
- const float newz = (float) vZ * param->z_spacing - param->vol_oz;
- const float newy = (float) vY - param->vol_oy;
- const float newx = (float) vX - param->vol_ox;
- /* Vary alpha depending on the z height */
- const float alpha = -3 * PI_2 + param->theta;
- const float cf = cos(param->phi);
- const float ct = cos(alpha);
- const float cg = cos(param->psi);
- const float sf = sin(param->phi);
- const float st = sin(alpha);
- const float sg = sin(param->psi);
- const float mat_0 = cg * cf - sg * st * sf;
- const float mat_1 = -cg * sf - sg * st * cf;
- const float mat_2 = -sg * ct;
- const float mat_3 = sg * cf + cg * st * sf;
- const float mat_4 = -sg * sf + cg * st * cf;
- const float mat_5 = cg * ct;
- /* Force newz values to be equal to 0 */
- const float oldx = newx * mat_0 + newy * mat_1 /* + newz * mat_2 */ + param->proj_ox + newz / param->vol_sz * param->proj_ox_variation;
- const float oldy = newx * mat_3 + newy * mat_4 /* + newz * mat_5 */ + param->proj_oy;
- volume[idx] += interpolate (proj, param, oldx, oldy);
- }
- 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;
- 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;
- }
|