Archive: Refraction (or something that looks like refraction :)


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);

29th February 2004 01:33 UTC

lol at first I didn't realize you could look around with the mouse and I wasn't impressed at all but after I figured that I have to say that makes it a lot better. The refraction doesn't look quite right to me but then again maybe I'm wrong. Nice work though!


29th February 2004 01:58 UTC

i'm 100% sure that it's correct for a refraction index of 1 *g*
but as is said the formula is guessed it looks like refraction to me and it bends the rays when they hit the surface in the correct direction but since it produces some garbage for refraction indixes < 1, I don't think my formula is correct... but it looks nice though doesn't it? :)
I hope someone with better maths skills than me can solve that problem :)


29th February 2004 12:44 UTC

Refraction goes like this: if ***952;1 is the angle between the incoming ray and the surface normal, and ***952;2 is the angle between the outgoing/refracted ray and the surface normal, then:


sin ***952;2 = n1
sin ***952;1 n2


Where n1 and n2 are the refraction indices of the first and second material.

29th February 2004 13:08 UTC

I found that formula, too, but how do i rotate a vector around an axis that is not the x, y or z axis?

[edit]
That axis should be the normal vector of the plane containing the incoming ray and the sphere surface normal vector right? :)
[/edit]


29th February 2004 16:11 UTC

Well, I guess it works pretty close to correct. I don't really know what you wanted.


29th February 2004 19:03 UTC

http://www.siggraph.org/education/ma...ran/3drota.htm

Check you rotation around an arbitrary axis. That should help.