nuclear@0: uniform sampler3D volume;
nuclear@0: uniform sampler2D ray_tex;
nuclear@4: uniform sampler1D xfer_tex;
nuclear@7: uniform float ray_step;
nuclear@7: uniform float zclip;
nuclear@0: 
nuclear@0: struct Ray {
nuclear@0: 	vec3 origin, dir;
nuclear@0: };
nuclear@0: 
nuclear@4: struct AABBox {
nuclear@4: 	vec3 min, max;
nuclear@4: };
nuclear@4: 
nuclear@0: struct ISect {
nuclear@0: 	bool hit;
nuclear@7: 	float t0, t1;
nuclear@0: 	vec3 pos;
nuclear@0: 	vec3 normal;
nuclear@0: };
nuclear@0: 
nuclear@1: vec3 sky(Ray ray);
nuclear@7: vec3 ray_march(Ray ray, float t0, float t1);
nuclear@26: vec3 shade(Ray ray, vec3 pos, vec3 norm, vec3 ldir);
nuclear@0: Ray get_primary_ray();
nuclear@5: ISect intersect_aabb(Ray ray, AABBox aabb);
nuclear@25: ISect intersect_sphere(Ray ray, float rad);
nuclear@25: 
nuclear@25: //#define USE_AABB
nuclear@0: 
nuclear@0: void main()
nuclear@0: {
nuclear@0: 	Ray ray = get_primary_ray();
nuclear@0: 
nuclear@4: 	vec3 color = vec3(0.0, 0.0, 0.0);
nuclear@4: 
nuclear@25: #ifdef USE_AABB
nuclear@25: 	const AABBox aabb = AABBox(vec3(-1.0, -1.0, -1.0), vec3(1.0, 1.0, 1.0));
nuclear@5: 	ISect res = intersect_aabb(ray, aabb);
nuclear@5: 	if(res.hit) {
nuclear@25: 		color += ray_march(ray, res.t0, res.t1);
nuclear@4: 	}
nuclear@25: #else
nuclear@25: 	ISect res = intersect_sphere(ray, 1.5);
nuclear@25: 	if(res.hit) {
nuclear@25: 		color += ray_march(ray, res.t0, res.t1);
nuclear@25: 	}
nuclear@25: #endif
nuclear@4: 
nuclear@4: 	gl_FragColor = vec4(color, 1.0);
nuclear@0: }
nuclear@0: 
nuclear@25: vec3 sky(Ray ray)
nuclear@25: {
nuclear@25: 	vec3 c0 = vec3(0.2, 0.5, 0.8);
nuclear@25: 	vec3 c1 = vec3(0.8, 0.5, 0.2);
nuclear@25: 	float t = smoothstep(-0.2, 0.2, ray.dir.z);
nuclear@25: 	return mix(c0, c1, t);
nuclear@25: }
nuclear@25: 
nuclear@20: float eval(vec3 pos, out vec3 grad)
nuclear@4: {
nuclear@5: 	vec3 tc = pos * 0.5 + 0.5;
nuclear@7: 
nuclear@7: 	if(tc.x < 0.0 || tc.y < 0.0 || tc.z < zclip || tc.x > 1.0 || tc.y > 1.0 || tc.z > 1.0) {
nuclear@20: 		grad = vec3(0.0, 0.0, 0.0);
nuclear@7: 		return 0.0;
nuclear@7: 	}
nuclear@7: 
nuclear@20: 	vec4 texel = texture3D(volume, tc);
nuclear@20: 	grad = texel.xyz;
nuclear@20: 
nuclear@20: 	return texture1D(xfer_tex, texel.a).x;
nuclear@4: }
nuclear@4: 
nuclear@7: vec3 ray_march(Ray ray, float t0, float t1)
nuclear@0: {
nuclear@0: 	float energy = 1.0;
nuclear@7: 	float t = t0;
nuclear@7: 	vec3 col = vec3(0.0, 0.0, 0.0);
nuclear@33: 	vec3 ldir = normalize(vec3(-1, 1, -4));
nuclear@0: 
nuclear@7: 
nuclear@7: 	while(t < t1) {
nuclear@7: 		vec3 pos = ray.origin + ray.dir * t;
nuclear@7: 		t += ray_step;
nuclear@7: 
nuclear@7: 		vec3 norm;
nuclear@20: 		float val = eval(pos, norm);
nuclear@7: 
nuclear@25: 		float energy_drop = exp(val * -ray_step);	// * scatter_coeff ?
nuclear@25: 		energy *= energy_drop;
nuclear@25: 
nuclear@26: 		vec3 irrad = shade(ray, pos, normalize(norm), ldir);
nuclear@25: 
nuclear@25: 		col += (1.0 - energy_drop) * energy * irrad;
nuclear@1: 		if(energy < 0.001) {
nuclear@1: 			break;
nuclear@1: 		}
nuclear@0: 	}
nuclear@0: 
nuclear@7: 	return col;
nuclear@0: }
nuclear@0: 
nuclear@26: vec3 shade(Ray ray, vec3 pos, vec3 norm, vec3 ldir)
nuclear@0: {
nuclear@26: 	vec3 vdir = -normalize(ray.dir);
nuclear@26: 	vec3 hdir = normalize(ldir + vdir);
nuclear@0: 
nuclear@25: 	float ndotl = abs(dot(ldir, norm));
nuclear@25: 	float ndoth = abs(dot(hdir, norm));
nuclear@0: 
nuclear@8: 	vec3 dcol = vec3(0.9, 0.9, 0.9) * max(ndotl, 0.0);
nuclear@8: 	vec3 scol = vec3(0.5, 0.5, 0.5) * pow(max(ndoth, 0.0), 50.0);
nuclear@0: 
nuclear@25: 	return dcol + scol;
nuclear@0: }
nuclear@0: 
nuclear@0: Ray get_primary_ray()
nuclear@0: {
nuclear@0: 	Ray ray;
nuclear@0: 	vec2 tc = gl_TexCoord[0].xy;
nuclear@0: 	ray.dir = gl_NormalMatrix * normalize(texture2D(ray_tex, tc).xyz);
nuclear@0: 	ray.origin = (gl_ModelViewMatrix * vec4(0.0, 0.0, 0.0, 1.0)).xyz;
nuclear@0: 	return ray;
nuclear@0: }
nuclear@4: 
nuclear@35: #ifdef USE_AABB
nuclear@5: ISect intersect_aabb(Ray ray, AABBox aabb)
nuclear@4: {
nuclear@5: 	ISect res;
nuclear@5: 	res.hit = false;
nuclear@7: 	res.t0 = res.t1 = 0.0;
nuclear@5: 
nuclear@4: 	if(ray.origin.x >= aabb.min.x && ray.origin.y >= aabb.min.y && ray.origin.z >= aabb.min.z &&
nuclear@4: 			ray.origin.x < aabb.max.x && ray.origin.y < aabb.max.y && ray.origin.z < aabb.max.z) {
nuclear@5: 		res.hit = true;
nuclear@5: 		return res;
nuclear@4: 	}
nuclear@4: 
nuclear@4: 	vec4 bbox[2];
nuclear@4: 	bbox[0] = vec4(aabb.min.x, aabb.min.y, aabb.min.z, 0);
nuclear@4: 	bbox[1] = vec4(aabb.max.x, aabb.max.y, aabb.max.z, 0);
nuclear@4: 
nuclear@4: 	int xsign = int(ray.dir.x < 0.0);
nuclear@4: 	float invdirx = 1.0 / ray.dir.x;
nuclear@4: 	float tmin = (bbox[xsign].x - ray.origin.x) * invdirx;
nuclear@4: 	float tmax = (bbox[1 - xsign].x - ray.origin.x) * invdirx;
nuclear@4: 
nuclear@4: 	int ysign = int(ray.dir.y < 0.0);
nuclear@4: 	float invdiry = 1.0 / ray.dir.y;
nuclear@4: 	float tymin = (bbox[ysign].y - ray.origin.y) * invdiry;
nuclear@4: 	float tymax = (bbox[1 - ysign].y - ray.origin.y) * invdiry;
nuclear@4: 
nuclear@4: 	if(tmin > tymax || tymin > tmax) {
nuclear@5: 		return res;
nuclear@4: 	}
nuclear@4: 
nuclear@4: 	if(tymin > tmin) tmin = tymin;
nuclear@4: 	if(tymax < tmax) tmax = tymax;
nuclear@4: 
nuclear@4: 	int zsign = int(ray.dir.z < 0.0);
nuclear@4: 	float invdirz = 1.0 / ray.dir.z;
nuclear@4: 	float tzmin = (bbox[zsign].z - ray.origin.z) * invdirz;
nuclear@4: 	float tzmax = (bbox[1 - zsign].z - ray.origin.z) * invdirz;
nuclear@4: 
nuclear@4: 	if(tmin > tzmax || tzmin > tmax) {
nuclear@5: 		return res;
nuclear@4: 	}
nuclear@4: 
nuclear@7: 	res.t0 = tmin;
nuclear@7: 	res.t1 = tmax;
nuclear@5: 	res.hit = true;
nuclear@5: 	return res;
nuclear@4: }
nuclear@35: #else
nuclear@25: 
nuclear@25: ISect intersect_sphere(Ray ray, float rad)
nuclear@25: {
nuclear@25: 	ISect res;
nuclear@25: 	res.hit = false;
nuclear@25: 	res.t0 = res.t1 = 0.0;
nuclear@25: 
nuclear@25: 	float a = dot(ray.dir, ray.dir);
nuclear@25: 	float b = 2.0 * dot(ray.dir, ray.origin);
nuclear@25: 	float c = dot(ray.origin, ray.origin) - rad * rad;
nuclear@25: 	float d = b * b - 4.0 * a * c;
nuclear@25: 
nuclear@25: 	if(d < 0.0) {
nuclear@25: 		return res;
nuclear@25: 	}
nuclear@25: 
nuclear@25: 	float sqrt_d = sqrt(d);
nuclear@25: 	float t0 = (-b + sqrt_d) / (2.0 * a);
nuclear@25: 	float t1 = (-b - sqrt_d) / (2.0 * a);
nuclear@25: 
nuclear@35: 	res.t0 = max(min(t0, t1), 0.0);
nuclear@35: 	res.t1 = max(max(t0, t1), 0.0);
nuclear@25: 	res.hit = true;
nuclear@25: 
nuclear@25: 	/*res.pos = ray.origin + ray.dir * res.t0;
nuclear@25: 	res.normal = normalize(res.pos);*/
nuclear@25: 	return res;
nuclear@25: }
nuclear@35: #endif