28th February 2004 22:07 UTC
Refraction (or something that looks like refraction :)
Hi there :) I experimented a bit with refraction and managed to write an AVS preset that displays a sphere refracting the image on a plane behind it.
I don't know whether my refraction maths are correct.
(or better: I *guessed* them *g*)
Maybe one of you math geeks can check that and/or give me a better formula :)
Here's the code of the DM:
//I didn't find a suitable translation for the
//German word "Auftreffpunkt" which means the
//point where something hits something.
//I translated it with "hit point".
//1st hit point is the point where the ray hits
//the outer bounds of the sphere.
//2nd hit point is the point where the ray hits
//the inner bounds of the sphere after being refracted.
//--- Calculate projection plane coordinates. ---
dx=x/zoom; dy=y/zoom; dz=1;
//--- Transform them. ---
dx1=dx*cz-dy*sz; dy1=dx*sz+dy*cz; dy=dy1*cx-sx; dz1=dy1*sx+dz*cx; dx=dx1*cy-dz1*sy; dz=dx1*sy+dz1*cy;
//--- Trace from camera to sphere surface ---
// r² = x² + y² + z²
// r² = (dx*k+ox-mx)² + (dy*k+oy-my)² + (dz*k+oz-mz)²
kdist=(sqr(dx) + sqr(dy) + sqr(dz));
kp = (2*dx*(ox-mx) + 2*dy*(oy-my) + 2*dz*(oz-mz))/kdist;
kq = (sqr(ox-mx)+sqr(oy-my)+sqr(oz-mz)-sqr(rad))/kdist;
kdet=sqr(kp/2)-kq;
k1=-kp/2+sqrt(kdet);
k2=-kp/2-sqrt(kdet);
k1=if(below(k1,0),INFIN,k1);
k2=if(below(k2,0),INFIN,k2);
ksphere=min(k1,k2);
//--- Trace from camera to plane ---
pa=0; pb=0; pc=1; pd=-20;
k1=(-pa*ox-pb*oy-pc*oz-pd)/(pa*dx + pb*dy + pc*dz);
kplane=if(below(k1,0),INFIN,k1);
//vec a1: o -> 1st hit point //vec ik1: <0> -> 1st hit point
a1x=dx*ksphere; ik1x=a1x+ox;
a1y=dy*ksphere; ik1y=a1y+oy;
a1z=dz*ksphere; ik1z=a1z+oz;
//--- Does the ray hit the sphere? ---
xdist=ik1x-mx; ydist=ik1y-my; zdist=ik1z-mz;
distkugel=sqrt(sqr(xdist)+sqr(ydist)+sqr(zdist));
isplane=above(abs(9-distkugel),epsilon);
k=if(isplane,kplane,ksphere);
ix1=dx*kplane+ox; iy1=dy*kplane+oy; iz1=dz*kplane+oz;
//The refraction method i used is a linear interpolation
//between the normal vector and the traced ray
//I'm not sure whether this method is correct. If you know a
//better method, let me know :)
//--- 1st Refraction ---
//Coefficients for recfraction
pr=1/refindex; ps=1-pr;
//normalise vec a1
a1r=invsqrt(sqr(a1x)+sqr(a1y)+sqr(a1z));
a1x=a1x*a1r; a1y=a1y*a1r; a1z=a1z*a1r;
//vec n: sphere surface normal at the 1st hit point
nx=ik1x-mx;
ny=ik1y-my;
nz=ik1z-mz;
//normalise vec n
nr=invsqrt(sqr(nx)+sqr(ny)+sqr(nz));
nx=nx*nr; ny=ny*nr; nz=nz*nr;
//vec b1: Ray after 1st refraction
b1x=pr*a1x-ps*nx;
b1y=pr*a1y-ps*ny;
b1z=pr*a1z-ps*nz;
//normalise vec b1
b1r=invsqrt(sqr(b1x) + sqr(b1y) + sqr(b1z));
b1x=b1x*b1r; b1y=b1y*b1r; b1z=b1z*b1r;
//--- Trace inside the sphere---
ka=sqr(b1x) + sqr(b1y) + sqr(b1z);
kp=(2*b1x*(ik1x-mx) + 2*b1y*(ik1y-my) + 2*b1z*(ik1z-mz))/ka;
kq=(sqr(ik1x-mx)+sqr(ik1y-my)+sqr(ik1z-mz)-sqr(rad))/ka;
kdet=sqr(kp/2)-kq;
k1=-kp/2+sqrt(kdet);
//k2=-kp/2-sqrt(kdet);
k1=if(below(k1,0),INFIN,k1);
//k2=if(below(k2,0),INFIN,k2);
ksphere2=k1;
//--- 2nd refraction ---
//Coefficients for recfraction
pr=1/refindex; ps=1-pr;
//vec a2: 1st hit point -> 2nd hit point // vec ik2: <0> -> 2nd hit point
a2x=b1x*ksphere2; ik2x=a2x+ik1x;
a2y=b1y*ksphere2; ik2y=a2y+ik1y;
a2z=b1z*ksphere2; ik2z=a2z+ik1z;
//normalise vec a2 (keep a2r for light dimming calculations)
a2r=sqrt(sqr(a2x)+sqr(a2y)+sqr(a2z));
a2x=a2x/a2r; a2y=a2y/a2r; a2z=a2z/a2r;
//vec n: sphere surface normal at the 2nd hit point
nx=ik2x-mx;
ny=ik2y-my;
nz=ik2z-mz;
//normalise vec n
nr=invsqrt(sqr(nx)+sqr(ny)+sqr(nz));
nx=nx*nr; ny=ny*nr; nz=nz*nr;
//vec b2: Ray after 2nd refraction
b2x=pr*a2x+ps*nx;
b2y=pr*a2y+ps*ny;
b2z=pr*a2z+ps*nz;
//--- Trace refracted ray to the plane ---
k1=(-pa*ik2x-pb*ik2y-pc*ik2z-pd)/(pa*b2x + pb*b2y + pc*b2z);
kplane2=if(below(k1,0),INFIN,k1);
ix2=b2x*kplane2+ik2x; iy2=b2y*kplane2+ik2y; iz2=b2z*kplane2+ik2z;
//--- Get coordinates of the point where the ray ends
ix=if(isplane,ix1,ix2);
iy=if(isplane,iy1,iy2);
iz=if(isplane,iz1,iz2);
//--- Determine texture coordinates ---
x=ix/40; y=iy/40;
//--- Determine Alpha ---
// Calculate distance to camera.
xdist=ix-ox; ydist=iy-oy; zdist=iz-oz;
dist=sqrt(sqr(xdist)+sqr(ydist)+sqr(zdist));
//if (plane) then (100% Alpha) else (calculate from distance between the 2 hitpoints)
alpha=if(isplane,1,exp(-a2r/50));
//if (dist>=viewdist) then (0% Alpha)
alpha=alpha*below(dist,viewdist);
alpha=if(equal(k,INFIN),0,alpha);