OpenCL image distortion

blackhole.opencl

// See http://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/sampler_t.html
const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_FILTER_NEAREST | CLK_ADDRESS_CLAMP_TO_EDGE;

constant int n = 2; // field exponent
constant float e = 1e-3; // smoothing factor
constant float c = 5e-4; // field multiplier
constant float br = 4e-5; // brightness factor

float z(float r, float k) {
	return k / (e + pow(r, n - 1));
}

float dz_dr(float r, float k) {
	//return k * -(n - 1) * pow(r, -n);
	float rn = pow(r, n);
	return k * -(n - 1) * rn / pow(e * r + rn, 2);
}

float2 field(float x, float y, float k) {
	float r = sqrt(x*x + y*y);
	float d = dz_dr(r, k);
	return (float2)(d * x / r, d * y / r);
}

float height(float x, float y, float k) {
	float r = sqrt(x*x + y*y);
	return dz_dr(r, k);
}


float2 ripple(float x, float y, float a, float w) {
	float m = sqrt(x*x + y*y) / 100;

	float s = sin(w*m);
	float x1 = a*y*m*s;
	float y1 = -a*x*m*s;

	return (float2)(x1, y1);
}

float2 displacement(float x, float y, float width) {
	/*return ripple(x - 0.5, aspect*(y - 0.5 + 0.05), 350, 100)
             - ripple(x - 0.5, aspect*(y - 0.5 - 0.05), 350, 100)
             + ripple(x - 0.5, aspect*(y - 0.5), 60, 1200);*/

	float k = c / width;
	return width * (150*field(x - 0.5, y - 0.5, k));
}

__kernel void test(
	read_only image2d_t inputImage,
	write_only image2d_t outputImage)
{
   // See http://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/get_image_dim.html
   int2 dimensions = get_image_dim(inputImage);
   int w = dimensions.x, h = dimensions.y;
   
   // See http://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/get_image_channel_data_type.html
   int channelDataType = get_image_channel_data_type(inputImage);
   
   // See http://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/get_image_channel_order.html
   int channelOrder = get_image_channel_order(inputImage);
   
   int xi = get_global_id(0), yi = get_global_id(1);

   float aspect = ((float) h) / ((float) w);
   float x = ((float) xi) / ((float) w), y = aspect*((float) yi) / ((float) h);
   int2 coordinates = (int2)(xi, yi);

   float2 disp = displacement(x, y, (float) w);
   int2 fromCoordinates = coordinates + convert_int2(disp);

   // See http://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/read_imagef2d.html
   //float4 pixel = read_imagef(inputImage, sampler, coordinates);
   float4 pixel = read_imagef(inputImage, sampler, fromCoordinates);
   float4 transformedPixel = pixel;
   float b = dot(disp, normalize((float2)(1, 0.664)));
   transformedPixel.x += b*br;
   transformedPixel.y += b*br;
   transformedPixel.z += b*br;
   //float4 transformedPixel = (float4)(0);//pixel;
   
   // See http://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/write_image.html
   write_imagef(outputImage, coordinates, transformedPixel);
}