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 r at height h along the z axis.

To describe partial disks, the user may specify a maximum phi 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, r Subscript normal i . In parametric form, it is described by

StartLayout 1st Row 1st Column phi 2nd Column equals u phi Subscript normal m normal a normal x Baseline 2nd Row 1st Column x 2nd Column equals left-parenthesis left-parenthesis 1 minus v right-parenthesis r plus v r Subscript normal i Baseline right-parenthesis cosine phi 3rd Row 1st Column y 2nd Column equals left-parenthesis left-parenthesis 1 minus v right-parenthesis r plus v r Subscript normal i Baseline right-parenthesis sine phi 4th Row 1st Column z 2nd Column equals h period EndLayout

Figure 6.9: Basic Setting for the Disk Shape. The disk has radius r and is located at height h along the z axis. A partial disk may be swept by specifying a maximum phi value and an inner radius r Subscript normal i .

Figure 6.10 is a rendered image of two disks.

Figure 6.10: Two Disks. A partial disk is on the left, and a complete disk is on the right.

<<Disk Definition>>= 
class Disk { public: <<Disk Public Methods>> 
Disk(const Transform *renderFromObject, const Transform *objectFromRender, bool reverseOrientation, Float height, Float radius, Float innerRadius, Float phiMax) : renderFromObject(renderFromObject), objectFromRender(objectFromRender), reverseOrientation(reverseOrientation), transformSwapsHandedness(renderFromObject->SwapsHandedness()), height(height), radius(radius), innerRadius(innerRadius), phiMax(Radians(Clamp(phiMax, 0, 360))) {} static Disk *Create(const Transform *renderFromObject, const Transform *objectFromRender, bool reverseOrientation, const ParameterDictionary &parameters, const FileLoc *loc, Allocator alloc); std::string ToString() const; Float Area() const { return phiMax * 0.5f * (Sqr(radius) - Sqr(innerRadius)); } PBRT_CPU_GPU Bounds3f Bounds() const; PBRT_CPU_GPU DirectionCone NormalBounds() const; pstd::optional<ShapeIntersection> Intersect(const Ray &ray, Float tMax = Infinity) const { pstd::optional<QuadricIntersection> isect = BasicIntersect(ray, tMax); if (!isect) return {}; SurfaceInteraction intr = InteractionFromIntersection(*isect, -ray.d, ray.time); return ShapeIntersection{intr, isect->tHit}; } 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 phi Subscript normal m normal a normal x >> 
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 phi value against phi Subscript normal m normal a normal x >> 
Float phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi; if (phi > phiMax) return {};
<<Return QuadricIntersection for disk intersection>> 
return QuadricIntersection{tShapeHit, pHit, phi};
} 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);
<<Return SurfaceInteraction for quadric intersection>> 
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 left-parenthesis u comma v right-parenthesis 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:

upper A equals StartFraction phi Subscript normal m normal a normal x Baseline Over 2 EndFraction left-parenthesis r squared minus r Subscript normal i Superscript 2 Baseline right-parenthesis period

<<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 z , with extent of radius in both the x and y directions.

<<Disk Method Definitions>>= 
Bounds3f Disk::Bounds() const { return (*renderFromObject)( Bounds3f(Point3f(-radius, -radius, height), Point3f( radius, radius, height))); }

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 z equals h 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 phi Subscript normal m normal a normal x >> 
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 phi value against phi Subscript normal m normal a normal x >> 
Float phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi; if (phi > phiMax) return {};
<<Return QuadricIntersection for disk intersection>> 
return QuadricIntersection{tShapeHit, pHit, phi};
}

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

h equals normal o Subscript z Baseline plus t bold d Subscript z

and so

t equals StartFraction h minus normal o Subscript z Baseline Over bold d Subscript z Baseline EndFraction period

The intersection method computes a t value and checks to see if it is inside the range of values left-parenthesis 0 comma monospace t monospace upper M monospace a monospace x right-parenthesis . 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 z 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 x and y coordinates of the center point left-parenthesis 0 comma 0 comma monospace h monospace e monospace i monospace g monospace h monospace t right-parenthesis are zero, and the z coordinate of pHit is equal to height.

<<See if hit point is inside disk radii and phi Subscript normal m normal a normal x >>= 
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 phi value against phi Subscript normal m normal a normal x >> 
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 phi value of the hit point is between zero and phi Subscript normal m normal a normal x , specified by the caller. Inverting the disk’s parameterization gives the same expression for phi 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 phi value against phi Subscript normal m normal a normal x >>= 
Float phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi; if (phi > phiMax) return {};

<<Return QuadricIntersection for disk intersection>>= 
return QuadricIntersection{tShapeHit, pHit, phi};

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);
<<Return SurfaceInteraction for quadric intersection>> 
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 phi Subscript normal m normal a normal x , 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 partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v are both trivially left-parenthesis 0 comma 0 comma 0 right-parenthesis .

<<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 z 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 2 pi . 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 left-parenthesis u comma v right-parenthesis 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 left-parenthesis u comma v right-parenthesis for the sampled point.

<<Compute left-parenthesis u comma v right-parenthesis 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.