## 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:

`Ray`to object space>>

`EFloat`ray coordinate values>>

`t`values>>

`t0`and

`t1`for nearest intersection>>

`SurfaceInteraction`from parametric information>>

`tHit`for quadric intersection>>

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

`Ray`to object space>>=

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 `Float`s, 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`.

`EFloat`ray coordinate values>>

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.

`EFloat`ray coordinate values>>=

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

`t`values>>=

`t0`and

`t1`for nearest intersection>>

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.

`t0`and

`t1`for nearest 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.

`SurfaceInteraction`from parametric information>>=

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.

`tHit`for quadric 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.

`Ray`to object space>>

`EFloat`ray coordinate values>>

`t`values>>

`t0`and

`t1`for nearest intersection>>

### 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