#include "common.h"

#include "colorspace.cl"

BL(const int row, const int col)

{

return (((row & 1) << 1) + (col & 1));

}

kernel void

rawprepare_1f(read_only image2d_t in, write_only image2d_t out,

const int width, const int height,

const int cx, const int cy,

global const float *sub, global const float *div,

const int rx, const int ry)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x < cx || y < cy) return;

if(x >= width + cx || y >= height + cy) return;

const float pixel = read_imageui(in, sampleri, (int2)(x, y)).x;

const int id = BL(ry+cy+y, rx+cx+x);

const float pixel_scaled = max(0.0f, (pixel - sub[id])) / div[id];

write_imagef(out, (int2)(x-cx, y-cy), (float4)(pixel_scaled, 0.0f, 0.0f, 0.0f));

}

kernel void

rawprepare_4f(read_only image2d_t in, write_only image2d_t out,

const int width, const int height,

const int cx, const int cy,

global const float *black, global const float *div)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x < cx || y < cy) return;

if(x >= width + cx || y >= height + cy) return;

float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

pixel.xyz = (pixel.xyz - black[0]) / div[0];

pixel.xyz = max(0.0f, pixel.xyz);

write_imagef(out, (int2)(x-cx, y-cy), pixel);

}

kernel void

invert_1f(read_only image2d_t in, write_only image2d_t out, const int width, const int height, global float *color,

const unsigned int filters, const int rx, const int ry)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

const float pixel = read_imagef(in, sampleri, (int2)(x, y)).x;

const float inv_pixel = color[FC(ry+y, rx+x, filters)] - pixel;

write_imagef (out, (int2)(x, y), (float4)(clamp(inv_pixel, 0.0f, 1.0f), 0.0f, 0.0f, 0.0f));

}

kernel void

invert_4f(read_only image2d_t in, write_only image2d_t out, const int width, const int height, global float *color,

const unsigned int filters, const int rx, const int ry)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

pixel.x = color[0] - pixel.x;

pixel.y = color[1] - pixel.y;

pixel.z = color[2] - pixel.z;

pixel.xyz = clamp(pixel.xyz, 0.0f, 1.0f);

write_imagef (out, (int2)(x, y), pixel);

}

kernel void

whitebalance_1f(read_only image2d_t in, write_only image2d_t out, const int width, const int height, global float *coeffs,

const unsigned int filters, const int rx, const int ry)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

const float pixel = read_imagef(in, sampleri, (int2)(x, y)).x;

write_imagef (out, (int2)(x, y), (float4)(pixel * coeffs[FC(ry+y, rx+x, filters)], 0.0f, 0.0f, 0.0f));

}

kernel void

whitebalance_4f(read_only image2d_t in, write_only image2d_t out, const int width, const int height, global float *coeffs,

const unsigned int filters, const int rx, const int ry)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

const float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

write_imagef (out, (int2)(x, y), (float4)(pixel.x * coeffs[0], pixel.y * coeffs[1], pixel.z * coeffs[2], pixel.w));

}

kernel void

exposure (read_only image2d_t in, write_only image2d_t out, const int width, const int height, const float black, const float scale)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

pixel.xyz = (pixel.xyz - black)*scale;

write_imagef (out, (int2)(x, y), pixel);

}

kernel void

highlights_4f (read_only image2d_t in, write_only image2d_t out, const int width, const int height,

const int mode, const float clip)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

// 4f/pixel means that this has been debayered already.

// it's thus hopeless to recover highlights here (this code path is just used for preview and non-raw images)

float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

// default: // 0, DT_IOP_HIGHLIGHTS_CLIP

pixel.x = fmin(clip, pixel.x);

pixel.y = fmin(clip, pixel.y);

pixel.z = fmin(clip, pixel.z);

write_imagef (out, (int2)(x, y), pixel);

}

kernel void

highlights_1f (read_only image2d_t in, write_only image2d_t out, const int width, const int height,

const int mode, const float clip, const int rx, const int ry, const int filters)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

// just carry over other 3 channels:

float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

const float near_clip = 0.96f*clip;

const float post_clip = 1.10f*clip;

float mean = 0.0f;

float blend = 0.0f;

switch(mode)

