1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/rt.cl	Mon Jul 12 10:38:07 2010 +0300
     1.3 @@ -0,0 +1,77 @@
     1.4 +struct RendInfo {
     1.5 +	int xsz, ysz;
     1.6 +	int num_sph;
     1.7 +	int max_iter;
     1.8 +};
     1.9 +
    1.10 +struct Sphere {
    1.11 +	float4 pos;
    1.12 +	float radius;
    1.13 +	float4 color;
    1.14 +};
    1.15 +
    1.16 +struct Ray {
    1.17 +	float4 origin, dir;
    1.18 +};
    1.19 +
    1.20 +struct SurfPoint {
    1.21 +	float t;
    1.22 +	float3 pos, norm;
    1.23 +};
    1.24 +
    1.25 +#define EPSILON 1e-6
    1.26 +
    1.27 +bool intersect(struct Ray ray, __global const struct Sphere *sph, struct SurfPoint *sp);
    1.28 +
    1.29 +__kernel void render(__global float4 *fb,
    1.30 +		__global const struct RendInfo *rinf,
    1.31 +		__global const struct Sphere *sphlist,
    1.32 +		__global const struct Ray *primrays)
    1.33 +{
    1.34 +	int idx = get_global_id(0);
    1.35 +
    1.36 +	struct Ray ray = primrays[idx];
    1.37 +	struct SurfPoint sp;
    1.38 +
    1.39 +	if(intersect(ray, sphlist, &sp)) {
    1.40 +		fb[idx] = (float4)(1, 0, 0, 1);
    1.41 +	} else {
    1.42 +		fb[idx] = (float4)(0, 0, 0, 1);
    1.43 +	}
    1.44 +}
    1.45 +
    1.46 +bool intersect(struct Ray ray,
    1.47 +		__global const struct Sphere *sph,
    1.48 +		struct SurfPoint *sp)
    1.49 +{
    1.50 +	float3 dir = ray.dir.xyz;
    1.51 +	float3 orig = ray.origin.xyz;
    1.52 +	float3 spos = sph->pos.xyz;
    1.53 +
    1.54 +	float a = dot(dir, dir);
    1.55 +	float b = 2.0 * dir.x * (orig.x - spos.x) +
    1.56 +		2.0 * dir.y * (orig.y - spos.y) +
    1.57 +		2.0 * dir.z * (orig.z - spos.z);
    1.58 +	float c = dot(spos, spos) + dot(orig, orig) +
    1.59 +		2.0 * dot(-spos, orig) - sph->radius * sph->radius;
    1.60 +
    1.61 +	float d = b * b - 4.0 * a * c;
    1.62 +	if(d < 0.0) return false;
    1.63 +
    1.64 +	float sqrt_d = sqrt(d);
    1.65 +	float t1 = (-b + sqrt_d) / (2.0 * a);
    1.66 +	float t2 = (-b - sqrt_d) / (2.0 * a);
    1.67 +
    1.68 +	if(t1 < EPSILON) t1 = t2;
    1.69 +	if(t2 < EPSILON) t2 = t1;
    1.70 +	float t = t1 < t2 ? t1 : t2;
    1.71 +
    1.72 +	if(t < EPSILON || t > 1.0) {
    1.73 +		return false;
    1.74 +	}
    1.75 +
    1.76 +	sp->t = t;
    1.77 +	sp->pos = orig + dir * sp->t;
    1.78 +	sp->norm = (sp->pos - spos) / sph->radius;
    1.79 +	return true;
    1.80 +}