## 6.3 Cylinders

Another useful quadric is the cylinder; pbrt provides a cylinder Shape that is centered around the axis. The user can supply a minimum and maximum value for the cylinder, as well as a radius and maximum sweep value (Figure 6.7).

<<Cylinder Definition>>=
class Cylinder { public: <<Cylinder Public Methods>>
Cylinder(const Transform *renderFromObj, const Transform *objFromRender, bool reverseOrientation, Float radius, Float zMin, Float zMax, Float phiMax); static Cylinder *Create(const Transform *renderFromObject, const Transform *objectFromRender, bool reverseOrientation, const ParameterDictionary &parameters, const FileLoc *loc, Allocator alloc); PBRT_CPU_GPU Bounds3f Bounds() const; std::string ToString() const; Float Area() const { return (zMax - zMin) * radius * phiMax; } DirectionCone NormalBounds() const { return DirectionCone::EntireSphere(); } 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 { Float phi; Point3f pHit; <<Transform Ray origin and direction to object space>>  <<Solve quadratic equation to find cylinder t0 and t1 values>>
Interval t0, t1; <<Compute cylinder quadratic coefficients>>
Interval a = Sqr(di.x) + Sqr(di.y); Interval b = 2 * (di.x * oi.x + di.y * oi.y); Interval c = Sqr(oi.x) + Sqr(oi.y) - Sqr(Interval(radius));
Interval f = b / (2 * a); Interval vx = oi.x - f * di.x, vy = oi.y - f * di.y; Interval length = Sqrt(Sqr(vx) + Sqr(vy)); Interval discrim = 4 * a * (Interval(radius) + length) * (Interval(radius) - length); if (discrim.LowerBound() < 0) return {};
Interval rootDiscrim = Sqrt(discrim); Interval q; if ((Float)b < 0) q = -.5f * (b - rootDiscrim); else q = -.5f * (b + rootDiscrim); t0 = q / a; t1 = c / q; <<Swap quadratic values so that t0 is the lesser>>
if (t0.LowerBound() > t1.LowerBound()) pstd::swap(t0, t1);
<<Check quadric shape t0 and t1 for nearest intersection>>
if (t0.UpperBound() > tMax || t1.LowerBound() <= 0) return {}; Interval tShapeHit = t0; if (tShapeHit.LowerBound() <= 0) { tShapeHit = t1; if (tShapeHit.UpperBound() > tMax) return {}; }
<<Compute cylinder hit point and >>
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>>
phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi;
<<Test cylinder intersection against clipping parameters>>
if (pHit.z < zMin || pHit.z > zMax || phi > phiMax) { if (tShapeHit == t1) return {}; tShapeHit = t1; if (t1.UpperBound() > tMax) return {}; <<Compute cylinder hit point and >>
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>>
phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi;
if (pHit.z < zMin || pHit.z > zMax || phi > phiMax) return {}; }
} bool IntersectP(const Ray &r, Float tMax = Infinity) const { return BasicIntersect(r, tMax).has_value(); } SurfaceInteraction InteractionFromIntersection( const QuadricIntersection &isect, Vector3f wo, Float time) const { Point3f pHit = isect.pObj; Float phi = isect.phi; <<Find parametric representation of cylinder hit>>
Float u = phi / phiMax; Float v = (pHit.z - zMin) / (zMax - zMin); <<Compute cylinder and >>
Vector3f dpdu(-phiMax * pHit.y, phiMax * pHit.x, 0); Vector3f dpdv(0, 0, zMax - zMin);
<<Compute cylinder and >>
Vector3f d2Pduu = -phiMax * phiMax * Vector3f(pHit.x, pHit.y, 0); Vector3f d2Pduv(0, 0, 0), d2Pdvv(0, 0, 0); <<Compute coefficients for fundamental forms>>
Float E = Dot(dpdu, dpdu), F = Dot(dpdu, dpdv), G = Dot(dpdv, dpdv); Vector3f n = Normalize(Cross(dpdu, dpdv)); Float e = Dot(n, d2Pduu), f = Dot(n, d2Pduv), g = Dot(n, d2Pdvv);
<<Compute and from fundamental form coefficients>>
Float EGF2 = DifferenceOfProducts(E, G, F, F); Float invEGF2 = (EGF2 == 0) ? Float(0) : 1 / EGF2; Normal3f dndu = Normal3f((f * F - e * G) * invEGF2 * dpdu + (e * F - f * E) * invEGF2 * dpdv); Normal3f dndv = Normal3f((g * F - f * G) * invEGF2 * dpdu + (f * F - g * E) * invEGF2 * dpdv);
<<Compute error bounds for cylinder intersection>>
Vector3f pError = gamma(3) * Abs(Vector3f(pHit.x, pHit.y, 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));
} pstd::optional<ShapeSample> Sample(Point2f u) const { Float z = Lerp(u[0], zMin, zMax); Float phi = u[1] * phiMax; <<Compute cylinder sample position pi and normal n from and >>
Point3f pObj = Point3f(radius * std::cos(phi), radius * std::sin(phi), z); <<Reproject pObj to cylinder surface and compute pObjError>>
Point3fi pi = (*renderFromObject)(Point3fi(pObj, pObjError)); Normal3f n = Normalize((*renderFromObject)(Normal3f(pObj.x, pObj.y, 0))); if (reverseOrientation) n *= -1;
Point2f uv(phi / phiMax, (pObj.z - zMin) / (zMax - zMin)); 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; } 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: <<Cylinder Private Members>>
const Transform *renderFromObject, *objectFromRender; bool reverseOrientation, transformSwapsHandedness; Float radius, zMin, zMax, phiMax;
};

