const sampler_t sampleri = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;

const sampler_t samplerf = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_LINEAR;

#ifndef M_PI

#define M_PI 3.14159265358979323846 // should be defined by the OpenCL compiler acc. to standard

FC(const int row, const int col, const unsigned int filters)

{

return filters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3;

}

kernel void

whitebalance_1ui(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 uint4 pixel = read_imageui(in, sampleri, (int2)(x, y));

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

}

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(rx+y, ry+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;

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

switch(mode)

{

case 1: // DT_IOP_HIGHLIGHTS_LCH

if(pixel.x > clip || pixel.y > clip || pixel.z > clip)

pixel.x = pixel.y = pixel.z = 0.299f * pixel.x + 0.587f * pixel.y + 0.144f*pixel.z;

break;

default: // 0, DT_IOP_HIGHLIGHTS_CLIP

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

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

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

break;

}

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));

float lum[3] = {0.299f, 0.587f, 0.144f};

float accum[3] = {0.0f, 0.0f, 0.0f};

int cnt[3] = {0, 0, 0};

switch(mode)

{

case 1: // DT_IOP_HIGHLIGHTS_LCH

if(pixel.x > clip)

{

// go through all 9 neighbours

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

{

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

{

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

if(px > clip)

{

int c = FC(ry+y+jj, rx+x+ii, filters);

accum[c] += lum[c]*px;

cnt[c] ++;

}

}

}

if(cnt[0] && cnt[1] && cnt[2])

{

pixel.x = 0.0f;

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

pixel.x += accum[c]/cnt[c];

}

else pixel.x = clip;

}

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 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(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);

}

XYZ_to_Lab(float *xyz, float *lab)

{

xyz[0] *= (1.0f/0.9642f);

xyz[2] *= (1.0f/0.8242f);

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

xyz[c] = xyz[c] > 0.008856f ? native_powr(xyz[c], 1.0f/3.0f) : 7.787f*xyz[c] + 16.0f/116.0f;

lab[0] = 116.0f * xyz[1] - 16.0f;

lab[1] = 500.0f * (xyz[0] - xyz[1]);

lab[2] = 200.0f * (xyz[1] - xyz[2]);

}

kernel void

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

global float *mat, read_only image2d_t lutr, read_only image2d_t lutg, read_only image2d_t lutb,

const int map_blues,

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], Lab[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);

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

if(map_blues && 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] += mat[3*j+i] * cam[i];

}

XYZ_to_Lab(XYZ, (float *)&pixel);

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, 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));

const float L_in = pixel.x/100.0f;

// use lut or extrapolation:

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

if (autoscale_ab == 0)

{

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

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

pixel.y = lookup(table_a, a_in);

pixel.z = lookup(table_b, b_in);

}

else if(pixel.x > 0.01f)

{

pixel.y *= L/pixel.x;

pixel.z *= L/pixel.x;

}

else

{

pixel.y *= L/0.01f;

pixel.z *= L/0.01f;

}

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);

}

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*t)*native_sin(M_PI*t/width)/(M_PI*M_PI*t*t));

return r;

}

sinf_fast(float t)

{

const float a = 4/(M_PI*M_PI);

const float p = 0.225f;

t = a*t*(M_PI - 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*r)*sinf_fast(M_PI*t/width))/(DT_LANCZOS_EPSILON + M_PI*M_PI*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 int2 roi_out, const float scale_in, const float scale_out,

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

{

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 + scale_out * ci.x + x + 0.5f;

pi.y = roi_out.y + scale_out * ci.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;

po.x -= roi_in.x;

po.y -= roi_in.y;

const int ii = (int)po.x;

const int jj = (int)po.y;

float4 o;

if (ii >=0 && jj >= 0 && ii <= in_width-2 && jj <= in_height-2)

o = read_imagef(in, samplerf, po);

else

o = (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 int2 roi_out, const float scale_in, const float scale_out,

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

{

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 + scale_out * ci.x + x + 0.5f;

pi.y = roi_out.y + scale_out * ci.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;

po.x -= roi_in.x;

po.y -= roi_in.y;

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++)

{

float wx = interpolation_func_bicubic((float)(tx + ii) - po.x);

float wy = interpolation_func_bicubic((float)(ty + jj) - po.y);

float w = wx * wy;

pixel += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)) * w;

weight += w;

}

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

}

__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 int2 roi_out, const float scale_in, const float scale_out,

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

{

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 + scale_out * ci.x + x + 0.5f;

pi.y = roi_out.y + scale_out * ci.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;

po.x -= roi_in.x;

po.y -= roi_in.y;

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++)

