POV-Ray: Documentation: 1.3.10.7 Ray-sphere intersection


POV-Ray : Documentation : 1.3.10.7 Ray-sphere intersection

1.3.10.6 Initializing the raytracer

1.3.10.8 The Trace macro

1.3.10.7 Ray-sphere intersection

#macro calcRaySphereIntersection(P, D, sphereInd)
  #local V = P-Coord[sphereInd][0];
  #local R = Coord[sphereInd][1].x;

  #local DV = vdot(D, V);
  #local D2 = vdot(D, D);
  #local SQ = DV*DV-D2*(vdot(V, V)-R*R);
  #if(SQ < 0) #local Result = -1;
  #else
    #local SQ = sqrt(SQ);
    #local T1 = (-DV+SQ)/D2;
    #local T2 = (-DV-SQ)/D2;
    #local Result = (T1<T2 ? T1 : T2);
  #end
  Result
#end

This is the core of the whole raytracing process.

First let's see how a macro works (if you know it, just skip the following section):

1.3.10.7.1 Inner workings of a #macro

A macro works like a substitution command (similar to the #define macros in the C programming language). The body of the macro is in practice inserted in the place where the macro is called. For example you can use a macro like this:

#macro UnitSphere()
  sphere { 0,1 }
#end

object { UnitSphere() pigment { rgb 1 } }

The result of this code is, in effect, as if you had written:

object { sphere { 0,1 } pigment { rgb 1 } }

Of course there is no reason in making this, as you could have just #declared the UnitSphere as a sphere of radius 1. However, the power of macros kick in when you start using macro parameters. For example:

#macro Sphere(Radius)
  sphere { 0, Radius }
#end

object { Sphere(3) pigment { rgb 1 } }

Now you can use the macro Sphere to create a sphere with the specified radius. Of course this does not make much sense either, as you could just write the sphere primitive directly because it is so short, but this example is intentionally short to show how it works; the macros become very handy when they create something much more complicated than just a sphere.

There is one important difference between macros in POV-Ray and real substitution macros: Any #local statement inside the macro definition will be parsed at the visibility level of the macro only, that is, it will have no effect on the environment where the macro was called from. The following example shows what I am talking about:

#macro Sphere(Radius)
  #local Color = <1,1,1>;
  sphere { 0, Radius pigment { rgb Color } }
#end

#declare Color = <1,0,0>;
object { Sphere(3) }
   // 'Color' is still <1,0,0> here, 
   // thus the following box will be red:
box { -1,1 pigment { rgb Color } }

In the example above, although the macro creates a local identifier called Color and there is an identifier with the same name at global level, the local definition does not affect the global one. Also even if there was not any global definition of Color, the one inside the macro is not seen outside it.

There is one important exception to this, and this is one of the most powerful features of macros (thanks to this they can be used as if they were functions): If an identifier (be it local or global) appears alone in the body of a macro (usually at the end), its value will be passed outside the macro (as if it was a return value). The following example shows how this works:

#macro Factorial(N)
  #local Result = 1;
  #local Ind = 2;
  #while(Ind <= N)
    #local Result = Result*Ind;
    #local Ind = Ind+1;
  #end
  Result
#end

#declare Value = Factorial(5);

Although the identifier Result is local to the macro, its value is passed as if it was a return value because of the last line of the macro (where Result appears alone) and thus the identifier Value will be set to the factorial of 5.

1.3.10.7.2 The ray-sphere intersection macro

Here is again the macro at the beginning of the page so that you do not have to scroll so much in order to see it:

#macro calcRaySphereIntersection(P, D, sphereInd)
  #local V = P-Coord[sphereInd][0];
  #local R = Coord[sphereInd][1].x;

  #local DV = vdot(D, V);
  #local D2 = vdot(D, D);
  #local SQ = DV*DV-D2*(vdot(V, V)-R*R);
  #if(SQ < 0) #local Result = -1;
  #else
    #local SQ = sqrt(SQ);
    #local T1 = (-DV+SQ)/D2;
    #local T2 = (-DV-SQ)/D2;
    #local Result = (T1<T2 ? T1 : T2);
  #end
  Result
#end

The idea behind this macro is that it takes a starting point (ie. the starting point of the ray) a direction vector (the direction where the ray is shot) and an index to the sphere definition array defined previously. The macro returns a factor value; this value expresses how much we have to multiply the direction vector in order to hit the sphere.

This means that if the ray hits the specified sphere, the intersection point will be located at:
StartingPoint + Result*Direction

The return value can be negative, which means that the intersection point was actually behind the starting point. A negative value will be just ignored, as if the ray did not hit anything. We can use this to make a little trick (which may seem obvious when said, but not so obvious when you have to figure it out for yourself): If the ray actually does not hit the sphere, we return just a negative value (does not really matter which).

And how does the macro do it? What is the theory behind those complicated-looking mathematical expressions?

I will use a syntax similar to POV-Ray syntax to express mathematical formulas here since that is probably the easiest way of doing it.

Let's use the following letters:

P = Starting point of the ray
D = Direction of the ray
C = Center of the sphere
R = Radius of the sphere

The theory behind the macro is that we have to see what is the value T for which holds that:

vlength(P+T*D-C) = R

This means: The length of the vector between the center of the sphere (C) and the intersection point (P+T*D) is equal to the radius (R).

If we use an additional letter so that:

V = P-C

then the formula is reduced to:

vlength(T*D+V) = R

which makes our life easier. This formula can be opened as:

(T*D_x+V_x)² + (T*D_y+V_y)² + (T*D_z+V_z)² - R² = 0

Solving T from that is rather trivial math. We get a 2nd order polynomial which has two solutions (I will use the "·" symbol to represent the dot-product of two vectors):

T = (-D·V ± sqrt((D·V)² - D²(V²-R²))) / D²

Note: D² means actually D·D)

When the discriminant (ie. the expression inside the square root) is negative, the ray does not hit the sphere and thus we can return a negative value (the macro returns -1). We must check this in order to avoid the square root of a negative number error; as it has a very logical meaning in this case, the checking is natural.

If the value is positive, there are two solutions (or just one if the value is zero, but that does not really matter here), which corresponds to the two intersection points of the ray with the sphere.

As we get two values, we have to return the one which is smaller (the closest intersection). This is what this portion of the code does:

    #local Result = (T1<T2 ? T1 : T2);

Note: this is an incomplete algorithm: If one value is negative and the other positive (this happens when the starting point is inside the sphere), we would have to return the positive one. The way it is now results in that we will not see the inner surface of the sphere if we put the camera inside one.

For our simple scene this is enough as we do not put our camera inside a sphere nor we have transparent spheres. We could add a check there which looks if one of the values is positive and the other negative and returns the positive one. However, this has an odd and very annoying result (you can try it if you want). This is most probably caused by the inaccuracy of floating point numbers and happens when calculating reflections (the starting point is exactly on the surface of the sphere). We could correct these problems by using epsilon values to get rid of accuracy problems, but in our simple scene this will not be necessary.

1.3.10.6 Initializing the raytracer

1.3.10.8 The Trace macro