## 3.3 Cylinders

Another useful quadric is the cylinder; pbrt provides cylinder Shapes that are centered around the axis. The implementation is in the files shapes/cylinder.h and shapes/cylinder.cpp. The user supplies a minimum and maximum value for the cylinder, as well as a radius and maximum sweep value (Figure 3.6).

<<Cylinder Declarations>>=
class Cylinder : public Shape { public: <<Cylinder Public Methods>>
Cylinder(const Transform *ObjectToWorld, const Transform *WorldToObject, bool reverseOrientation, Float radius, Float zMin, Float zMax, Float phiMax) : Shape(ObjectToWorld, WorldToObject, reverseOrientation), radius(radius), zMin(std::min(zMin, zMax)), zMax(std::max(zMin, zMax)), phiMax(Radians(Clamp(phiMax, 0, 360))) { } Bounds3f ObjectBound() const; bool Intersect(const Ray &ray, Float *tHit, SurfaceInteraction *isect, bool testAlphaTexture) const; bool IntersectP(const Ray &ray, bool testAlphaTexture) const; Float Area() const; Interaction Sample(const Point2f &u) const;
protected: <<Cylinder Private Data>>
const Float radius, zMin, zMax, phiMax;
};

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

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

<<Cylinder Public Methods>>=
Cylinder(const Transform *ObjectToWorld, const Transform *WorldToObject, bool reverseOrientation, Float radius, Float zMin, Float zMax, Float phiMax) : Shape(ObjectToWorld, WorldToObject, reverseOrientation), radius(radius), zMin(std::min(zMin, zMax)), zMax(std::max(zMin, zMax)), phiMax(Radians(Clamp(phiMax, 0, 360))) { }

<<Cylinder Private Data>>=
const Float radius, zMin, zMax, phiMax;

### 3.3.1 Bounding

As was done with the sphere, the cylinder bounding method computes a conservative bounding box using the range but without taking into account the maximum .

<<Cylinder Method Definitions>>=

### 3.3.2 Intersection Tests

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

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

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

<<Initialize EFloat ray coordinate values>>
EFloat ox(ray.o.x, oErr.x), oy(ray.o.y, oErr.y), oz(ray.o.z, oErr.z); EFloat dx(ray.d.x, dErr.x), dy(ray.d.y, dErr.y), dz(ray.d.z, dErr.z);
EFloat a = dx * dx + dy * dy; EFloat b = 2 * (dx * ox + dy * oy); EFloat c = ox * ox + oy * oy - EFloat(radius) * EFloat(radius);

The solution process for the quadratic equation is similar for all quadric shapes, so some fragments from the Sphere intersection method will be reused in the following.

<<Cylinder Method Definitions>>+=
bool Cylinder::Intersect(const Ray &r, Float *tHit, SurfaceInteraction *isect, bool testAlphaTexture) const { Float phi; Point3f pHit; <<Transform Ray to object space>>
Vector3f oErr, dErr; Ray ray = (*WorldToObject)(r, &oErr, &dErr);
<<Initialize EFloat ray coordinate values>>
EFloat ox(ray.o.x, oErr.x), oy(ray.o.y, oErr.y), oz(ray.o.z, oErr.z); EFloat dx(ray.d.x, dErr.x), dy(ray.d.y, dErr.y), dz(ray.d.z, dErr.z);
EFloat a = dx * dx + dy * dy; EFloat b = 2 * (dx * ox + dy * oy); EFloat c = ox * ox + oy * oy - EFloat(radius) * EFloat(radius);
<<Solve quadratic equation for t values>>
EFloat t0, t1; if (!Quadratic(a, b, c, &t0, &t1)) return false; <<Check quadric shape t0 and t1 for nearest intersection>>
if (t0.UpperBound() > ray.tMax || t1.LowerBound() <= 0) return false; EFloat tShapeHit = t0; if (tShapeHit.LowerBound() <= 0) { tShapeHit = t1; if (tShapeHit.UpperBound() > ray.tMax) return false; }
<<Compute cylinder hit point and >>
pHit = ray((Float)tShapeHit); <<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 false; tShapeHit = t1; if (t1.UpperBound() > ray.tMax) return false; <<Compute cylinder hit point and >>
pHit = ray((Float)tShapeHit); <<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 false; }
<<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); Float F = Dot(dpdu, dpdv); Float G = Dot(dpdv, dpdv); Vector3f N = Normalize(Cross(dpdu, dpdv)); Float e = Dot(N, d2Pduu); Float f = Dot(N, d2Pduv); Float g = Dot(N, d2Pdvv);
<<Compute and from fundamental form coefficients>>
Float invEGF2 = 1 / (E * G - F * F); 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));
<<Initialize SurfaceInteraction from parametric information>>
*isect = (*ObjectToWorld)( SurfaceInteraction(pHit, pError, Point2f(u, v), -ray.d, dpdu, dpdv, dndu, dndv, ray.time, this));
*tHit = (Float)tShapeHit;
return true; }

As with spheres, the implementation here refines the computed intersection point to ameliorate the effect of accumulated rounding error in the point computed by evaluating the ray equation; see Section 3.9.4. We can then 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 = ray((Float)tShapeHit); <<Refine cylinder intersection point>>
phi = std::atan2(pHit.y, pHit.x); if (phi < 0) phi += 2 * Pi;

The next part of 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 tried—this resembles 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 false; tShapeHit = t1; if (t1.UpperBound() > ray.tMax) return false; <<Compute cylinder hit point and >>
pHit = ray((Float)tShapeHit); <<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 false; }

Again the value is computed by scaling to lie between 0 and 1. Straightforward 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); Float F = Dot(dpdu, dpdv); Float G = Dot(dpdv, dpdv); Vector3f N = Normalize(Cross(dpdu, dpdv)); Float e = Dot(N, d2Pduu); Float f = Dot(N, d2Pduv); Float g = Dot(N, d2Pdvv);
<<Compute and from fundamental form coefficients>>
Float invEGF2 = 1 / (E * G - F * F); 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 quite 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); Float F = Dot(dpdu, dpdv); Float G = Dot(dpdv, dpdv); Vector3f N = Normalize(Cross(dpdu, dpdv)); Float e = Dot(N, d2Pduu); Float f = Dot(N, d2Pduv); Float g = Dot(N, d2Pdvv);
<<Compute and from fundamental form coefficients>>
Float invEGF2 = 1 / (E * G - F * F); 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);

### 3.3.3 Surface Area

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

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