{

case 1: // DT_IOP_HIGHLIGHTS_LCH

{

// go through a bayer block of 2x2

for(int jj=0;jj<=1;jj++)

{

for(int ii=0;ii<=1;ii++)

{

float px = read_imagef(in, sampleri, (int2)(x+ii, y+jj)).x;

mean += px*.25f;

blend += (fmin(post_clip, px) - near_clip)/(post_clip-near_clip);

}

}

blend = clamp(blend, 0.0f, 1.0f);

if(blend > 0.0f)

pixel.x = blend*mean + (1.0f-blend)*pixel.x;

break;

}

default: // 0, DT_IOP_HIGHLIGHTS_CLIP

pixel.x = fmin(clip, pixel.x);

break;

}

write_imagef (out, (int2)(x, y), pixel);

}

lookup_unbounded(read_only image2d_t lut, const float x, global const float *a)

{

// in case the tone curve is marked as linear, return the fast

// path to linear unbounded (does not clip x at 1)

if(a[0] >= 0.0f)

{

if(x < 1.0f/a[0])

{

const int xi = clamp(x*65535.0f, 0.0f, 65535.0f);

const int2 p = (int2)((xi & 0xff), (xi >> 8));

return read_imagef(lut, sampleri, p).x;

}

else return a[1] * native_powr(x*a[0], a[2]);

}

else return x;

}

lookup_unbounded_twosided(read_only image2d_t lut, const float x, global const float *a)

{

// in case the tone curve is marked as linear, return the fast

// path to linear unbounded (does not clip x at 1)

if(a[0] >= 0.0f)

{

const float ar = 1.0f/a[0];

const float al = 1.0f - 1.0f/a[3];

if(x < ar && x >= al)

{

// lut lookup

const int xi = clamp(x*65535.0f, 0.0f, 65535.0f);

const int2 p = (int2)((xi & 0xff), (xi >> 8));

return read_imagef(lut, sampleri, p).x;

}

else

{

// two-sided extrapolation (with inverted x-axis for left side)

const float xx = (x >= ar) ? x : 1.0f - x;

global const float *aa = (x >= ar) ? a : a + 3;

return aa[1] * native_powr(xx*aa[0], aa[2]);

}

}

else return x;

}

lookup(read_only image2d_t lut, const float x)

{

int xi = clamp(x*65535.0f, 0.0f, 65535.0f);

int2 p = (int2)((xi & 0xff), (xi >> 8));

return read_imagef(lut, sampleri, p).x;

}

kernel void

basecurve (read_only image2d_t in, write_only image2d_t out, const int width, const int height,

read_only image2d_t table, global float *a)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

// use lut or extrapolation:

pixel.x = lookup_unbounded(table, pixel.x, a);

pixel.y = lookup_unbounded(table, pixel.y, a);

pixel.z = lookup_unbounded(table, pixel.z, a);

write_imagef (out, (int2)(x, y), pixel);

}

kernel void

colorin_unbound (read_only image2d_t in, write_only image2d_t out, const int width, const int height,

global float *cmat, global float *lmat,

read_only image2d_t lutr, read_only image2d_t lutg, read_only image2d_t lutb,

const int blue_mapping, global float *a)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

float cam[3], XYZ[3];

cam[0] = lookup_unbounded(lutr, pixel.x, a);

cam[1] = lookup_unbounded(lutg, pixel.y, a+3);

cam[2] = lookup_unbounded(lutb, pixel.z, a+6);

if(blue_mapping)

{

const float YY = cam[0] + cam[1] + cam[2];

if(YY > 0.0f)

{

// manual gamut mapping. these values cause trouble when converting back from Lab to sRGB:

const float zz = cam[2] / YY;

// lower amount and higher bound_z make the effect smaller.

// the effect is weakened the darker input values are, saturating at bound_Y

const float bound_z = 0.5f, bound_Y = 0.8f;

const float amount = 0.11f;

if (zz > bound_z)

{

const float t = (zz - bound_z) / (1.0f - bound_z) * fmin(1.0f, YY / bound_Y);

cam[1] += t * amount;

cam[2] -= t * amount;

}

}

}

// now convert camera to XYZ using the color matrix

for(int j=0;j<3;j++)

{

XYZ[j] = 0.0f;

for(int i=0;i<3;i++) XYZ[j] += cmat[3*j+i] * cam[i];

}

float4 xyz = (float4)(XYZ[0], XYZ[1], XYZ[2], 0.0f);

