6.3 Cylinders

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

Figure 6.7: Basic Setting for the Cylinder Shape. The cylinder has a radius of  r and covers a range along the z axis. A partial cylinder may be swept by specifying a maximum phi value.

<<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));
<<Compute cylinder quadratic discriminant discrim>> 
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 {};
<<Compute quadratic t values>> 
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 t 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 phi >> 
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>> 
Float hitRad = std::sqrt(Sqr(pHit.x) + Sqr(pHit.y)); pHit.x *= radius / hitRad; pHit.y *= radius / hitRad;
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 phi >> 
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>> 
Float hitRad = std::sqrt(Sqr(pHit.x) + Sqr(pHit.y)); pHit.x *= radius / hitRad; pHit.y *= radius / hitRad;
phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi;
if (pHit.z < zMin || pHit.z > zMax || phi > phiMax) return {}; }
<<Return QuadricIntersection for cylinder intersection>> 
return QuadricIntersection{Float(tShapeHit), pHit, phi};
} 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 partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v >> 
Vector3f dpdu(-phiMax * pHit.y, phiMax * pHit.x, 0); Vector3f dpdv(0, 0, zMax - zMin);
<<Compute cylinder partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v >> 
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 partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v 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));
<<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));
} 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 z and phi >> 
Point3f pObj = Point3f(radius * std::cos(phi), radius * std::sin(phi), z); <<Reproject pObj to cylinder surface and compute pObjError>> 
Float hitRad = std::sqrt(Sqr(pObj.x) + Sqr(pObj.y)); pObj.x *= radius / hitRad; pObj.y *= radius / hitRad; Vector3f pObjError = gamma(3) * Abs(Vector3f(pObj.x, pObj.y, 0));
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:

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 r cosine phi 3rd Row 1st Column y 2nd Column equals r sine phi 4th Row 1st Column z 2nd Column equals z Subscript normal m normal i normal n Baseline plus v left-parenthesis z Subscript normal m normal a normal x Baseline minus z Subscript normal m normal i normal n Baseline right-parenthesis period EndLayout

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 phi Subscript normal m normal a normal x value less than  2 pi .

Figure 6.8: Two Cylinders. A partial cylinder is on the left, and a complete cylinder is on the right.

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 z Subscript normal m normal a normal x Baseline minus z Subscript normal m normal i normal n , and its width is r phi Subscript normal m normal a normal x :

<<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 z range but does not take into account the maximum phi .

<<Cylinder Method Definitions>>= 
Bounds3f Cylinder::Bounds() const { return (*renderFromObject)(Bounds3f({-radius, -radius, zMin}, { radius, radius, zMax})); }

Its surface normal bounding function is conservative in two ways: not only does it not account for phi Subscript normal m normal a normal x Baseline less-than 2 pi , 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));
<<Compute cylinder quadratic discriminant discrim>> 
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 {};
<<Compute quadratic t values>> 
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 t 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 phi >> 
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>> 
Float hitRad = std::sqrt(Sqr(pHit.x) + Sqr(pHit.y)); pHit.x *= radius / hitRad; pHit.y *= radius / hitRad;
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 phi >> 
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>> 
Float hitRad = std::sqrt(Sqr(pHit.x) + Sqr(pHit.y)); pHit.x *= radius / hitRad; pHit.y *= radius / hitRad;
phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi;
if (pHit.z < zMin || pHit.z > zMax || phi > phiMax) return {}; }
<<Return QuadricIntersection for cylinder intersection>> 
return QuadricIntersection{Float(tShapeHit), pHit, phi};
}

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));
<<Compute cylinder quadratic discriminant discrim>> 
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 {};
<<Compute quadratic t values>> 
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 t 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 z axis with radius r  is

x squared plus y squared minus r squared equals 0 period

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

left-parenthesis normal o Subscript x Baseline plus t bold d Subscript x Baseline right-parenthesis squared plus left-parenthesis normal o Subscript y Baseline plus t bold d Subscript y Baseline right-parenthesis squared equals r squared period

When we expand this equation and find the coefficients of the quadratic equation a t squared plus b t plus c equals 0 , we have