{

float wx = interpolation_func_lanczos(2, (float)(tx + ii) - po.x);

float wy = interpolation_func_lanczos(2, (float)(ty + jj) - po.y);

float w = wx * wy;

pixel += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)) * w;

weight += w;

}

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

}

__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 int2 roi_out, const float scale_in, const float scale_out,

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

{

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 + scale_out * ci.x + x + 0.5f;

pi.y = roi_out.y + scale_out * ci.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;

po.x -= roi_in.x;

po.y -= roi_in.y;

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++)

{

float wx = interpolation_func_lanczos(3, (float)(tx + ii) - po.x);

float wy = interpolation_func_lanczos(3, (float)(ty + jj) - po.y);

float w = wx * wy;

pixel += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)) * w;

weight += w;

}

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

}

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 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);

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 : 0.0f;

rx = ppi[2] - roi_in_x;

ry = ppi[3] - roi_in_y;

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

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

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 : 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 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);

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++)

{

float wx = interpolation_func_bicubic((float)(tx + ii) - rx);

float wy = interpolation_func_bicubic((float)(ty + jj) - ry);

float w = wx * wy;

sum += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)).x * w;

weight += w;

}

pixel.x = (tx >= 0 && ty >= 0 && tx <= iwidth - 1 && ty <= iheight - 1) ? sum/weight : 0.0f;

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++)

{

float wx = interpolation_func_bicubic((float)(tx + ii) - rx);

float wy = interpolation_func_bicubic((float)(ty + jj) - ry);

float w = wx * wy;

sum2 += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)).yw * w;

weight += w;

}

pixel.yw = (tx >= 0 && ty >= 0 && tx <= iwidth - 1 && ty <= iheight - 1) ? sum2/weight : (float2)0.0f;

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

ry = ppi[5] - (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++)

{

float wx = interpolation_func_bicubic((float)(tx + ii) - rx);

float wy = interpolation_func_bicubic((float)(ty + jj) - ry);

float w = wx * wy;

sum += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)).z * w;

weight += w;

}

pixel.z = (tx >= 0 && ty >= 0 && tx <= iwidth - 1 && ty <= iheight - 1) ? sum/weight : 0.0f;

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

}

kernel void

lens_distort_lanczos2 (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 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);

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++)

{

float wx = interpolation_func_lanczos(2, (float)(tx + ii) - rx);

float wy = interpolation_func_lanczos(2, (float)(ty + jj) - ry);

float w = wx * wy;

sum += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)).x * w;

weight += w;

}

pixel.x = (tx >= 0 && ty >= 0 && tx <= iwidth - 1 && ty <= iheight - 1) ? sum/weight : 0.0f;

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++)

{

float wx = interpolation_func_lanczos(2, (float)(tx + ii) - rx);

float wy = interpolation_func_lanczos(2, (float)(ty + jj) - ry);

float w = wx * wy;

sum2 += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)).yw * w;

weight += w;

}

pixel.yw = (tx >= 0 && ty >= 0 && tx <= iwidth - 1 && ty <= iheight - 1) ? sum2/weight : (float2)0.0f;

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

ry = ppi[5] - (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++)

{

float wx = interpolation_func_lanczos(2, (float)(tx + ii) - rx);

float wy = interpolation_func_lanczos(2, (float)(ty + jj) - ry);

float w = wx * wy;

sum += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)).z * w;

weight += w;

}

pixel.z = (tx >= 0 && ty >= 0 && tx <= iwidth - 1 && ty <= iheight - 1) ? sum/weight : 0.0f;

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

}

kernel void

lens_distort_lanczos3 (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 x = get_global_id(0);

const int y = get_global_id(1);

const int kwidth = 3;

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);

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++)

{

float wx = interpolation_func_lanczos(3, (float)(tx + ii) - rx);

float wy = interpolation_func_lanczos(3, (float)(ty + jj) - ry);

float w = wx * wy;

sum += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)).x * w;

weight += w;

}

pixel.x = (tx >= 0 && ty >= 0 && tx <= iwidth - 1 && ty <= iheight - 1) ? sum/weight : 0.0f;

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++)

{

float wx = interpolation_func_lanczos(3, (float)(tx + ii) - rx);

float wy = interpolation_func_lanczos(3, (float)(ty + jj) - ry);

float w = wx * wy;

sum2 += read_imagef(in, sampleri, (int2)(tx + ii, ty + jj)).yw * w;

weight += w;

}

pixel.yw = (tx >= 0 && ty >= 0 && tx <= iwidth - 1 && ty <= iheight - 1) ? sum2/weight : (float2)0.0f;

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

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

tx = rx;

ty = ry;

sum = 0.0f;

weight = 0.0f;