pixel.xyz = XYZ_to_Lab(xyz).xyz;

write_imagef (out, (int2)(x, y), pixel);

}

kernel void

colorin_clipping (read_only image2d_t in, write_only image2d_t out, const int width, const int height,

global float *cmat, global float *lmat,

read_only image2d_t lutr, read_only image2d_t lutg, read_only image2d_t lutb,

const int blue_mapping, global float *a)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

float cam[3], RGB[3], XYZ[3];

cam[0] = lookup_unbounded(lutr, pixel.x, a);

cam[1] = lookup_unbounded(lutg, pixel.y, a+3);

cam[2] = lookup_unbounded(lutb, pixel.z, a+6);

if(blue_mapping)

{

const float YY = cam[0] + cam[1] + cam[2];

if(YY > 0.0f)

{

// manual gamut mapping. these values cause trouble when converting back from Lab to sRGB:

const float zz = cam[2] / YY;

// lower amount and higher bound_z make the effect smaller.

// the effect is weakened the darker input values are, saturating at bound_Y

const float bound_z = 0.5f, bound_Y = 0.8f;

const float amount = 0.11f;

if (zz > bound_z)

{

const float t = (zz - bound_z) / (1.0f - bound_z) * fmin(1.0f, YY / bound_Y);

cam[1] += t * amount;

cam[2] -= t * amount;

}

}

}

// convert camera to RGB using the first color matrix

for(int j=0;j<3;j++)

{

RGB[j] = 0.0f;

for(int i=0;i<3;i++) RGB[j] += cmat[3*j+i] * cam[i];

}

// clamp at this stage

for(int i=0; i<3; i++) RGB[i] = clamp(RGB[i], 0.0f, 1.0f);

// convert clipped RGB to XYZ

for(int j=0;j<3;j++)

{

XYZ[j] = 0.0f;

for(int i=0;i<3;i++) XYZ[j] += lmat[3*j+i] * RGB[i];

}

float4 xyz = (float4)(XYZ[0], XYZ[1], XYZ[2], 0.0f);

pixel.xyz = XYZ_to_Lab(xyz).xyz;

write_imagef (out, (int2)(x, y), pixel);

}

kernel void

tonecurve (read_only image2d_t in, write_only image2d_t out, const int width, const int height,

read_only image2d_t table_L, read_only image2d_t table_a, read_only image2d_t table_b,

const int autoscale_ab, const int unbound_ab, global float *coeffs_L, global float *coeffs_ab,

const float low_approximation)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

const float L_in = pixel.x/100.0f;

// use lut or extrapolation:

const float L = lookup_unbounded(table_L, L_in, coeffs_L);

if (autoscale_ab == 0)

{

const float a_in = (pixel.y + 128.0f) / 256.0f;

const float b_in = (pixel.z + 128.0f) / 256.0f;

if (unbound_ab == 0)

{

pixel.y = lookup(table_a, a_in);

pixel.z = lookup(table_b, b_in);

}

else

{

// use lut or two-sided extrapolation

pixel.y = lookup_unbounded_twosided(table_a, a_in, coeffs_ab);

pixel.z = lookup_unbounded_twosided(table_b, b_in, coeffs_ab + 6);

}

}

else

{

if(pixel.x > 0.01f)

{

pixel.y *= L/pixel.x;

pixel.z *= L/pixel.x;

}

else

{

pixel.y *= low_approximation;

pixel.z *= low_approximation;

}

}

pixel.x = L;

write_imagef (out, (int2)(x, y), pixel);

}

__kernel void

colorcorrection (read_only image2d_t in, write_only image2d_t out, const int width, const int height,

const float saturation, const float a_scale, const float a_base,

const float b_scale, const float b_base)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

float4 pixel = read_imagef(in, sampleri, (int2)(x, y));

pixel.y = saturation*(pixel.y + pixel.x * a_scale + a_base);

pixel.z = saturation*(pixel.z + pixel.x * b_scale + b_base);

write_imagef (out, (int2)(x, y), pixel);

}

mul_mat_vec_2(const float4 m, const float2 *p, float2 *o)

{

(*o).x = (*p).x*m.x + (*p).y*m.y;

(*o).y = (*p).x*m.z + (*p).y*m.w;

}

backtransform(float2 *p, float2 *o, const float4 m, const float2 t)

