nuclear@2: struct RendInfo {
nuclear@2: 	int xsz, ysz;
nuclear@2: 	int num_sph;
nuclear@2: 	int max_iter;
nuclear@2: };
nuclear@2: 
nuclear@2: struct Sphere {
nuclear@2: 	float4 pos;
nuclear@2: 	float radius;
nuclear@2: 	float4 color;
nuclear@2: };
nuclear@2: 
nuclear@2: struct Ray {
nuclear@2: 	float4 origin, dir;
nuclear@2: };
nuclear@2: 
nuclear@2: struct SurfPoint {
nuclear@2: 	float t;
nuclear@2: 	float3 pos, norm;
nuclear@2: };
nuclear@2: 
nuclear@2: #define EPSILON 1e-6
nuclear@2: 
nuclear@2: bool intersect(struct Ray ray, __global const struct Sphere *sph, struct SurfPoint *sp);
nuclear@2: 
nuclear@2: __kernel void render(__global float4 *fb,
nuclear@2: 		__global const struct RendInfo *rinf,
nuclear@2: 		__global const struct Sphere *sphlist,
nuclear@2: 		__global const struct Ray *primrays)
nuclear@2: {
nuclear@2: 	int idx = get_global_id(0);
nuclear@2: 
nuclear@2: 	struct Ray ray = primrays[idx];
nuclear@2: 	struct SurfPoint sp;
nuclear@2: 
nuclear@2: 	if(intersect(ray, sphlist, &sp)) {
nuclear@2: 		fb[idx] = (float4)(1, 0, 0, 1);
nuclear@2: 	} else {
nuclear@2: 		fb[idx] = (float4)(0, 0, 0, 1);
nuclear@2: 	}
nuclear@2: }
nuclear@2: 
nuclear@2: bool intersect(struct Ray ray,
nuclear@2: 		__global const struct Sphere *sph,
nuclear@2: 		struct SurfPoint *sp)
nuclear@2: {
nuclear@2: 	float3 dir = ray.dir.xyz;
nuclear@2: 	float3 orig = ray.origin.xyz;
nuclear@2: 	float3 spos = sph->pos.xyz;
nuclear@2: 
nuclear@2: 	float a = dot(dir, dir);
nuclear@2: 	float b = 2.0 * dir.x * (orig.x - spos.x) +
nuclear@2: 		2.0 * dir.y * (orig.y - spos.y) +
nuclear@2: 		2.0 * dir.z * (orig.z - spos.z);
nuclear@2: 	float c = dot(spos, spos) + dot(orig, orig) +
nuclear@2: 		2.0 * dot(-spos, orig) - sph->radius * sph->radius;
nuclear@2: 
nuclear@2: 	float d = b * b - 4.0 * a * c;
nuclear@2: 	if(d < 0.0) return false;
nuclear@2: 
nuclear@2: 	float sqrt_d = sqrt(d);
nuclear@2: 	float t1 = (-b + sqrt_d) / (2.0 * a);
nuclear@2: 	float t2 = (-b - sqrt_d) / (2.0 * a);
nuclear@2: 
nuclear@2: 	if(t1 < EPSILON) t1 = t2;
nuclear@2: 	if(t2 < EPSILON) t2 = t1;
nuclear@2: 	float t = t1 < t2 ? t1 : t2;
nuclear@2: 
nuclear@2: 	if(t < EPSILON || t > 1.0) {
nuclear@2: 		return false;
nuclear@2: 	}
nuclear@2: 
nuclear@2: 	sp->t = t;
nuclear@2: 	sp->pos = orig + dir * sp->t;
nuclear@2: 	sp->norm = (sp->pos - spos) / sph->radius;
nuclear@2: 	return true;
nuclear@2: }