3.2 Spheres

Spheres are a special case of a general type of surface called quadrics—surfaces described by quadratic polynomials in , , and . They are the simplest type of curved surface that is useful to a ray tracer and are a good starting point for general ray intersection routines. pbrt supports six types of quadrics: spheres, cones, disks (a special case of a cone), cylinders, hyperboloids, and paraboloids.

Many surfaces can be described in one of two main ways: in implicit form and in parametric form. An implicit function describes a 3D surface as

The set of all points (, , ) that fulfill this condition defines the surface. For a unit sphere at the origin, the familiar implicit equation is . Only the set of points one unit from the origin satisfies this constraint, giving the unit sphere’s surface.

Many surfaces can also be described parametrically using a function to map 2D points to 3D points on the surface. For example, a sphere of radius can be described as a function of 2D spherical coordinates , where ranges from to and ranges from to (Figure 3.4):

We can transform this function into a function over and also generalize it slightly to allow partial spheres that only sweep out and with the substitution

This form is particularly useful for texture mapping, where it can be directly used to map a texture defined over to the sphere. Figure 3.5 shows an image of two spheres; a grid image map has been used to show the parameterization.

As we describe the implementation of the sphere shape, we will make use of both the implicit and parametric descriptions of the shape, depending on which is a more natural way to approach the particular problem we’re facing.

The Sphere class represents a sphere that is centered at the origin in object space. Its implementation is in the files shapes/sphere.h and shapes/sphere.cpp.

To place a sphere elsewhere in the scene, the user must apply an appropriate transformation when specifying the sphere in the input file. It takes both the object-to-world and world-to-object transformations as parameters to the constructor, passing them along to the parent Shape constructor.

The radius of the sphere can have an arbitrary positive value, and the sphere’s extent can be truncated in two different ways. First, minimum and maximum values may be set; the parts of the sphere below and above these planes, respectively, are cut off. Second, considering the parameterization of the sphere in spherical coordinates, a maximum value can be set. The sphere sweeps out values from 0 to the given such that the section of the sphere with spherical values above is also removed.

3.2.1 Bounding

Computing an object space bounding box for a sphere is straightforward. The implementation here uses the values of and provided by the user to tighten up the bound when less than an entire sphere is being rendered. However, it doesn’t do the extra work to compute a tighter bounding box when is less than . This improvement is left as an exercise.

3.2.2 Intersection Tests

The task of deriving a ray–sphere intersection test is simplified by the fact that the sphere is centered at the origin. However, if the sphere has been transformed to another position in world space, then it is necessary to transform rays to object space before intersecting them with the sphere, using the world-to-object transformation. Given a ray in object space, the intersection computation can be performed in object space instead. ^†

The following fragment shows the entire intersection method:

First, the given world space ray is transformed to the sphere’s object space. The remainder of the intersection test will take place in that coordinate system. The oErr and dErr variables respectively bound the floating-point round-off error in the transformed ray’s origin and direction that was introduced by applying the transformation. (See Section 3.9 for more information about floating-point arithmetic and its implications for accurate ray intersection calculations.)

If a sphere is centered at the origin with radius , its implicit representation is

By substituting the parametric representation of the ray from Equation (2.3) into the implicit sphere equation, we have

Note that all elements of this equation besides are known values. The values where the equation holds give the parametric positions along the ray where the implicit sphere equation holds and thus the points along the ray where it intersects the sphere. We can expand this equation and gather the coefficients for a general quadratic equation in ,

where ^†

This result directly translates to this fragment of source code. Note that in this code, instances of the EFloat class, not Floats, are used to represent floating-point values. EFloat tracks accumulated floating-point rounding error; its use is discussed in Section 3.9. For now, it can just be read as being equivalent to Float.

The ray origin and direction values used in the intersection test are initialized with the floating-point error bounds from transforming the ray to object space.

There are two possible solutions to the quadratic equation, giving zero, one, or two nonimaginary values where the ray intersects the sphere.

The Quadratic() utility function solves a quadratic equation, returning false if there are no real solutions and returning true and setting t0 and t1 appropriately if there are solutions. It is defined later in Section 3.9.4, where we discuss how to implement it robustly using floating-point arithmetic.

The computed parametric distances t0 and t1 track uncertainty due to errors in the original ray parameters and errors accrued in Quadratic(); the lower and upper range of the uncertainty interval can be queried using the methods EFloat::LowerBound() and EFloat::UpperBound().

The fragment <<Check quadric shape t0 and t1 for nearest intersection>> takes the two intersection values and determines which, if any, is the closest valid intersection. For an intersection to be valid, its value must be greater than zero and less than ray.tMax. The following code uses the error intervals provided by the EFloat class and only accepts intersections that are unequivocally in the range .