StartLayout 1st Row 1st Column a 2nd Column equals bold d Subscript x Superscript 2 Baseline plus bold d Subscript y Superscript 2 Baseline 2nd Row 1st Column b 2nd Column equals 2 left-parenthesis bold d Subscript x Baseline normal o Subscript x Baseline plus bold d Subscript y Baseline normal o Subscript y Baseline right-parenthesis 3rd Row 1st Column c 2nd Column equals normal o Subscript x Superscript 2 Baseline plus normal o Subscript y Superscript 2 Baseline minus r squared period EndLayout

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

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 phi from x and y ; it turns out that the result is the same as for the sphere.

<<Compute cylinder hit point and phi >>= 
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>> 
Float hitRad = std::sqrt(Sqr(pHit.x) + Sqr(pHit.y)); pHit.x *= radius / hitRad; pHit.y *= radius / hitRad;
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 z range and that the angle phi is acceptable. If not, it rejects the hit and checks t 1 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 phi >> 
pHit = Point3f(oi) + (Float)tShapeHit * Vector3f(di); <<Refine cylinder intersection point>> 
Float hitRad = std::sqrt(Sqr(pHit.x) + Sqr(pHit.y)); pHit.x *= radius / hitRad; pHit.y *= radius / hitRad;
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.

<<Return QuadricIntersection for cylinder intersection>>= 
return QuadricIntersection{Float(tShapeHit), pHit, phi};

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 partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v >> 
Vector3f dpdu(-phiMax * pHit.y, phiMax * pHit.x, 0); Vector3f dpdv(0, 0, zMax - zMin);
<<Compute cylinder partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v >> 
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 partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v 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));
<<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));
}

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

<<Find parametric representation of cylinder hit>>= 
Float u = phi / phiMax; Float v = (pHit.z - zMin) / (zMax - zMin); <<Compute cylinder partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v >> 
Vector3f dpdu(-phiMax * pHit.y, phiMax * pHit.x, 0); Vector3f dpdv(0, 0, zMax - zMin);
<<Compute cylinder partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v >> 
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 partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v 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:

StartLayout 1st Row 1st Column StartFraction partial-differential normal p Over partial-differential u EndFraction 2nd Column equals left-parenthesis minus phi Subscript normal m normal a normal x Baseline y comma phi Subscript normal m normal a normal x Baseline x comma 0 right-parenthesis 2nd Row 1st Column StartFraction partial-differential normal p Over partial-differential v EndFraction 2nd Column equals left-parenthesis 0 comma 0 comma z Subscript normal m normal a normal x Baseline minus z Subscript normal m normal i normal n Baseline right-parenthesis period EndLayout

<<Compute cylinder partial-differential normal p slash partial-differential u and partial-differential normal p slash partial-differential v >>= 
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

StartLayout 1st Row 1st Column StartFraction partial-differential squared normal p Over partial-differential u squared EndFraction 2nd Column equals minus phi Subscript normal m normal a normal x Superscript 2 Baseline left-parenthesis x comma y comma 0 right-parenthesis 2nd Row 1st Column StartFraction partial-differential squared normal p Over partial-differential u partial-differential v EndFraction 2nd Column equals left-parenthesis 0 comma 0 comma 0 right-parenthesis 3rd Row 1st Column StartFraction partial-differential squared normal p Over partial-differential v squared EndFraction 2nd Column equals left-parenthesis 0 comma 0 comma 0 right-parenthesis period EndLayout

<<Compute cylinder partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v >>= 
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 partial-differential bold n Subscript slash partial-differential u and partial-differential bold n Subscript slash partial-differential v 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 phi 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 z and phi >> 
Point3f pObj = Point3f(radius * std::cos(phi), radius * std::sin(phi), z); <<Reproject pObj to cylinder surface and compute pObjError>> 
Float hitRad = std::sqrt(Sqr(pObj.x) + Sqr(pObj.y)); pObj.x *= radius / hitRad; pObj.y *= radius / hitRad; Vector3f pObjError = gamma(3) * Abs(Vector3f(pObj.x, pObj.y, 0));
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 z and phi , the corresponding object-space position and normal are easily found.

<<Compute cylinder sample position pi and normal n from z and phi >>= 
Point3f pObj = Point3f(radius * std::cos(phi), radius * std::sin(phi), z); <<Reproject pObj to cylinder surface and compute pObjError>> 
Float hitRad = std::sqrt(Sqr(pObj.x) + Sqr(pObj.y)); pObj.x *= radius / hitRad; pObj.y *= radius / hitRad; Vector3f pObjError = gamma(3) * Abs(Vector3f(pObj.x, pObj.y, 0));
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; }