Simplifying ray tracing with affine transformation

For our ray tracing lab, we need to trace objects of varying parameters.  For example, we need to know how to trace a sphere of radius 10, centered at (1,1,1),but also how to trace the sphere of radius 1 centered at (0,0,0).  The key concept to take away from the ray tracing methods in our text “An Integrated Introduction..” is that if we know how to trace a single “canonical” sphere, then we know how to trace all other spheres since they are just transformed version of the canonical sphere.  

Without loss of generality, we can take this canonical sphere to be the unit sphere centered at the origin ((Subscript[S, orig])).

Create a sphere just for illustration.  We won’t be using it for intersection.

step=45;
verts=Flatten[Table[{Sin[a]Cos[b],Cos[a]Cos[b],Sin[b]},{a,0,2 Pi,((2Pi)/(step-1))},{b,0,2 Pi,((2 Pi)/(step-1))}],1];
inds=Flatten[Table[{(i-1)+(j-1)*step+1,(i)+(j-1)*step+1,(i)+(j)*step+1,(i-1)+(j)*step+1},{i,step-1},{j,step-1}],1];

Graphics3D[GraphicsComplex[verts,Polygon[inds]]]

Define two ellipsoids (axis0, axis1, axis2, center)

ep1={{.8,0,0},{0,1,.2},{0,0,1},{0,.5,-2}};
ep0={{.2,0,0},{0,.2,0},{0,0,.2},{0,0,0}};

This is how you do the OpenGL rendering, construct the transformation matrix from the axes and center, and use it as the model-view matrix.

ellipsoid[{a0_,a1_,a2_,center_}]:=
Module[{m,av0=Append[a0,0],av1=Append[a1,0],av2=Append[a2,0],c=Append[center,1]},
m=Transpose[{av0,av1,av2,c}];
Map[Most[m.Append[#,1]]&,verts]];

This is an OpenGL rendering of our scene.  Our final goal is to produce a ray traced version of this from some view point.

Graphics3D[{
EdgeForm[],
Red,
GraphicsComplex[ellipsoid[ep0],Polygon[inds]],
Blue,
GraphicsComplex[ellipsoid[ep1],Polygon[inds]]}]

Here we start the ray tracing

Intersect ray against unit sphere.  Element[ret,Reals] checks the solution of the ray-ellipsoid intersection.

isectSphere[{pt_,dir_}]:=Module[{res,ret}, 
res=First[Solve[(pt+t*dir).(pt+t*dir)==1,t]]; (*overkill, could write it as simple linear equations*)
ret=(t/.First[res]);
If[Element[ret,Reals],ret,Infinity]];

intersection against a sphere is first transform the ray into the sphere-space.  The transformation is the inverse of the transformation as defined by the three axes and a point.

isectEllipsoid[{a0_,a1_,a2_,c_}][{pt_,dir_}]:=Module[{m,im,np,ndir},
m=Transpose[{a0,a1,a2,c}];
im=Inverse[m];
np=Most[im.Append[pt,1]];
ndir=Most[im.Append[dir,0]];
isectSphere[{np,ndir}]];

For every pixel on screen, ray trace the ellpsoids.

px=Module[{eye={0,0,4},pixs,dir,episect0,episect1,t0,t1},
(* turn vectors/points into homogenious coordinates before passing it to isectEllipsoid*)
episect0=isectEllipsoid[Transpose[Append[Transpose[ep0],{0,0,0,1}]]]; 
episect1=isectEllipsoid[Transpose[Append[Transpose[ep1],{0,0,0,1}]]];
pixs=Table[{i,j,1},{i,-1,1,.01},{j,-1,1,.01}];
Reverse[Table[
dir = Normalize[pixs[[i,j]]-eye];
t0=episect0[{eye,dir}];
t1=episect1[{eye,dir}];
If[t0<t1,
If[t0<Infinity,{1,0,0},{1,1,1}],
If[t1<Infinity,{0,0,1},{1,1,1}]],
{j,Length[pixs]},{i,Length[pixs]}]]];

Result

Graphics[Raster[px]]