{

(*p).y /= (1.0f + (*p).x*t.x);

(*p).x /= (1.0f + (*p).y*t.y);

mul_mat_vec_2(m, p, o);

}

keystone_backtransform(float2 *i, const float4 k_space, const float2 ka, const float4 ma, const float2 mb)

{

float xx = (*i).x - k_space.x;

float yy = (*i).y - k_space.y;

/*float u = ka.x-kb.x+kc.x-kd.x;

float div = ((ma.z*xx-ma.x*yy)*mb.y+(ma.y*yy-ma.w*xx)*mb.x+ma.x*ma.w-ma.y*ma.z);

(*i).x = (ma.w*xx-ma.y*yy)/div + ka.x;

(*i).y =-(ma.z*xx-ma.x*yy)/div + ka.y;

}

interpolation_func_bicubic(float t)

{

float r;

t = fabs(t);

r = (t >= 2.0f) ? 0.0f : ((t > 1.0f) ? (0.5f*(t*(-t*t + 5.0f*t - 8.0f) + 4.0f)) : (0.5f*(t*(3.0f*t*t - 5.0f*t) + 2.0f)));

return r;

}

#define DT_LANCZOS_EPSILON (1e-9f)

#if 0

interpolation_func_lanczos(float width, float t)

{

float ta = fabs(t);

float r = (ta > width) ? 0.0f : ((ta < DT_LANCZOS_EPSILON) ? 1.0f : width*native_sin(M_PI_F*t)*native_sin(M_PI_F*t/width)/(M_PI_F*M_PI_F*t*t));

return r;

}

sinf_fast(float t)

{

const float a = 4.0f/(M_PI_F*M_PI_F);

const float p = 0.225f;

t = a*t*(M_PI_F - fabs(t));

return p*(t*fabs(t) - t) + t;

}

interpolation_func_lanczos(float width, float t)

{

/* Compute a value for sinf(pi.t) in [-pi pi] for which the value will be

int a = (int)t;

float r = t - (float)a;

// Compute the correct sign for sinf(pi.r)

union { float f; unsigned int i; } sign;

sign.i = ((a&1)<<31) | 0x3f800000;

return (DT_LANCZOS_EPSILON + width*sign.f*sinf_fast(M_PI_F*r)*sinf_fast(M_PI_F*t/width))/(DT_LANCZOS_EPSILON + M_PI_F*M_PI_F*t*t);

}

__kernel void

clip_rotate_bilinear(read_only image2d_t in, write_only image2d_t out, const int width, const int height,

const int in_width, const int in_height,

const int2 roi_in, const float2 roi_out, const float scale_in, const float scale_out,

const int flip, const float2 t, const float2 k, const float4 mat,

const float4 k_space, const float2 ka, const float4 ma, const float2 mb)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

float2 pi, po;

pi.x = roi_out.x + x + 0.5f;

pi.y = roi_out.y + y + 0.5f;

pi.x -= flip ? t.y * scale_out : t.x * scale_out;

pi.y -= flip ? t.x * scale_out : t.y * scale_out;

pi /= scale_out;

backtransform(&pi, &po, mat, k);

po *= scale_in;

po.x += t.x * scale_in;

po.y += t.y * scale_in;

if (k_space.z > 0.0f) keystone_backtransform(&po,k_space,ka,ma,mb);

po.x -= roi_in.x + 0.5f;

po.y -= roi_in.y + 0.5f;

const int ii = (int)po.x;

const int jj = (int)po.y;

float4 o = (ii >=0 && jj >= 0 && ii <= in_width-2 && jj <= in_height-2) ? read_imagef(in, samplerf, po) : (float4)0.0f;

write_imagef (out, (int2)(x, y), o);

}

__kernel void

clip_rotate_bicubic(read_only image2d_t in, write_only image2d_t out, const int width, const int height,

const int in_width, const int in_height,

const int2 roi_in, const float2 roi_out, const float scale_in, const float scale_out,

const int flip, const float2 t, const float2 k, const float4 mat,

const float4 k_space, const float2 ka, const float4 ma, const float2 mb)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

const int kwidth = 2;

if(x >= width || y >= height) return;

float2 pi, po;

pi.x = roi_out.x + x + 0.5f;

pi.y = roi_out.y + y + 0.5f;

pi.x -= flip ? t.y * scale_out : t.x * scale_out;

pi.y -= flip ? t.x * scale_out : t.y * scale_out;

pi /= scale_out;

backtransform(&pi, &po, mat, k);

po *= scale_in;

po.x += t.x * scale_in;

po.y += t.y * scale_in;

if (k_space.z > 0.0f) keystone_backtransform(&po,k_space,ka,ma,mb);

po.x -= roi_in.x + 0.5f;

po.y -= roi_in.y + 0.5f;

int tx = po.x;

int ty = po.y;

float4 pixel = (float4)0.0f;

float weight = 0.0f;

for(int jj = 1 - kwidth; jj <= kwidth; jj++)

for(int ii= 1 - kwidth; ii <= kwidth; ii++)

{

const int i = tx + ii;

const int j = ty + jj;

float wx = interpolation_func_bicubic((float)i - po.x);

float wy = interpolation_func_bicubic((float)j - po.y);

float w = wx * wy;

pixel += read_imagef(in, samplerc, (int2)(i, j)) * w;

weight += w;

}

pixel = weight > 0.0f ? pixel / weight : (float4)0.0f;

write_imagef (out, (int2)(x, y), pixel);

}

__kernel void

clip_rotate_lanczos2(read_only image2d_t in, write_only image2d_t out, const int width, const int height,

const int in_width, const int in_height,

const int2 roi_in, const float2 roi_out, const float scale_in, const float scale_out,

const int flip, const float2 t, const float2 k, const float4 mat,

const float4 k_space, const float2 ka, const float4 ma, const float2 mb)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

const int kwidth = 2;

if(x >= width || y >= height) return;

float2 pi, po;

pi.x = roi_out.x + x + 0.5f;

pi.y = roi_out.y + y + 0.5f;

pi.x -= flip ? t.y * scale_out : t.x * scale_out;

pi.y -= flip ? t.x * scale_out : t.y * scale_out;

pi /= scale_out;

backtransform(&pi, &po, mat, k);

po *= scale_in;

po.x += t.x * scale_in;

po.y += t.y * scale_in;

if (k_space.z > 0.0f) keystone_backtransform(&po,k_space,ka,ma,mb);

po.x -= roi_in.x + 0.5f;

po.y -= roi_in.y + 0.5f;

int tx = po.x;

int ty = po.y;

float4 pixel = (float4)0.0f;

float weight = 0.0f;

for(int jj = 1 - kwidth; jj <= kwidth; jj++)

for(int ii= 1 - kwidth; ii <= kwidth; ii++)

{

const int i = tx + ii;

const int j = ty + jj;

float wx = interpolation_func_lanczos(2, (float)i - po.x);

float wy = interpolation_func_lanczos(2, (float)j - po.y);

float w = wx * wy;

pixel += read_imagef(in, samplerc, (int2)(i, j)) * w;

weight += w;

}

pixel = weight > 0.0f ? pixel / weight : (float4)0.0f;

write_imagef (out, (int2)(x, y), pixel);

}

__kernel void

clip_rotate_lanczos3(read_only image2d_t in, write_only image2d_t out, const int width, const int height,

const int in_width, const int in_height,

const int2 roi_in, const float2 roi_out, const float scale_in, const float scale_out,

const int flip, const float2 t, const float2 k, const float4 mat,

const float4 k_space, const float2 ka, const float4 ma, const float2 mb)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

const int kwidth = 3;

if(x >= width || y >= height) return;

float2 pi, po;

pi.x = roi_out.x + x + 0.5f;

pi.y = roi_out.y + y + 0.5f;

pi.x -= flip ? t.y * scale_out : t.x * scale_out;

pi.y -= flip ? t.x * scale_out : t.y * scale_out;

pi /= scale_out;

backtransform(&pi, &po, mat, k);

po *= scale_in;

po.x += t.x * scale_in;

po.y += t.y * scale_in;

if (k_space.z > 0.0f) keystone_backtransform(&po,k_space,ka,ma,mb);

po.x -= roi_in.x + 0.5f;

po.y -= roi_in.y + 0.5f;

int tx = (int)po.x;

int ty = (int)po.y;

float4 pixel = (float4)0.0f;

float weight = 0.0f;

for(int jj = 1 - kwidth; jj <= kwidth; jj++)

for(int ii= 1 - kwidth; ii <= kwidth; ii++)

{

const int i = tx + ii;

const int j = ty + jj;

float wx = interpolation_func_lanczos(3, (float)i - po.x);

float wy = interpolation_func_lanczos(3, (float)j - po.y);

float w = wx * wy;

pixel += read_imagef(in, samplerc, (int2)(i, j)) * w;

weight += w;

}

pixel = weight > 0.0f ? pixel / weight : (float4)0.0f;

write_imagef (out, (int2)(x, y), pixel);

}

kernel void

lens_distort_bilinear (read_only image2d_t in, write_only image2d_t out, const int width, const int height,

const int iwidth, const int iheight, const int roi_in_x, const int roi_in_y, global float *pi,

const int do_nan_checks)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

if(x >= width || y >= height) return;

float4 pixel;

float rx, ry;

const int piwidth = 2*3*width;

global float *ppi = pi + mad24(y, piwidth, 2*3*x);

if(do_nan_checks)

{

bool valid = true;

for(int i = 0; i < 6; i++) valid = valid && isfinite(ppi[i]);

if(!valid)

{

pixel = (float4)0.0f;

write_imagef (out, (int2)(x, y), pixel);

return;

}

}

rx = ppi[0] - roi_in_x;

ry = ppi[1] - roi_in_y;

pixel.x = (rx >= 0 && ry >= 0 && rx <= iwidth - 1 && ry <= iheight - 1) ? read_imagef(in, samplerf, (float2)(rx, ry)).x : NAN;

rx = ppi[2] - roi_in_x;

ry = ppi[3] - roi_in_y;

pixel.yw = (rx >= 0 && ry >= 0 && rx <= iwidth - 1 && ry <= iheight - 1) ? read_imagef(in, samplerf, (float2)(rx, ry)).yw : (float2)NAN;

rx = ppi[4] - roi_in_x;

ry = ppi[5] - roi_in_y;

pixel.z = (rx >= 0 && ry >= 0 && rx <= iwidth - 1 && ry <= iheight - 1) ? read_imagef(in, samplerf, (float2)(rx, ry)).z : NAN;

pixel = all(isfinite(pixel.xyz)) ? pixel : (float4)0.0f;

write_imagef (out, (int2)(x, y), pixel);

}

kernel void

lens_distort_bicubic (read_only image2d_t in, write_only image2d_t out, const int width, const int height,

const int iwidth, const int iheight, const int roi_in_x, const int roi_in_y, global float *pi,

const int do_nan_checks)

{

const int x = get_global_id(0);

const int y = get_global_id(1);

const int kwidth = 2;

if(x >= width || y >= height) return;

float4 pixel = (float4)0.0f;

float rx, ry;

int tx, ty;

float sum, weight;

float2 sum2;

const int piwidth = 2*3*width;

global float *ppi = pi + mad24(y, piwidth, 2*3*x);

if(do_nan_checks)

{

bool valid = true;

for(int i = 0; i < 6; i++) valid = valid && isfinite(ppi[i]);

if(!valid)

{

pixel = (float4)0.0f;

write_imagef (out, (int2)(x, y), pixel);

return;

}

}

rx = ppi[0] - (float)roi_in_x;

ry = ppi[1] - (float)roi_in_y;

tx = rx;

ty = ry;

sum = 0.0f;

weight = 0.0f;

for(int jj = 1 - kwidth; jj <= kwidth; jj++)

for(int ii= 1 - kwidth; ii <= kwidth; ii++)

{

const int i = tx + ii;

const int j = ty + jj;

float wx = interpolation_func_bicubic((float)i - rx);

float wy = interpolation_func_bicubic((float)j - ry);

float w = (i < 0 || j < 0 || i >= iwidth || j >= iheight) ? 0.0f : wx * wy;

sum += read_imagef(in, samplerc, (int2)(i, j)).x * w;

weight += w;

}

pixel.x = sum/weight;

rx = ppi[2] - (float)roi_in_x;

ry = ppi[3] - (float)roi_in_y;

tx = rx;

ty = ry;

sum2 = (float2)0.0f;

weight = 0.0f;

for(int jj = 1 - kwidth; jj <= kwidth; jj++)

for(int ii= 1 - kwidth; ii <= kwidth; ii++)

{

const int i = tx + ii;

const int j = ty + jj;

float wx = interpolation_func_bicubic((float)i - rx);