In parametric form, a cylinder is described by the following equations:

Figure 6.8 shows a rendered image of two cylinders. Like the sphere image, the right cylinder is a complete cylinder, while the left one is a partial cylinder because it has a value less than .

Similar to the Sphere constructor, the Cylinder constructor takes transformations that define its object space and the parameters that define the cylinder itself. Its constructor just initializes the corresponding member variables, so we will not include it here.

<<Cylinder Public Methods>>=
Cylinder(const Transform *renderFromObj, const Transform *objFromRender, bool reverseOrientation, Float radius, Float zMin, Float zMax, Float phiMax);

<<Cylinder Private Members>>=
const Transform *renderFromObject, *objectFromRender; bool reverseOrientation, transformSwapsHandedness; Float radius, zMin, zMax, phiMax;

### 6.3.1 Area and Bounding

A cylinder is a rolled-up rectangle. If you unroll the rectangle, its height is , and its width is :

<<Cylinder Public Methods>>+=
Float Area() const { return (zMax - zMin) * radius * phiMax; }

As was done with the sphere, the cylinder’s spatial bounding method computes a conservative bounding box using the range but does not take into account the maximum .

<<Cylinder Method Definitions>>=

Its surface normal bounding function is conservative in two ways: not only does it not account for , but the actual set of normals of a cylinder can be described by a circle on the sphere of all directions. However, DirectionCone’s representation is not able to bound such a distribution more tightly than with the entire sphere of directions, and so that is the bound that is returned.

<<Cylinder Public Methods>>+=
DirectionCone NormalBounds() const { return DirectionCone::EntireSphere(); }

### 6.3.2 Intersection Tests

Also similar to the sphere (and for similar reasons), Cylinder provides a BasicIntersect() method that returns a QuadricIntersection as well as an InteractionFromIntersection() method that converts that to a full SurfaceInteraction. Given these, the Intersect() method is again a simple composition of them. (If pbrt used virtual functions, a design alternative would be to have a QuadricShape class that provided a default Intersect() method and left BasicIntersect() and InteractionFromIntersection() as pure virtual functions for subclasses to implement.)

<<Cylinder Public Methods>>+=
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}; }

The form of the BasicIntersect() method also parallels the sphere’s, computing appropriate quadratic coefficients, solving the quadratic equation, and then handling the various cases for partial cylinders. A number of fragments can be reused from the Sphere’s implementation.

