## 6.4 Disks

The disk is an interesting quadric since it has a particularly straightforward intersection routine that avoids solving the quadratic equation. In pbrt, a Disk is a circular disk of radius at height along the axis.

To describe partial disks, the user may specify a maximum value beyond which the disk is cut off (Figure 6.9). The disk can also be generalized to an annulus by specifying an inner radius, . In parametric form, it is described by

Figure 6.10 is a rendered image of two disks. <<Disk Definition>>=
class Disk { public: <<Disk Public Methods>>
<<Reject disk intersections for rays parallel to the disk’s plane>>
if (Float(di.z) == 0) return {};
Float tShapeHit = (height - Float(oi.z)) / Float(di.z); if (tShapeHit <= 0 || tShapeHit >= tMax) return {};
<<See if hit point is inside disk radii and >>
Point3f pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); Float dist2 = Sqr(pHit.x) + Sqr(pHit.y); if (dist2 > Sqr(radius) || dist2 < Sqr(innerRadius)) return {}; <<Test disk value against >>
Float phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi; if (phi > phiMax) return {};
} SurfaceInteraction InteractionFromIntersection( const QuadricIntersection &isect, Vector3f wo, Float time) const { Point3f pHit = isect.pObj; Float phi = isect.phi; <<Find parametric representation of disk hit>>
Float u = phi / phiMax; Float rHit = std::sqrt(Sqr(pHit.x) + Sqr(pHit.y)); Float v = (radius - rHit) / (radius - innerRadius); Vector3f dpdu(-phiMax * pHit.y, phiMax * pHit.x, 0); Vector3f dpdv = Vector3f(pHit.x, pHit.y, 0) * (innerRadius - radius) / rHit; Normal3f dndu(0, 0, 0), dndv(0, 0, 0);
<<Refine disk intersection point>>
pHit.z = height;
<<Compute error bounds for disk intersection>>
Vector3f pError(0, 0, 0);
bool flipNormal = reverseOrientation ^ transformSwapsHandedness; Vector3f woObject = (*objectFromRender)(wo); return (*renderFromObject)( SurfaceInteraction(Point3fi(pHit, pError), Point2f(u, v), woObject, dpdu, dpdv, dndu, dndv, time, flipNormal));
} bool IntersectP(const Ray &r, Float tMax = Infinity) const { return BasicIntersect(r, tMax).has_value(); } pstd::optional<ShapeSample> Sample(Point2f u) const { Point2f pd = SampleUniformDiskConcentric(u); Point3f pObj(pd.x * radius, pd.y * radius, height); Point3fi pi = (*renderFromObject)(Point3fi(pObj)); Normal3f n = Normalize((*renderFromObject)(Normal3f(0, 0, 1))); if (reverseOrientation) n *= -1; <<Compute for sampled point on disk>>
Float phi = std::atan2(pd.y, pd.x); if (phi < 0) phi += 2 * Pi; Float radiusSample = std::sqrt(Sqr(pObj.x) + Sqr(pObj.y)); Point2f uv(phi / phiMax, (radius - radiusSample) / (radius - innerRadius));
return ShapeSample{Interaction(pi, n, uv), 1 / Area()}; } Float PDF(const Interaction &) const { return 1 / Area(); } pstd::optional<ShapeSample> Sample(const ShapeSampleContext &ctx, Point2f u) const { <<Sample shape by area and compute incident direction wi>>
pstd::optional<ShapeSample> ss = Sample(u); ss->intr.time = ctx.time; Vector3f wi = ss->intr.p() - ctx.p(); if (LengthSquared(wi) == 0) return {}; wi = Normalize(wi);
<<Convert area sampling PDF in ss to solid angle measure>>
ss->pdf /= AbsDot(ss->intr.n, -wi) / DistanceSquared(ctx.p(), ss->intr.p()); if (IsInf(ss->pdf)) return {};
return ss; } PBRT_CPU_GPU Float PDF(const ShapeSampleContext &ctx, Vector3f wi) const { <<Intersect sample ray with shape geometry>>
Ray ray = ctx.SpawnRay(wi); pstd::optional<ShapeIntersection> isect = Intersect(ray); if (!isect) return 0;
<<Compute PDF in solid angle measure from shape intersection point>>
Float pdf = (1 / Area()) / (AbsDot(isect->intr.n, -wi) / DistanceSquared(ctx.p(), isect->intr.p())); if (IsInf(pdf)) pdf = 0;
return pdf; }
private: <<Disk Private Members>>
const Transform *renderFromObject, *objectFromRender; bool reverseOrientation, transformSwapsHandedness; Float height, radius, innerRadius, phiMax;
};

