1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/sdr/julia.p.glsl	Thu Oct 21 23:39:26 2010 +0300
     1.3 @@ -0,0 +1,289 @@
     1.4 +/* vim: set ft=glsl:ts=4:sw=4 */
     1.5 +uniform vec4 seed;
     1.6 +uniform sampler2D ray_tex;
     1.7 +uniform float err_thres;
     1.8 +uniform int iter;
     1.9 +
    1.10 +#define quat(s, x, y, z)	vec4(x, y, z, s)
    1.11 +#define quat_identity()		vec4(0.0, 0.0, 0.0, 1.0)
    1.12 +
    1.13 +#define vec2quat(v)		(v).wxyz
    1.14 +
    1.15 +struct Ray {
    1.16 +	vec3 origin;
    1.17 +	vec3 dir;
    1.18 +};
    1.19 +
    1.20 +struct Julia {
    1.21 +	bool inside;
    1.22 +	vec4 q;
    1.23 +	vec4 qprime;
    1.24 +};
    1.25 +
    1.26 +struct ISect {
    1.27 +	bool hit;
    1.28 +	float t;
    1.29 +	vec3 pos;
    1.30 +	vec3 normal;
    1.31 +	vec3 color;
    1.32 +	float kr;
    1.33 +};
    1.34 +
    1.35 +ISect find_intersection(Ray ray);
    1.36 +vec3 shade(Ray ray, ISect isect);
    1.37 +float amboc(ISect isect);
    1.38 +vec3 sky(Ray ray);
    1.39 +Julia julia(vec4 q, vec4 c);
    1.40 +float julia_dist(vec4 z);
    1.41 +vec3 julia_grad(vec4 z);
    1.42 +vec4 quat_mul(vec4 q1, vec4 q2);
    1.43 +vec4 quat_sq(vec4 q);
    1.44 +float quat_length_sq(vec4 q);
    1.45 +ISect ray_julia(Ray ray);
    1.46 +ISect ray_sphere(Ray ray, float rad);
    1.47 +ISect ray_floor(Ray ray);
    1.48 +Ray get_primary_ray();
    1.49 +
    1.50 +
    1.51 +void main()
    1.52 +{
    1.53 +	Ray ray = get_primary_ray();
    1.54 +
    1.55 +	float energy = 1.0;
    1.56 +	vec3 color = vec3(0.0, 0.0, 0.0);
    1.57 +
    1.58 +	while(energy > 0.001) {
    1.59 +		ISect res = find_intersection(ray);
    1.60 +
    1.61 +		if(res.hit) {
    1.62 +			color += shade(ray, res) * energy;
    1.63 +			energy *= res.kr;
    1.64 +
    1.65 +			ray.origin = res.pos;
    1.66 +			ray.dir = reflect(ray.dir, res.normal);
    1.67 +		} else {
    1.68 +			color += sky(ray) * energy;
    1.69 +			break;
    1.70 +		}
    1.71 +	}
    1.72 +
    1.73 +	gl_FragColor = vec4(color, 1.0);
    1.74 +}
    1.75 +
    1.76 +
    1.77 +ISect find_intersection(Ray ray)
    1.78 +{
    1.79 +	ISect res;
    1.80 +	res.hit = false;
    1.81 +	
    1.82 +	ISect bhit = ray_sphere(ray, 2.0);
    1.83 +	if(bhit.hit) {
    1.84 +		ray.origin = bhit.pos;
    1.85 +		res = ray_julia(ray);
    1.86 +	}
    1.87 +
    1.88 +	if(!res.hit) {
    1.89 +		res = ray_floor(ray);
    1.90 +	}
    1.91 +	return res;
    1.92 +}
    1.93 +
    1.94 +vec3 shade(Ray ray, ISect isect)
    1.95 +{
    1.96 +	vec3 ldir = normalize(vec3(10.0, 10.0, -10.0) - isect.pos);
    1.97 +	vec3 vdir = -ray.dir;
    1.98 +	vec3 hdir = normalize(ldir + vdir);
    1.99 +
   1.100 +	float ndotl = dot(ldir, isect.normal);
   1.101 +	float ndoth = dot(hdir, isect.normal);
   1.102 +
   1.103 +	vec3 dcol = /*isect.color * max(ndotl, 0.0) */ amboc(isect);
   1.104 +	vec3 scol = vec3(1.0, 1.0, 1.0) * pow(max(ndoth, 0.0), 40.0);
   1.105 +
   1.106 +	return /*vec3(0.05, 0.05, 0.05) + */dcol;// + scol;
   1.107 +}
   1.108 +
   1.109 +#define AO_STEP		0.04
   1.110 +#define AO_MAGIC	8.0
   1.111 +float amboc(ISect isect)
   1.112 +{
   1.113 +	float sum = 0.0;
   1.114 +
   1.115 +	for(float fi=0.0; fi<5.0; fi+=1.0) {
   1.116 +		float sample_dist = fi * AO_STEP;
   1.117 +		vec3 pt = isect.pos + isect.normal * sample_dist;
   1.118 +		float jdist = julia_dist(quat(pt.x, pt.y, pt.z, 0.0));
   1.119 +
   1.120 +		sum += 1.0 / pow(2.0, fi) * (sample_dist - jdist);
   1.121 +	}
   1.122 +	
   1.123 +	return 1.0 - AO_MAGIC * sum;
   1.124 +}
   1.125 +
   1.126 +vec3 sky(Ray ray)
   1.127 +{
   1.128 +	vec3 col1 = vec3(0.75, 0.78, 0.8);
   1.129 +	vec3 col2 = vec3(0.56, 0.7, 1.0);
   1.130 +
   1.131 +	float t = max(ray.dir.y, -0.5);
   1.132 +	return mix(col1, col2, t);
   1.133 +}
   1.134 +
   1.135 +Julia julia(vec4 q, vec4 c)
   1.136 +{
   1.137 +	Julia res;
   1.138 +	res.inside = true;
   1.139 +
   1.140 +	res.q = q;
   1.141 +	res.qprime = quat_identity();
   1.142 +
   1.143 +	for(int i=0; i<iter; i++) {
   1.144 +		res.qprime = 2.0 * quat_mul(res.q, res.qprime);
   1.145 +		res.q = quat_sq(res.q) + c;
   1.146 +
   1.147 +		if(dot(res.q, res.q) > 8.0) {
   1.148 +			res.inside = false;
   1.149 +			break;
   1.150 +		}
   1.151 +	}
   1.152 +	return res;
   1.153 +}
   1.154 +
   1.155 +float julia_dist(vec4 z)
   1.156 +{
   1.157 +	Julia jres = julia(z, seed);
   1.158 +
   1.159 +	float lenq = length(jres.q);
   1.160 +	float lenqprime = length(jres.qprime);
   1.161 +
   1.162 +	return 0.5 * lenq * log(lenq) / lenqprime;
   1.163 +}
   1.164 +
   1.165 +#define OFFS 1e-4
   1.166 +vec3 julia_grad(vec4 z)
   1.167 +{
   1.168 +	vec3 grad;
   1.169 +	grad.x = julia_dist(z + quat(OFFS, 0.0, 0.0, 0.0)) - julia_dist(z - quat(OFFS, 0.0, 0.0, 0.0));
   1.170 +	grad.y = julia_dist(z + quat(0.0, OFFS, 0.0, 0.0)) - julia_dist(z - quat(0.0, OFFS, 0.0, 0.0));
   1.171 +	grad.z = julia_dist(z + quat(0.0, 0.0, OFFS, 0.0)) - julia_dist(z - quat(0.0, 0.0, OFFS, 0.0));
   1.172 +	return grad;
   1.173 +}
   1.174 +
   1.175 +vec4 quat_mul(vec4 q1, vec4 q2)
   1.176 +{
   1.177 +	vec4 res;
   1.178 +	res.w = q1.w * q2.w - dot(q1.xyz, q2.xyz);
   1.179 +	res.xyz = q1.w * q2.xyz + q2.w * q1.xyz + cross(q1.xyz, q2.xyz);
   1.180 +	return res;
   1.181 +}
   1.182 +
   1.183 +vec4 quat_sq(vec4 q)
   1.184 +{
   1.185 +	vec4 res;
   1.186 +	res.w = q.w * q.w - dot(q.xyz, q.xyz);
   1.187 +	res.xyz = 2.0 * q.w * q.xyz;
   1.188 +	return res;
   1.189 +}
   1.190 +
   1.191 +ISect ray_julia(Ray inray)
   1.192 +{
   1.193 +	float dist_acc = 0.0;
   1.194 +	Ray ray = inray;
   1.195 +	ISect res;
   1.196 +
   1.197 +	for(float fi=0.0; ; fi+=0.1) {
   1.198 +		vec4 q = quat(ray.origin.x, ray.origin.y, ray.origin.z, 0.0);
   1.199 +		
   1.200 +		float dist = julia_dist(q);
   1.201 +
   1.202 +		ray.origin += ray.dir * dist;
   1.203 +		dist_acc += dist;
   1.204 +
   1.205 +		if(dist < err_thres) {
   1.206 +			res.hit = true;
   1.207 +			res.t = dist_acc;
   1.208 +			res.pos = ray.origin;
   1.209 +			res.normal = normalize(julia_grad(quat(res.pos.x, res.pos.y, res.pos.z, 0.0)));
   1.210 +			res.color = vec3(0.75, 0.8, 0.9);//abs(res.normal) * 0.2;
   1.211 +			//res.kr = 0.6;
   1.212 +			res.kr = 0.0;
   1.213 +			break;
   1.214 +		}
   1.215 +
   1.216 +		if(dot(ray.origin, ray.origin) > 100.0) {
   1.217 +			res.hit = false;
   1.218 +			break;
   1.219 +		}
   1.220 +	}
   1.221 +
   1.222 +	return res;
   1.223 +}
   1.224 +
   1.225 +ISect ray_sphere(Ray ray, float rad)
   1.226 +{
   1.227 +	ISect res;
   1.228 +	res.hit = false;
   1.229 +
   1.230 +	float a = dot(ray.dir, ray.dir);
   1.231 +	float b = 2.0 * dot(ray.dir, ray.origin);
   1.232 +	float c = dot(ray.origin, ray.origin) - rad * rad;
   1.233 +
   1.234 +	float d = b * b - 4.0 * a * c;
   1.235 +	if(d < 0.0) return res;
   1.236 +
   1.237 +	float sqrt_d = sqrt(d);
   1.238 +	float t1 = (-b + sqrt_d) / (2.0 * a);
   1.239 +	float t2 = (-b - sqrt_d) / (2.0 * a);
   1.240 +
   1.241 +	if((t1 >= 0.0 || t2 >= 0.0)) {
   1.242 +		if(t1 < 0.0) t1 = t2;
   1.243 +		if(t2 < 0.0) t2 = t1;
   1.244 +		
   1.245 +		res.hit = true;
   1.246 +		res.t = min(t1, t2);
   1.247 +		res.pos = ray.origin + ray.dir * res.t;
   1.248 +		res.color = vec3(1.0, 0.3, 0.2);
   1.249 +		res.normal = res.pos / rad;
   1.250 +	}
   1.251 +
   1.252 +	return res;
   1.253 +}
   1.254 +
   1.255 +#define FLOOR_HEIGHT	(-2.0)
   1.256 +
   1.257 +ISect ray_floor(Ray ray)
   1.258 +{
   1.259 +	ISect res;
   1.260 +	res.hit = false;
   1.261 +
   1.262 +	if(ray.origin.y < FLOOR_HEIGHT || ray.dir.y >= 0.0) {
   1.263 +		return res;
   1.264 +	}
   1.265 +
   1.266 +	res.normal = vec3(0.0, 1.0, 0.0);
   1.267 +	float ndotdir = dot(res.normal, ray.dir);
   1.268 +
   1.269 +	float t = (FLOOR_HEIGHT - ray.origin.y) / ndotdir;
   1.270 +	res.pos = ray.origin + ray.dir * t;
   1.271 +
   1.272 +	if(abs(res.pos.x) > 8.0 || abs(res.pos.z) > 8.0) {
   1.273 +		res.hit = false;
   1.274 +	} else {
   1.275 +		res.hit = true;
   1.276 +	
   1.277 +		float chess = mod(floor(res.pos.x) + floor(res.pos.z), 2.0);
   1.278 +		res.color = mix(vec3(0.498, 0.165, 0.149), vec3(0.776, 0.851, 0.847), chess);
   1.279 +		res.kr = 0.0;
   1.280 +	}
   1.281 +	return res;
   1.282 +}
   1.283 +
   1.284 +Ray get_primary_ray()
   1.285 +{
   1.286 +	Ray ray;
   1.287 +	vec2 tc = gl_TexCoord[0].xy;
   1.288 +	ray.dir = gl_NormalMatrix * normalize(texture2D(ray_tex, tc).xyz);
   1.289 +	ray.origin = (gl_ModelViewMatrix * vec4(0.0, 0.0, 0.0, 1.0)).xyz;
   1.290 +	return ray;
   1.291 +}
   1.292 +