Since is guaranteed to be less than or equal to (and is less than tMax), then if is greater than tMax or is less than , it is certain that both intersections are out of the range of interest. Otherwise, is the tentative hit value. It may be less than , however, in which case we ignore it and try . If that is also out of range, we have no valid intersection. If there is an intersection, then tShapeHit is initialized to hold the parametric value for the intersection.

Given the parametric distance along the ray to the intersection with a full sphere, the intersection point pHit can be computed as that offset along the ray.

It is next necessary to handle partial spheres with clipped or ranges—intersections that are in clipped areas must be ignored. The implementation starts by computing the value for the hit point. Using the parametric representation of the sphere,

so . It is necessary to remap the result of the standard library’s std::atan() function to a value between and , to match the sphere’s original definition.

Due to floating-point precision limitations, this computed intersection point pHit may lie a bit to one side of the actual sphere surface; the <<Refine sphere intersection point>> fragment, which is defined in Section 3.9.4, improves the precision of this value.

The hit point can now be tested against the specified minima and maxima for and . One subtlety is that it’s important to skip the tests if the range includes the entire sphere; the computed pHit.z value may be slightly out of the range due to floating-point round-off, so we should only perform this test when the user expects the sphere to be partially incomplete. If the intersection wasn’t actually valid, the routine tries again with .

At this point in the routine, it is certain that the ray hits the sphere. The method next computes and values by scaling the previously computed value for the hit to lie between 0 and 1 and by computing a value between 0 and 1 for the hit point, based on the range of values for the given sphere. Then it finds the parametric partial derivatives of position and and surface normal and .

Computing the partial derivatives of a point on the sphere is a short exercise in algebra. Here we will show how the component of , , is calculated; the other components are found similarly. Using the parametric definition of the sphere, we have

Using a substitution based on the parametric definition of the sphere’s coordinate, this simplifies to

Similarly,

and

A similar process gives . The complete result is

3.2.3 Partial Derivatives of Normal Vectors

It is also useful to determine how the normal changes as we move along the surface in the and directions. For example, the antialiasing techniques in Chapter 10 are dependent on this information to antialias textures on objects that are seen reflected in curved surfaces. The differential changes in normal and are given by the Weingarten equations from differential geometry:

where , , and are coefficients of the first fundamental form and are given by

These are easily computed with the and values found earlier. The , , and are coefficients of the second fundamental form,

The two fundamental forms capture elementary metric properties of a surface, including notions of distance, angle, and curvature; see a differential geometry textbook such as Gray (1993) for details. To find , , and , it is necessary to compute the second-order partial derivatives and so on.

For spheres, a little more algebra gives the second derivatives:

3.2.4 SurfaceInteraction Initialization

Having computed the surface parameterization and all the relevant partial derivatives, the SurfaceInteraction structure can be initialized with the geometric information for this intersection. The pError value passed to the SurfaceInteraction constructor bounds the rounding error in the computed pHit point. It is initialized in the fragment <<Compute error bounds for sphere intersection>>, which is defined later, in Section 3.9.4.

Since there is an intersection, the tHit parameter to the Intersect() method is updated with the parametric hit distance along the ray, which was stored in tShapeHit. Updating *tHit allows subsequent intersection tests to terminate early if the potential hit would be farther away than the existing intersection.

A natural question to ask at this point is, “What effect does the world-to-object transformation have on the correct parametric distance to return?” Indeed, the intersection method has found a parametric distance to the intersection for the object space ray, which may have been translated, rotated, scaled, or worse when it was transformed from world space. However, it can be shown that the parametric distance to an intersection in object space is exactly the same as it would have been if the ray was left in world space and the intersection had been done there and, thus, tHit can be set directly. Note that if the object space ray’s direction had been normalized after the transformation, then this would no longer be the case and a correction factor related to the unnormalized ray’s length would be needed. This is another motivation for not normalizing the object space ray’s direction vector after transformation.

The Sphere::IntersectP() routine is almost identical to Sphere::Intersect(), but it does not initialize the SurfaceInteraction structure. Because the Intersect() and IntersectP() methods are always so closely related, in the following we will not show implementations of IntersectP() for the remaining shapes.

3.2.5 Surface Area

To compute the surface area of quadrics, we use a standard formula from integral calculus. If a curve from to is revolved around the axis, the surface area of the resulting swept surface is

where denotes the derivative . ^† Since most of our surfaces of revolution are only partially swept around the axis, we will instead use the formula

The sphere is a surface of revolution of a circular arc. The function that defines the profile curve along the axis of the sphere is

and its derivative is

Recall that the sphere is clipped at and . The surface area is therefore

For the full sphere , , and , so we have the standard formula , confirming that the formula makes sense.