The Disk constructor directly initializes its various member variables from the values passed to it. We have omitted it here because it is trivial.

<<Disk Private Members>>=
const Transform *renderFromObject, *objectFromRender; bool reverseOrientation, transformSwapsHandedness; Float height, radius, innerRadius, phiMax;

### 6.4.1 Area and Bounding

Disks have easily computed surface area, since they are just portions of an annulus:

<<Disk Public Methods>>=
Float Area() const { return phiMax * 0.5f * (Sqr(radius) - Sqr(innerRadius)); }

The bounding method is also quite straightforward; it computes a bounding box centered at the height of the disk along , with extent of radius in both the and directions.

<<Disk Method Definitions>>=

A disk has a single surface normal.

<<Disk Method Definitions>>+=
DirectionCone Disk::NormalBounds() const { Normal3f n = (*renderFromObject)(Normal3f(0, 0, 1)); if (reverseOrientation) n = -n; return DirectionCone(Vector3f(n)); }

### 6.4.2 Intersection Tests

The Disk intersection test methods follow the same form as the earlier quadrics. We omit Intersect(), as it is exactly the same as Sphere::Intersect() and Cylinder::Intersect(), with calls to BasicIntersect() and then InteractionFromIntersection().

The basic intersection test for a ray with a disk is easy. The intersection of the ray with the plane that the disk lies in is found and then the intersection point is checked to see if it lies inside the disk.

<<Disk Public Methods>>+=
pstd::optional<QuadricIntersection> BasicIntersect(const Ray &r, Float tMax) const { <<Transform Ray origin and direction to object space>>  <<Compute plane intersection for disk>>
<<Reject disk intersections for rays parallel to the disk’s plane>>
if (Float(di.z) == 0) return {};
Float tShapeHit = (height - Float(oi.z)) / Float(di.z); if (tShapeHit <= 0 || tShapeHit >= tMax) return {};
<<See if hit point is inside disk radii and >>
Point3f pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); Float dist2 = Sqr(pHit.x) + Sqr(pHit.y); if (dist2 > Sqr(radius) || dist2 < Sqr(innerRadius)) return {}; <<Test disk value against >>
Float phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi; if (phi > phiMax) return {};
}

The first step is to compute the parametric value where the ray intersects the plane that the disk lies in. We want to find such that the component of the ray’s position is equal to the height of the disk. Thus,

and so

The intersection method computes a value and checks to see if it is inside the range of values . If not, the routine can report that there is no intersection.

<<Compute plane intersection for disk>>=
<<Reject disk intersections for rays parallel to the disk’s plane>>
if (Float(di.z) == 0) return {};
Float tShapeHit = (height - Float(oi.z)) / Float(di.z); if (tShapeHit <= 0 || tShapeHit >= tMax) return {};

If the ray is parallel to the disk’s plane (i.e., the component of its direction is zero), no intersection is reported. The case where a ray is both parallel to the disk’s plane and lies within the plane is somewhat ambiguous, but it is most reasonable to define intersecting a disk edge-on as “no intersection.” This case must be handled explicitly so that not-a-number floating-point values are not generated by the following code.

<<Reject disk intersections for rays parallel to the disk’s plane>>=
if (Float(di.z) == 0) return {};

Now the intersection method can compute the point pHit where the ray intersects the plane. Once the plane intersection is known, an invalid intersection is returned if the distance from the hit to the center of the disk is more than Disk::radius or less than Disk::innerRadius. This check can be optimized by computing the squared distance to the center, taking advantage of the fact that the and coordinates of the center point are zero, and the coordinate of pHit is equal to height.

<<See if hit point is inside disk radii and >>=
Point3f pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); Float dist2 = Sqr(pHit.x) + Sqr(pHit.y); if (dist2 > Sqr(radius) || dist2 < Sqr(innerRadius)) return {}; <<Test disk value against >>
Float phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi; if (phi > phiMax) return {};

If the distance check passes, a final test makes sure that the value of the hit point is between zero and , specified by the caller. Inverting the disk’s parameterization gives the same expression for as the other quadric shapes. Because a ray can only intersect a disk once, there is no need to consider a second intersection if this test fails, as was the case with the two earlier quadrics.