<<Cylinder Public Methods>>+=
pstd::optional<QuadricIntersection> BasicIntersect(const Ray &r, Float tMax) const { Float phi; Point3f pHit; <<Transform Ray origin and direction to object space>>  <<Solve quadratic equation to find cylinder t0 and t1 values>>
Interval t0, t1; <<Compute cylinder quadratic coefficients>>
Interval a = Sqr(di.x) + Sqr(di.y); Interval b = 2 * (di.x * oi.x + di.y * oi.y); Interval c = Sqr(oi.x) + Sqr(oi.y) - Sqr(Interval(radius));
Interval f = b / (2 * a); Interval vx = oi.x - f * di.x, vy = oi.y - f * di.y; Interval length = Sqrt(Sqr(vx) + Sqr(vy)); Interval discrim = 4 * a * (Interval(radius) + length) * (Interval(radius) - length); if (discrim.LowerBound() < 0) return {};
Interval rootDiscrim = Sqrt(discrim); Interval q; if ((Float)b < 0) q = -.5f * (b - rootDiscrim); else q = -.5f * (b + rootDiscrim); t0 = q / a; t1 = c / q; <<Swap quadratic values so that t0 is the lesser>>
if (t0.LowerBound() > t1.LowerBound()) pstd::swap(t0, t1);
<<Check quadric shape t0 and t1 for nearest intersection>>
if (t0.UpperBound() > tMax || t1.LowerBound() <= 0) return {}; Interval tShapeHit = t0; if (tShapeHit.LowerBound() <= 0) { tShapeHit = t1; if (tShapeHit.UpperBound() > tMax) return {}; }
<<Compute cylinder hit point and >>
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>>
phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi;
<<Test cylinder intersection against clipping parameters>>
if (pHit.z < zMin || pHit.z > zMax || phi > phiMax) { if (tShapeHit == t1) return {}; tShapeHit = t1; if (t1.UpperBound() > tMax) return {}; <<Compute cylinder hit point and >>
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>>
phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi;
if (pHit.z < zMin || pHit.z > zMax || phi > phiMax) return {}; }
}

As before, the fragment that computes the quadratic discriminant, <<Compute cylinder quadratic discriminant discrim>>, is defined in Section 6.8.3 after topics related to floating-point accuracy have been discussed.

<<Solve quadratic equation to find cylinder t0 and t1 values>>=
Interval t0, t1; <<Compute cylinder quadratic coefficients>>
Interval a = Sqr(di.x) + Sqr(di.y); Interval b = 2 * (di.x * oi.x + di.y * oi.y); Interval c = Sqr(oi.x) + Sqr(oi.y) - Sqr(Interval(radius));
Interval f = b / (2 * a); Interval vx = oi.x - f * di.x, vy = oi.y - f * di.y; Interval length = Sqrt(Sqr(vx) + Sqr(vy)); Interval discrim = 4 * a * (Interval(radius) + length) * (Interval(radius) - length); if (discrim.LowerBound() < 0) return {};
Interval rootDiscrim = Sqrt(discrim); Interval q; if ((Float)b < 0) q = -.5f * (b - rootDiscrim); else q = -.5f * (b + rootDiscrim); t0 = q / a; t1 = c / q; <<Swap quadratic values so that t0 is the lesser>>
if (t0.LowerBound() > t1.LowerBound()) pstd::swap(t0, t1);

As with spheres, the ray–cylinder intersection formula can be found by substituting the ray equation into the cylinder’s implicit equation. The implicit equation for an infinitely long cylinder centered on the axis with radius  is

Substituting the ray equation, Equation (3.4), we have

When we expand this equation and find the coefficients of the quadratic equation , we have

Interval a = Sqr(di.x) + Sqr(di.y); Interval b = 2 * (di.x * oi.x + di.y * oi.y); Interval c = Sqr(oi.x) + Sqr(oi.y) - Sqr(Interval(radius));