<<Test disk value against >>=
Float phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi; if (phi > phiMax) return {};

Finding the SurfaceInteraction corresponding to a disk intersection follows the same process of inverting the parametric representation we have seen before.

<<Disk Public Methods>>+=
SurfaceInteraction InteractionFromIntersection( const QuadricIntersection &isect, Vector3f wo, Float time) const { Point3f pHit = isect.pObj; Float phi = isect.phi; <<Find parametric representation of disk hit>>
Float u = phi / phiMax; Float rHit = std::sqrt(Sqr(pHit.x) + Sqr(pHit.y)); Float v = (radius - rHit) / (radius - innerRadius); Vector3f dpdu(-phiMax * pHit.y, phiMax * pHit.x, 0); Vector3f dpdv = Vector3f(pHit.x, pHit.y, 0) * (innerRadius - radius) / rHit; Normal3f dndu(0, 0, 0), dndv(0, 0, 0);
<<Refine disk intersection point>>
pHit.z = height;
<<Compute error bounds for disk intersection>>
Vector3f pError(0, 0, 0);
bool flipNormal = reverseOrientation ^ transformSwapsHandedness; Vector3f woObject = (*objectFromRender)(wo); return (*renderFromObject)( SurfaceInteraction(Point3fi(pHit, pError), Point2f(u, v), woObject, dpdu, dpdv, dndu, dndv, time, flipNormal));
}

The parameter u is first scaled to reflect the partial disk specified by , and v is computed by inverting the parametric equation. The equations for the partial derivatives at the hit point can be derived with a process similar to that used for the previous quadrics. Because the normal of a disk is the same everywhere, the partial derivatives and are both trivially .

<<Find parametric representation of disk hit>>=
Float u = phi / phiMax; Float rHit = std::sqrt(Sqr(pHit.x) + Sqr(pHit.y)); Float v = (radius - rHit) / (radius - innerRadius); Vector3f dpdu(-phiMax * pHit.y, phiMax * pHit.x, 0); Vector3f dpdv = Vector3f(pHit.x, pHit.y, 0) * (innerRadius - radius) / rHit; Normal3f dndu(0, 0, 0), dndv(0, 0, 0);

As usual, the implementation of IntersectP() is straightforward.

<<Disk Public Methods>>+=
bool IntersectP(const Ray &r, Float tMax = Infinity) const { return BasicIntersect(r, tMax).has_value(); }

### 6.4.3 Sampling

The Disk area sampling method uses a utility routine, SampleUniformDiskConcentric(), that uniformly samples a unit disk. (It is defined in Section A.5.1.) The point that it returns is then scaled by the radius and offset in so that it lies on the disk of a given radius and height. Note that our implementation here does not account for partial disks due to Disk::innerRadius being nonzero or Disk::phiMax being less than . Fixing this bug is left for an exercise at the end of the chapter.

<<Disk Public Methods>>+=
pstd::optional<ShapeSample> Sample(Point2f u) const { Point2f pd = SampleUniformDiskConcentric(u); Point3f pObj(pd.x * radius, pd.y * radius, height); Point3fi pi = (*renderFromObject)(Point3fi(pObj)); Normal3f n = Normalize((*renderFromObject)(Normal3f(0, 0, 1))); if (reverseOrientation) n *= -1; <<Compute for sampled point on disk>>
Float phi = std::atan2(pd.y, pd.x); if (phi < 0) phi += 2 * Pi; Float radiusSample = std::sqrt(Sqr(pObj.x) + Sqr(pObj.y)); Point2f uv(phi / phiMax, (radius - radiusSample) / (radius - innerRadius));
return ShapeSample{Interaction(pi, n, uv), 1 / Area()}; }

The same computation as in the Intersect() method gives the parametric for the sampled point.

<<Compute for sampled point on disk>>=
Float phi = std::atan2(pd.y, pd.x); if (phi < 0) phi += 2 * Pi; Float radiusSample = std::sqrt(Sqr(pObj.x) + Sqr(pObj.y)); Point2f uv(phi / phiMax, (radius - radiusSample) / (radius - innerRadius));

<<Disk Public Methods>>+=
Float PDF(const Interaction &) const { return 1 / Area(); }

We do not provide a specialized solid angle sampling method for disks, but follow the same approach that we did for cylinders, sampling uniformly by area and then computing the probability density to be with respect to solid angle. The implementations of those methods are not included here, as they are the same as they were for cylinders.