As with spheres, the implementation refines the computed intersection point to reduce the rounding error in the point computed by evaluating the ray equation; see Section 6.8.5. Afterward, we invert the parametric description of the cylinder to compute from and ; it turns out that the result is the same as for the sphere.

<<Compute cylinder hit point and >>=
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>>
phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi;

The next step in the intersection method makes sure that the hit is in the specified range and that the angle is acceptable. If not, it rejects the hit and checks if it has not already been considered—these tests resemble the conditional logic in Sphere::Intersect().

<<Test cylinder intersection against clipping parameters>>=
if (pHit.z < zMin || pHit.z > zMax || phi > phiMax) { if (tShapeHit == t1) return {}; tShapeHit = t1; if (t1.UpperBound() > tMax) return {}; <<Compute cylinder hit point and >>
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>>
phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi;
if (pHit.z < zMin || pHit.z > zMax || phi > phiMax) return {}; }

For a successful intersection, the same three values suffice to provide enough information to later compute the corresponding SurfaceInteraction.

As with the sphere, IntersectP()’s implementation is a simple wrapper around BasicIntersect().

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

InteractionFromIntersection() computes all the quantities needed to initialize a SurfaceInteraction from a cylinder’s QuadricIntersection.

<<Cylinder Public Methods>>+=
SurfaceInteraction InteractionFromIntersection( const QuadricIntersection &isect, Vector3f wo, Float time) const { Point3f pHit = isect.pObj; Float phi = isect.phi; <<Find parametric representation of cylinder hit>>
Float u = phi / phiMax; Float v = (pHit.z - zMin) / (zMax - zMin); <<Compute cylinder and >>
Vector3f dpdu(-phiMax * pHit.y, phiMax * pHit.x, 0); Vector3f dpdv(0, 0, zMax - zMin);
<<Compute cylinder and >>
Vector3f d2Pduu = -phiMax * phiMax * Vector3f(pHit.x, pHit.y, 0); Vector3f d2Pduv(0, 0, 0), d2Pdvv(0, 0, 0); <<Compute coefficients for fundamental forms>>
Float E = Dot(dpdu, dpdu), F = Dot(dpdu, dpdv), G = Dot(dpdv, dpdv); Vector3f n = Normalize(Cross(dpdu, dpdv)); Float e = Dot(n, d2Pduu), f = Dot(n, d2Pduv), g = Dot(n, d2Pdvv);
<<Compute and from fundamental form coefficients>>
Float EGF2 = DifferenceOfProducts(E, G, F, F); Float invEGF2 = (EGF2 == 0) ? Float(0) : 1 / EGF2; Normal3f dndu = Normal3f((f * F - e * G) * invEGF2 * dpdu + (e * F - f * E) * invEGF2 * dpdv); Normal3f dndv = Normal3f((g * F - f * G) * invEGF2 * dpdu + (f * F - g * E) * invEGF2 * dpdv);
<<Compute error bounds for cylinder intersection>>
Vector3f pError = gamma(3) * Abs(Vector3f(pHit.x, pHit.y, 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));
}

Again the parametric value is computed by scaling to lie between 0 and 1. Inversion of the parametric equation for the cylinder’s value gives the parametric coordinate.

<<Find parametric representation of cylinder hit>>=
Float u = phi / phiMax; Float v = (pHit.z - zMin) / (zMax - zMin); <<Compute cylinder and >>
Vector3f dpdu(-phiMax * pHit.y, phiMax * pHit.x, 0); Vector3f dpdv(0, 0, zMax - zMin);
<<Compute cylinder and >>
Vector3f d2Pduu = -phiMax * phiMax * Vector3f(pHit.x, pHit.y, 0); Vector3f d2Pduv(0, 0, 0), d2Pdvv(0, 0, 0); <<Compute coefficients for fundamental forms>>
Float E = Dot(dpdu, dpdu), F = Dot(dpdu, dpdv), G = Dot(dpdv, dpdv); Vector3f n = Normalize(Cross(dpdu, dpdv)); Float e = Dot(n, d2Pduu), f = Dot(n, d2Pduv), g = Dot(n, d2Pdvv);
<<Compute and from fundamental form coefficients>>
Float EGF2 = DifferenceOfProducts(E, G, F, F); Float invEGF2 = (EGF2 == 0) ? Float(0) : 1 / EGF2; Normal3f dndu = Normal3f((f * F - e * G) * invEGF2 * dpdu + (e * F - f * E) * invEGF2 * dpdv); Normal3f dndv = Normal3f((g * F - f * G) * invEGF2 * dpdu + (f * F - g * E) * invEGF2 * dpdv);

The partial derivatives for a cylinder are easy to derive:

<<Compute cylinder and >>=
Vector3f dpdu(-phiMax * pHit.y, phiMax * pHit.x, 0); Vector3f dpdv(0, 0, zMax - zMin);

We again use the Weingarten equations to compute the parametric partial derivatives of the cylinder normal. The relevant partial derivatives are

<<Compute cylinder and >>=
Vector3f d2Pduu = -phiMax * phiMax * Vector3f(pHit.x, pHit.y, 0); Vector3f d2Pduv(0, 0, 0), d2Pdvv(0, 0, 0); <<Compute coefficients for fundamental forms>>
Float E = Dot(dpdu, dpdu), F = Dot(dpdu, dpdv), G = Dot(dpdv, dpdv); Vector3f n = Normalize(Cross(dpdu, dpdv)); Float e = Dot(n, d2Pduu), f = Dot(n, d2Pduv), g = Dot(n, d2Pdvv);
<<Compute and from fundamental form coefficients>>
Float EGF2 = DifferenceOfProducts(E, G, F, F); Float invEGF2 = (EGF2 == 0) ? Float(0) : 1 / EGF2; Normal3f dndu = Normal3f((f * F - e * G) * invEGF2 * dpdu + (e * F - f * E) * invEGF2 * dpdv); Normal3f dndv = Normal3f((g * F - f * G) * invEGF2 * dpdu + (f * F - g * E) * invEGF2 * dpdv);

### 6.3.3 Sampling

Uniformly sampling the surface area of a cylinder is straightforward: uniform sampling of the height and give uniform area sampling. Intuitively, it can be understood that this approach works because a cylinder is just a rolled-up rectangle.

<<Cylinder Public Methods>>+=
pstd::optional<ShapeSample> Sample(Point2f u) const { Float z = Lerp(u[0], zMin, zMax); Float phi = u[1] * phiMax; <<Compute cylinder sample position pi and normal n from and >>
Point3f pObj = Point3f(radius * std::cos(phi), radius * std::sin(phi), z); <<Reproject pObj to cylinder surface and compute pObjError>>
Point3fi pi = (*renderFromObject)(Point3fi(pObj, pObjError)); Normal3f n = Normalize((*renderFromObject)(Normal3f(pObj.x, pObj.y, 0))); if (reverseOrientation) n *= -1;
Point2f uv(phi / phiMax, (pObj.z - zMin) / (zMax - zMin)); return ShapeSample{Interaction(pi, n, uv), 1 / Area()}; }

Given and , the corresponding object-space position and normal are easily found.

<<Compute cylinder sample position pi and normal n from and >>=
Point3f pObj = Point3f(radius * std::cos(phi), radius * std::sin(phi), z); <<Reproject pObj to cylinder surface and compute pObjError>>
Point3fi pi = (*renderFromObject)(Point3fi(pObj, pObjError)); Normal3f n = Normalize((*renderFromObject)(Normal3f(pObj.x, pObj.y, 0))); if (reverseOrientation) n *= -1;

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

Unlike the Sphere, pbrt’s Cylinder does not have a specialized solid angle sampling method. Instead, it samples a point on the cylinder uniformly by area without making use of the reference point before converting the area density for that point to a solid angle density before returning it. Both the Sample() and PDF() methods can be implemented using the same fragments that were used for solid angle sampling of reference points inside spheres.

<<Cylinder Public Methods>>+=
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; }

<<Cylinder Public Methods>>+=
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; }