Book Corrections for Geometric Tools for Computer Graphics


Book Corrections Organized by Date of Change

28 November 2021, page 301-302. The pseudocode for CircleThroughPointTangentToLineGivenRadius has several typographical errors and code errors. On page 301, the "if (cPrime < 0)" block assigns the negative line direction to (a,b,c). The block also needs to negate cPrime. The declaration of tmp2 has a sum of squares, so it can never be less than -epsilon. It should be "tmp2 = r * r - tmp1 * tmp1". There are several occurrences of "r" in the pseudocode. These should be "radius". On page 302, the last two assignments of the center values should be "center[0] = tmpPt + tmp2 * D" and "center[1] = tmpPt - tmp2 * D" where "Point2D D = { b, -a };"

17 December 2019, page 272. In the GetExtremeIndex function, the input direction D is listed as type Point but should be type Vector. The call to GetMiddleIndex has only two parameters but should have three parameters: mid = GetMiddleIndex(i0,i1,C.N);

18 November 2019, page 213. In the next to last paragraph, the equation for a solid cone is given as Dot(hat(a),X-V) >= cos(theta) but should be Dot(hat(a),X-V) >= cos(theta)*Length(X-V).

19 May 2018, pages 661-662. Yet another badly written topic.
  • Throughout the discussion, the torus point X is written as X = [Xx   Xy   Xz]. This is awkward notation when the torus point is simply (x,y,z).
  • On page 661, "cutting with the XZ plane" should be "cutting with the XY plane".
  • On page 661, the displayed equation for u has two occurrences of "arccos(θ)" that should be "arccos(Xx/ru)".
  • On page 661, "...if we want the parametric origin to be on the 'inside' of the torus...". In more mathematical terms: Given a circular cross section of radius r for the torus, we want the circle point closest to the origin to correspond to parameter v = 0.
  • On page 662, it is not clear in Figure 11.77 which angle is labeled with φ. Turns out that the right triangle has φ corresponding to the circle center. This angle starts at 0 for the circle point at distance R-r from the origin. It increase to π as you walk the circle counterclockwise to the circle point at distance R+r from the origin.
  • On page 662, the length rv is the same as length ru. Perhaps the distance from the origin to the projected point should have been named something simple like d.
  • On page 662, the displayed equation has "sin(φ) = Xz/rv" but should be "sin(φ) = Xz/r".
  • On page 662, the displayed equation specifies parameter u but should be parameter v. There are two occurrences of "arccos(φ)" that should be "arccos(-(rv-R)/r)".
  • The paragraphs starting with "The partial derivatives..." are mathematical gibberish. After a discussion of computing functions u(x,y,z) and v(x,y,z), one might assume that the partial derivatives of u and v with respect to x, y and z were to be computed. This is not the case. The notation ∂u and ∂v is not standard. The vector (X'-O) lives in 3D, so what does (X'-O) mean? What does it mean "the partial derivatives are computed in the local space of the torus"? The intent of the paragraphs is to motivate computing tangent vectors to the torus, presumably as a parametric surface P(u,v) =(x(u,v),y(u,v),z(u,v)). Linearly independent tangent vectors are ∂P/∂u and ∂P/∂v. Computing them requires knowing how x, y and z depend on u and v, but the previous discussion shows only how to compute u and v as functions of x, y and z.
A simpler presentation is the following.

A torus with origin (0,0,0), outer radius r0 and inner radius r1 (with r0 >= r1) is defined implicitly as follows. The point P0 = (x,y,z) is on the torus. Its projection onto the xy-plane is P1 = (x,y,0). The circular cross section of the torus that contains the projection has radius r0 and center P2 = r0*(x,y,0)/sqrt(x2+y2). The points triangle <P0,P1,P2> is a right triangle with right angle at P1. The hypotenuse <P0,P2> has length r1, leg <P1,P2> has length z and leg <P0,P1> has length |r0 - sqrt(x^2+y^2)|. The Pythagorean theorem says z2 + |r0 - sqrt(x2+y2)|2 = r12. This can be algebraically manipulated to (x2 + y2 + z2 + r02 - r12)2 - 4 * r02 * (x2 + y2) = 0 A parametric form is x = (r0 + r1 * cos(v)) * cos(u), y = (r0 + r1 * cos(v)) * sin(u), z = r1 * sin(v) for u in [0,2*pi) and v in [0,2*pi). The position is P = (x,y,z) and the partial derivatives are ∂P/∂u = (r0 + r1 * cos(v)) * (-sin(u), cos(u), 0) and ∂P/∂v = r1(-cos(u) * sin(v), -sin(u) * sin(v), cos(v)). It is simple to show that the partial derivatives are perpendicular. To compute (u,v) from (x,y,z), u = atan2(y,x) and v = atan2(z,sqrt(x2+y2)-r0).

31 March 2018, pages 442-446. Section 10.9.2 for computing the distance between a line and a rectangle in 3D is filled with errors--what an atrocious description. As a replacement, please use Distance from Line to Rectangle in 3D, which describes a robust and efficient algorithm for computing the distance.
  • On page 442, Figure 10.42 is misleading because the line is drawn as if it is neither parallel nor perpendicular to the rectangle, yet the closest point Q is shown to be on the rectangle. This is not possible for a non-parallel/non-perpendicular line.
  • On page 442, the displayed |R(u,v) - L(t)| is missing a power 2.
  • On page 442, in the second displayed equation involving Q(u,v,t) it would have been helpful to use Δ for V-P in order to simplify reading the equations. Of course, things are worse in that the equation is wrong. It should be
    (e0 • e0)u2 + (e1 • e1)v2 + (d • d)t2 + 2(e0 • e1)u v - 2(e0 • d)u t - 2(e1 • d)v t
     + 2(e0 • Δ)u + 2(e1 • Δ)v - 2(d • Δ)t + |Δ|2
    The definitions for the quadratic coefficients on page 443 must be adjusted accordingly.
    a00 = (e0 • e0), a01 = 2(e0 • e1), a11 = (e1 • e1), a02 = -2(e0 • d), a12 = -2(e1 • d),
    a22 = (d • d), b0 = 2(e0 • Δ), b1 = 2(e1 • Δ), b2 = -2(d • Δ), c = |Δ|2 
    The book's a12 is listed twice; the first occurrence is wrong. The book's b2 is also wrong. The book also failed to define c in that list.
  • On page 443, "...then the solution (bar(u),bar(v),bar(t))..." does not say what problem is being solved. It is the solution to the linear system of equations mentioned later. Perhaps "...then the minimizer (bar(u),bar(v),bar(t))..." is better said.
  • On page 443, there is an expression |T(u,v,t)-L(t)| that should be |R(u,v)-L(t)|.
  • On page 443, "Otherwise, the minimum of ∇ Q = (bar(u),bar(v),bar(t))..." is incorrect terminology. The minimum we seek is for Q. For a convex quadratic with unconstrained domain, the global minimum occurs when ∇ Q is (0,0,0), which leads to a minimizer (bar(u),bar(v),bar(t)). In the constrained domain, the minimum can occur on a face of the convex domain, but any gradient is computed in the domain of that face. Better is "Otherwise, the minimizer (bar(u),bar(v),bar(t)) of Q will be on a face...".
  • Speaking of that linear system at the bottom of page 443, it is incorrectly written. The variables are written on the right-hand side and the b-vector on the left. These are swapped, the b-vector needs to be negated and the matrix of coefficients should be written to show the symmetry.
  • On page 444, a01, a02, and a12> each needs an additional factor of 2. On page 445, b0 and b1 each needs an additional factor of, and b2 needs an additional factor of -2.
  • On page 445, The block of pseudocode that flips signs when the determinant is negative does not make sense. Theoretically, the determinant is nonnegative, and it can be computed with an algorithm different from that of the pseudocode to guarantee that it is nonnegative even with floating-point rounding errors.
  • On page 445, the test for "(near) parallelism" is one way to try to avoid floating-point rounding errors when computing the distance. But as the PDF mentioned previously points out, subtractive cancellation can cause gross errors in the result. I recommend using the constrained conjugate gradient algorithm of that PDF.
  • On page 445, there are terms "rectangle[i]" for i in {0,1,2,3}. Previously in the pseudocode, the notation was "rectangle.V[i]".
  • On page 446, the pseudocode first computes u and v but then the conditional tests involves s and t. The s and t are meant to be u and v, respectively.

30 August 2017, page 771. In the SegmentInCode listing, the case for the "vertex is convex" has the inequalities reversed. The line of code should be
                      return Kross(diff, edgeR) < 0 and Kross(diff, edgeL) > 0
                  

02 July 2016, page 941. All the product formulas are incorrect. They should be
                      sin(a)*sin(b) = (cos(a-b) - cos(a+b))/2
                      sin(a)*cos(b) = (sin(a+b) + sin(a-b))/2
                      cos(a)*cos(b) = (cos(a-b) + cos(a+b))/2
                  
The last formula in the book appears to be an attempt to cover all the cases, but the first three equations are all the cases.

07 July 2014, pages 157-158. For simplicity of notation here, define N = hat(n), N^T is the transpose of N, and V = hat(v). The last row, first column entry of the general shearing matrix should be -(tan(theta))*Dot(N,Q)*V. The first line on page 158 makes an appeal to applying the point transformation T to a vector N, which makes no sense mathematically. Even if it did, it is not clear why the author claims T causes N to change at all. Let P = P1 + h*N, where h = Dot(N,P-Q). The shearing is P" = P + t*V, where from trigonometry, t = h*tan(theta), so P" = P + (tan(theta))*Dot(N,P-Q)*V and P"-Q = [I + (tan(theta))*V*N^T](P-Q). Add Q to both sides and factor to P" = [I + (tan(theta))*V*N^T]*P - (tan(theta))*V*N^T*Q = [I + (tan(theta))*V*N^T]*P - (tan(theta))*Dot(N,Q)*V.

14 December 2013, pages 208-210. The algorithm uses barycentric coordinates to determine closest edge, but this does not work for a triangle with an obtuse angle. Thus, the SquaredDistance pseudocode is incorrect. The online source code does not use this algorithm; it uses the previously discussed (correct) algorithm.

07 October 2013, page 229. Equation (6.24) has 5 cases (top to bottom). The final inequality in case 2 should be Dot(d0,delta+T0*d0) ≤ 0. The final inequality in case 4 should be Dot(d0,delta+T0*d0-T1*d1) ≤ 0.

21 September 2013, page 15. The "Matrix" row of Table 2.1 has pseudocode "Matrix3x3 m, m;" that should be "Matrix3x3 m, n;".

01 March 2013, page 225. The denominator of Equation (6.17) should be dot(perp(d0),d1) rather than dot(perp(d1),d0); this has the effect of negating the right-hand side of the equation.

03 November 2012, page 78. Figure 3.19 is confusing and has one typographical error. It is an attempt to describe the linear transformation T(x) = 2*x. In terms of basis vectors u and v: T(u) = 2*u and T(v) = 2*v. The vector s in the figure should be more clearly written:
  s = T(2*v)
    = 2*T(v)   // because T is linear
    = 2*(2*v)  // by definition of T
Call the top vector of the figure r, in which case
  r = T(3*u)
    = 3*T(u)   // because T is linear
    = 3*(2*u)  // by definition of T
Thus, the current figure should have "3T(u) = 2(3u)".

03 November 2012, page 79. Figure 3.20 is confusing and has two typographical errors. It is an attempt to describe the linear transformation that is defined on the basis vectors u and v: T(u) = 2*u and T(v) = 1.5*v. Given a vector x = a*u+b*v, the transformation is T(x) = T(a*u+b*v) = a*T(u)+b*T(v) = a*(2*u)+b*(1.5*v). The vector s in the figure is incorrect and should be
  s = T(2*v)
    = 2*T(v)   // because T is linear
    = 2*(1.5*v)  // by definition of T
Call the top vector of the figure r, in which case
  r = T(3*u)
    = 3*T(u)   // because T is linear
    = 3*(2*u)  // by definition of T
Thus, the current figure should have "3T(u) = 2(3u)".

10 September 2012, page 511. Lines 4 and 5 are missing right parentheses. The pseudocode should be
  t[0] = (-b + sqrt(discrm))/(2*a);
  t[1] = (-b - sqrt(discrm))/(2*a);

10 September 2012, page 518. The page has 7 references to array elements "[3]". These should be "[2]".

10 September 2012, page 511. Lines 4 and 5 are missing right parentheses. The pseudocode should be
  t[0] = (-b + sqrt(discrm))/(2*a);
  t[1] = (-b - sqrt(discrm))/(2*a);

08 July 2012, page 548. In matrix M, entry 6 of row 2 says "2(ux vy + vx vy)". This should be "2(ux vy + vx uy)".

24 June 2012, page 83. The sentence starting with "We've just done two forbidden operations on points--scaling by a vector and directly adding them." This should say "We've just done two forbidden operations on points--multiplying points by scalars and directly adding them."

24 June 2012, page 84. The first displayed equation has a term "alpha3*(P3-P2)". This term should be "alpha3*(P3-P1)".

20 February 2010, pages 564-569. The first paragraph of Intersection Detection mentions a base point B and refers you to Figure 11.14, which shows the base point labeled as C. Later in Intersection Detection, there is a section Finite Cone that refers to the base point using name B. This point is as shown in Figure 11.14 and is not the same as the apex A.

On page 567, a vector hat(w) is defined. The author used the hat notation to denote unit-length vectors. However, the vector Cross(hat(d),Cross(hat(n),hat(d))) is not necessarily unit length. For example, if hat(n) is nearly parallel to hat(d), Cross(hat(n),hat(d)) is nearly the zero vector, in which case hat(w) is nearly the zero vector. The correct equation is
    hat(w) = Cross(hat(d),Cross(hat(n),hat(d)))/Length(Cross(hat(n),hat(d)))
On pages 567, 568, and 569, the author has occurrences of arrow(w), but these should be hat(w).

The variable a is the distance between points B and Ia, so the last equation on page 567 should be
    a = Length(Ia - B)
That is, the subtraction of h is incorrect. Just to be clear, the variable b is the distance between points B and Ic, and the variable c is the distance between points Ia and Ic. On pages 568 and 569, any occurrence of Length(Ia-B)-h should be replaced by Length(Ia-B).

Observe that with the correct hat(w) equation, it must be that Dot(hat(n),hat(w)) = Length(Cross(hat(n),hat(d)), so it is not necessary to explicitly compute hat(w) in the source code. The equation for b2 is
  b2 = a2 [1/Length(Cross(hat(n),hat(d)))2 - 1]
However, notice that this equation has numerical problems when hat(n) and hat(d) are nearly parallel. The length of the cross product is nearly zero, so computing the reciprocal length is a problem. Instead, you want to code up the test as
  v = Cross(hat(n),hat(d));
  lensqr = Dot(v,v);
  a*a*(1 - lensqr) ≤ lensqr*r*r;
Here is an alternative for the intersection test of plane and finite cone. Let the plane have unit-length normal N and contain point P. The plane equation is Dot(N,X-P) = 0. The idea is to compute the interval of projection of the finite cone onto the line containing P with direction N. This inteval has endpoints P + t0*N and P + t1*N, where t0 < t1. If the interval [t0,t1] contains 0, then the cone and plane intersect.

The projection interval endpoints are determined by one of the following.
  1. C and two points on the circle that terminates the cone wall. The two points project to the same interval endpoint. This happens when N and A are parallel.

  2. Two points on the circle.

  3. C and one point on the circle.
The cone has vertex C, unit-length axis direction D, height h, and angle A measured from the axis to the cone wall. The circle that terminates the cone wall is parameterized by
    X = C + h*D + h*tan(A)*(cos(B)*U + sin(B)*V)
for B in [0,2*pi). The vectors U and V are unit length and {U,V,D} is an orthonormal set (vectors are unit length and mutually perpendicular). The projection of X onto the normal line is P + t(B)*N, where
   t(B) = Dot(N,X-P) = Dot(N,C-P) + h*Dot(N,D) + h*tan(A)*(cos(B)*Dot(N,U) + sin(B)*Dot(N,V))
I have indicated that t is a function of B (different circle points project to different normal line points). The projection of C onto the normal line is Dot(N,C-P).

The minimum value of t(B) is t0 and the maximum value of t(B) is t1. Define d = Dot(N,C-P) and
 r = sqrt(Dot(N,U)^2 + Dot(N,V)^2) = sqrt(1 - Dot(N,D)^2)
t0 is the minimum, and t0 is the maximum, of the set
    {d, d + h*Dot(N,D) - h*tan(A)*r, d + h*Dot(N,D) + h*tan(A)*r}
The pseudocode is
    bool ConeIntersectsPlane (Cone cone, Plane plane)
    {
        Vector CmP = cone.C - plane.P;
        float d = Dot(plane.N, CmP);
        float NdotD = Dot(plane.N, cone.D);
        float r = sqrt(1 - NdotD*NdotD);
        float tmp0 = d + cone.h*NdotD;
        float tmp1 = cone.h*cone.tanAngle*r;
        float vm = tmp0 - tmp1;
        float vp = tmp0 + tmp1;
        float t0 = min(d, min(vm, vp));
        float t1 = max(d, max(vm, vp));
        return t0 <= 0 && 0 <= t1;
    }
This algorithm does not require you to explicitly determine which of the book-mentioned cases 1, 2, 3a, or 3b applies.

03 June 2009, page 454. The pseudocode after the "Closest point on the box" comment should be
  qPrime.x = box.extents.x;
  qPrime.y = line.origin.y;
  qPrime.z = line.origin.z;

29 January 2009, page 43. Bullet 'v' in Section 2.6.3 is a repeat of bullet 'ii' and should be removed.

29 January 2009, page 44. The example combination of vectors in a displayed equation is
  2v1 + 3v3 +-1v3
but should be
  2v1 + 3v2 +-1v3
That is, the subscript on the middle term needed to be fixed.

29 January 2009, page 46. The sentence: "We can map real numbers to real numbers: ..." should say "We can map nonnegative real numbers fo nonnegative real numbers: ..."

29 January 2009, page 53. The last sentence should say "... into Equations 2.4 and 2.5, ...".

29 January 2009, page 777. The first sentence of the last paragraph has a subject-verb mismatch. It should be "An edge sk has ...". In the same paragraph, the end of the second line has "s1.y" but should be "s1.y0".

06 December 2008, page 872. The first displayed equation is
  max{|a_0|, 1+|a_1|, ..., 1+|a_{n-1}|} = 1 + max{|a_0|, ..., |a_{n-1}|}
but the last equaltity should be an inequality
  max{|a_0|, 1+|a_1|, ..., 1+|a_{n-1}|} ≤ 1 + max{|a_0|, ..., |a_{n-1}|}

14 November 2008, page 428. The paragraph about the case of region 2 has "the gradient at the corner (1,1) cannot point inside the unit square". This should say "outside" rather than "inside".

17 August 2008, page 562. The pseudocode has a few issues. The first line of code that computes d should be a signed distance calculation. The first comparison of cos(theta) should have the abs about the cosine term, not about the epsilon term. If an epsilon comparison is made for when cos(theta) is nearly zero, to be consistent an epsilon comparison should also be made for when cos(theta) is nearly one (the circle case). Here is my own rewrite of the pseudocode without epsilon tests. The following data structures are assumed.
    struct Plane    { Point origin; Vector normal; };
    struct Cylinder { Point origin; Vector direction; float radius; };
    struct Line     { Point origin; Vector direction; };
    struct Circle   { Point center; Vector normal; float radius; };
    struct Ellipse  { Point center; Vector normal, major, minor; float majorLength, minorLength; };
    
    enum Type
    {
        EMPTY_SET,
        ONE_LINE,
        TWO_LINES,
        CIRCLE,
        ELLIPSE
    };
    struct Intersection
    {
        Type type;
        Line line;          // valid when type = ONE_LINE
        Line line1, line2;  // valid when type = TWO_LINES
        Circle circle;      // valid when type = CIRCLE
        Ellipse ellipse;    // valid when type = ELLIPSE
    };
The plane normal, cylinder direction, line direction, circle normal, ellipse normal, ellipse major direction, and ellipse minor direction are all required to be unit-length vectors.

The output is intr. The returned bool value is true when there is an intersection, in which case intr is valid (of the type specified in its Type member). The returned value is false when there is no intersection, in which case the Type member has value EMPTY_SET.
    bool PlaneCylinderIntersection (Plane plane, Cylinder cylinder, Intersection& intr)
    {
        float d = Dot(plane.normal, cylinder.origin - plane.origin);
        Point projectedCylinderOrigin = cylinder.origin - d * plane.normal;
        float cosTheta = Dot(cylinder.direction, plane.normal);
        if (abs(cosTheta) > 0)  // cylinder axis intersects plane in a unique point
        {
            if (abs(cosTheta) < 1)  // intersection is an ellipse
            {
                intr.type = ELLIPSE;
                intr.ellipse.normal = plane.normal;
                intr.ellipse.center = projectedCylinderOrigin - (d/cosTheta) * cylinder.direction;
                intr.ellipse.major = cylinder.direction - cosTheta * plane.normal;
                intr.ellipse.minor = Cross(plane.normal, intr.ellipse.major);
                intr.ellipse.majorLength = cylinder.radius / abs(cosTheta);
                intr.ellipse.minorLength = cylinder.radius;
                Normalize(intr.ellipse.major);
                Normalize(intr.ellipse.minor);
                return true;
            }
            else  // abs(cosTheta) == 1,  intersection is a circle
            {
                intr.type = CIRCLE;
                intr.circle.normal = plane.normal;
                intr.circle.center = projectedCylinderOrigin;
                intr.circle.radius = cylinder.radius;
                return true;
            }
        }
        else  // abs(cosTheta) == 0,  cylinder is parallel to plane
        {
            if (abs(d) < cylinder.radius)  // cylinder intersects plane in two lines
            {
                Vector offset = Cross(cylinder.direction, plane.normal);
                float e = sqrt(cylinder.radius * cylinder.radius - d * d);
                intr.type = TWO_LINES;
                intr.line1.origin = projectedCylinderOrigin - e * offset;
                intr.line1.direction = cylinder.direction;
                intr.line2.origin = projectedCylinderOrigin + e * offset;
                intr.line2.direction = cylinder.direction;
                return true;
            }
            else if (abs(d) == cylinder.radius)  // cylinder intersects plane in one line
            {
                intr.type = ONE_LINE;
                intr.line.origin = projectedCylinderOrigin;
                intr.line.direction = cylinder.direction;
                return true;
            }
            else // abs(d) > cylinder.radius,  cylinder does not intersect plane
            {
                intr.type = EMPTY_SET;
                return false;
            }
        }
    }

17 August 2008, page 550. The pseudocode computes b as an absolute value, but the absolute value should only be used in the comparison to the sphere radius on the line after computing b. The pseudocode on this page has had other corrections, so the revision is listed next. I do not like the book's formulation--the revision is more like how I write the actual code. This is the mathematically correct formulation, without concern for "epsilon" arithmetic to deal with numerical round-off errors. The plane normal is assumed to be a unit-length vector. The Circle3 data structure is for a circle in three dimensions.
    bool PlaneSphereIntersection (Plane plane, Sphere sphere, Circle3& circle)
    {
        float b = dotProd(plane.normal, sphere.center - plane.pointOnPlane);
        float absB = fabs(b);
        circle.center = sphere.center - b * plane.normal;
        circle.normal = plane.normal;
        if (absB <= sphere.radius)
        {
            // The sphere intersects the plane in a circle.  The circle is
            // degenerate when fabs(b) is equal to sphere.radius, in which
            // case the circle radius is zero.
            circle.radius = sqrt(sphere.radius * sphere.radius - b * b);
            return true;
        }
        return false;
    }
The circle parameter is an output of the function. It is valid only when the function's return value is true. The plane's normal vector is copied to the circle's normal so that you can actually reconstruct the circle. If the circle center is K, the normal is N, and the radius is R, then choose vectors U and V so that {U,V,N} is an orthonormal set (vectors are unit length and mutually perpendicular). A circle parameterization is
  X(t) = K + R*(cos(t)*U + sin(t)*V),  0 <= t < 2*pi
A numerically robust way to compute U and V is
    if (fabs(N.x) >= fabs(N.y))
    {
        // N.x or N.z is the largest magnitude component, swap them.
        invLength = 1/sqrt(N.x * N.x + N.z * N.z);
        U = Vector(-N.z * invLength, 0, N.x * invLength);
        V = Vector(N.y * U.z, N.z * U.x - N.x * U.z, -N.y * U.x);  // = Cross(N,U)
    }
    else
    {
        // N.y or N.z is the largest magnitude component, swap them.
        invLength = 1/sqrt(N.y * N.y + N.z * N.z);
        U = Vector(0, N.z * invLength, -N.y * invLength);
        V = Vector(N.y * U.z - N.z * U.y, -N.x * U.z, N.x * U.y);  // = Cross(N,U)
    }
If you want to allow for numerical round-off error, say, by requiring the intersection to be a single point when b is within epsilon of the radius, use
    bool PlaneSphereIntersection (Plane plane, Sphere sphere, Circle3& circle)
    {
        float b = dotProd(plane.normal, sphere.center - plane.pointOnPlane);
        float absB = fabs(b);
        circle.center = sphere.center - b * plane.normal;
        circle.normal = plane.normal;
        if (absB <= sphere.radius - epsilon)
        {
            // The sphere intersects the plane in a circle.  The circle is
            // degenerate when fabs(b) is equal to sphere.radius, in which
            // case the circle radius is zero.
            circle.radius = sqrt(sphere.radius * sphere.radius - b * b);
            return true;
        }
        else if (absB <= sphere.radius + epsilon)
        {
            // The sphere is _nearly_ tangent to the plane.  Choose the
            // intersection to be the point of tangency.
            circle.radius = 0;
            return true;
        }
        return false;
    }
The magnitude of the nonnegative number epsilon is up to the programmer. Observe that if you choose epsilon to be zero, you effectively have the pseudocode for the mathematically correct case.

14 August 2008, page 571. The end of the last sentence of the first paragraph of the subsection Nondegenerate Plane-Cone Intersections should be: cos(theta) = Dot(hat(a),hat(n)).

14 August 2008, page 575. The style of the pseudocode differs from other pages in this section. In the first block starting with ``float d = '', remove the 'float' types (to be consistent). In the section about the parabola, consistent notation would be
    parabola.axis = Normalize(cone.axis - cosTheta * plane.normal);
    parabola.perpendicularAxis = Cross(parabola.axis, plane.normal);
    parabola.vertex = cone.vertex + d* plane.normal + e*parabola.axis;
    parabola.focalLength = d/2 * cosTheta;
Also on this page, the first paragraph has a condition on the dot product of the plane normal hat(n) and the cone axis hat(a) in order that the two vectors be parallel. The condition should be: |Dot(hat(n),hat(a)| = 1 (for the true mathematical classification) or |Dot(hat(n),hat(a))| >= 1-epsilon (if you want to worry about numerical round-off errors).

14 August 2008, page 580. The pseudocode has a line that computes the minor axis. It should be
    minorAxis = Cross(plane.normal, majorAxis);
because nowhere is there defined ``ellipse.u''.

28 October 2007, page 888. The last displayed equation is
    E(hat(c)) = (sum_{i=0}^n hat(c) * vector(v))^2 = hat(c)^T M hat(c)
which is the square of a sum, but it should be a sum of squares,
    E(hat(c)) = sum_{i=0}^n (hat(c) * vector(v))^2 = hat(c)^T M hat(c)

28 October 2007, page 889. Line 5 of the first paragraph of section A.7.8 has
    M = sum_{i=0}^2 vector(v) vector(v)^T
The upper limit should be n,
    M = sum_{i=0}^n vector(v) vector(v)^T

22 September 2007, page 550. The last line of pseudocode on this page should be
    float radius = sqrt(sphere.radius^2 - b^2);
The pseudocode also rejects the "just touching" case as "not an intersection"; that is, "b < sphere.radius" should be "b ≤ sphere.radius". In fact, since the pseudocode allows a small amount of error ("radius < epsilon"), it might even be good to support an error term in the first test: "if (b ≤ sphere.radius + epsilon)".

05 February 2007, page 228. The symbols in Equation (6.23) on page 228 are not typeset correctly. In that equation, let us say that the first row of the equation is region R0 and the last row of the equation is R8 (the parameter plane is partitioned into 9 regions). The corrections are
    R2:  0 < tilde(t0) < T0 and bar(t1) >= T1
    R4:  0 < tilde(t1) < T1 and bar(t0) >= T0
    R6:  hat(t0) >= T0 and tilde(t1) <= 0
    R7:  tilde(t0) <= 0 and hat(t1) >= T1
    R8:  tilde(t0) >= T0 and tilde(t1) >= T1

28 January 2007, page 633. The pseudocode for the ray-OBB intersection is absurd. The first line of code is missing two right parentheses after "box.axis[i]". The inequality "-r + box.halfLength[i] > 0" should be "-r + box.halfLength[i] < 0". The rest of the pseudocode is worthless. The value s can very well be zero, leading to a division by zero when computing t0 and t1. If you find worthless pseudocode in the book, first let me know and, second, look at the code at my website. I wrote the real code (not the pseudocode), so if you find errors in the real code, please let me know.

23 January 2007, page 612. The last line of pseudocode is
    if (value < min) min = value; else max = value;
This should be
    if (value < min) min = value; else if (value > max) max = value;

31 October 2006, page 550. To be consistent with the other pseudocode, "if (b < r)" should be written as "if (b < sphere.radius)".

28 October 2006, page 548. Equation (11.20) has "r" when it should be "r^2" (r squared).

08 September 2006, page 384. The book-correction entry for page 635 (24 February 2004) was
  nx*(x1+x0)+ny*(y1+y0)+nz*(z1+z0)+2*d <= |nx*(x1-x0)|+|ny*(y1-y0)|+|nx*(z1-z0)|
but the left-hand side is missing absolute values. It should be
  |nx*(x1+x0)+ny*(y1+y0)+nz*(z1+z0)+2*d| <= |nx*(x1-x0)|+|ny*(y1-y0)|+|nx*(z1-z0)|
I have modified the older entries here.

06 September 2006, page 384. The first displayed equation has a Q that should be a P:
    Q' = P + s*e0 + t*e1

30 March 2006, page 387. The top line of pseudocode has "Dot(p,n)-d", but should be "Dot(p,n)+d".

28 January 2006, pages 417-426. These pages cover distance calculations between linear components. The pseudocode has major problems. Apparently, one of these blocks of code was written and then cut-paste-modified for the other blocks. They all exhibit the same ills. The difference between closest points needs to be of the form
  v = (base0 + s*direction0) - (base1 + t*direction1)
That is, the book parentheses are incorrect. Worse is the handling of the distinction between parallel and nonparallel components. If "det" is ever zero, the pseudocode will do a division by it. Also, the epsilon test is problematic when the segment direction vectors have large length. If D0 and D1 are the line segment lengths, each computed as a difference of the line segment end points, and if A is the angle between the directions, it is the case that
  det = |Cross(D0,D1)|^2 = |D0|^2 * |D1|^2 * (sin(A))^2
If the angle A is very small, but the lengths of the directions are very large, "det" can be quite large, causing you to classify the segments are nonparallel when they are nearly parallel.

My own implementation uses a line segment represention of C+s*D, where C is the midpoint of the segment, D is a unit-length direction, and |s| <= e where e is half the length of the segment. This leads to a test
  det = (sin(A))^2
so that the epsilon test is a relative error measurement, not an absolute error measurement. Section 10.8.3 of the book describes how I think of distance measurements for linear components, and I much prefer this approach based on actually writing code and seeing how it performs on real data.

28 January 2006, page 686. The pseudocode has two sets named CoPos and CoNeg. These should be CoSame and CoDiff, respectively.

13 January 2006, pages 767-768. The visibility graph does not contain an arc (i,i), where i is an index of a polygon vertex. The text should also state that adjacent polygon vertices are not considered to be visible to each other. Therefore, the visibility graph does not contain arcs (i,(i+1) modulo n), where n is the number of polygon vertices and 0 <= i <= n-1. On page 768, the text states that the adjacency matrix for a convex polygon has 3 zero diagonals, the rest of the entries are ones. However, the entries at (row,column) of (0,n-1) and (n-1,0) are also zero since these correspond to adjacent polygon vertices.

13 December 2005, page 120. In the middle of the page, I do not understand why the statement "with, as before, Vector(u) retaining its usual matrix representation" is made. This seems to imply that Vector(u) will be represented as a row matrix as shown at the top of the page. When Cross(Vector(v),Vector(u)) is represented using a skew-symmetric matrix for Vector(v), Vector(u) must be written as a column matrix. Moreover, the skew-symmetric matrix Tilde(v) is the negative of what it should be. Listing the rows from top to bottom, it should be Tilde(v) = ((0,-v3,+v2),(+v3,0,-v1),(-v2,+v1,0)).

22 November 2005, page 419. The pseudocode for c should be
  c = Dot(ray.direction, ray.direction);

22 November 2005, page 531. The pseudocode for b at the bottom of the page should be
  b = (s1 * n1n2dot - s2 * n1normsqr) / (n1n2dot^2 - n1normsqr * n2normsqr);
That is, the first occurrence of n2normsqr in the book should have been n1normsqr.

07 November 2005, page 356. The first equation for the torus is incorrect. It should be the following equation. I have factored it in a manner that leads to more efficient numerical computation:
  (x^2 + y^2 + z^2)^2 - 2*(r0^2 + r1^2)*(x^2 + y^2 + z^2) + 4*r0^2*z^2 + (r0^2 - r1^2)^2 = 0
My torus class (in Wm3Torus files) has the correct formula listed.

07 November 2005, page 660. The equation for c0 is missing a power 2. It should be
    c0 = (Dot(P,P) - (R^2+r^2))^2 - 4*R^2*(r^2 - P_z^2)

02 October 2005, pages 628-629. The pseudocode for ray-AABB intersection has a loop that illustrates clipping only against the planes x=xmin and x=xmax. This makes it unclear how the the full set of tests fit in with the rest of the loop logic. The following pseudocode is the replacement. The Point and Vector data structures are assumed to have operator[] functions to access the coordinates. Index 0 corresponds to x (P[0] and P.x are equivalent), index 1 corresponds to y, and index 2 corresponds to z. It is my opinion that this function should also return the interval of intersection (first point of intersection and last point of intersection). This change shows up in the last input to the function. It is guaranteed that tIntersect[0] <= tIntersect[1]. If the ray just grazes the AABB, we will return the same value in both parameters.
  boolean RayIntersectAABB (Point P, Vector d, AABB box, float tIntersect[2])
  {
      // The first part of this code tests for line-AABB intersection.
      tNear = -INFINITY;
      tFar = INFINITY;
      for (i = 0; i < 3; i++)
      {
          // Check if the line is parallel to the min[i] and max[i] planes.
          if (abs(d[i]) < epsilon)
          {
              // The line is parallel to the min[i] and max[i] planes.
              if (P[i] < box.min[i] || P[i] > box.max[i])
              {
                  // The line is outside the slab defined by min[i] and max[i].
                  return false;
              }
              // else:  The line is inside the slab defined by min[i] and max[i].
              // Just continue testing with the next slab.
              continue;
          }
          // The line is not parallel to the slab planes, so compute the
          // parameters of intersection.  The line segment of intersection is
          // P+t*d, where t0 <= t <= t1.
          t0 = (box.min[i] - P[i])/d[i];
          t1 = (box.max[i] - P[i])/d[i];
          if (t0 > t1)
          {
              tmp = t1;
              t1 = t0;
              t0 = tmp;
          }
          
          // Compare with current values.  The current slab may have increased
          // tNear and/or decreased tFar.
          if (t0 > tNear)
          {
              tNear = t0;
          }
          if (t1 < tFar)
          {
              tFar = t1;
          }
          
          // Check if the line misses the AABB entirely.
          if (tNear > tFar)
          {
              return false;
          }
      }
      
      // At this place in the code, the line intersects the AABB in the line
      // segment P+t*d, tNear <= t <= tFar.  If all you care about is line-AABB
      // intersection, then use the following three lines of code and omit the
      // remaining pseudocode.
      //
      // tIntersect[0] = tNear;
      // tIntersect[1] = tFar;
      // return true;
      // The ray is (P+t*d, t >= 0).  We need to compute the intersection of
      // the interval tNear <= t <= tFar with the interval t >= 0.
      if (tFar < 0)
      {
          return false;
      }
      if (tNear < 0)
      {
          tNear = 0;
      }
      tIntersect[0] = tNear;
      tIntersect[1] = tFar;
      return true;
  }

28 June 2005, page 516. The last line of pseudocode has Matrix4x4 when it should be Matrix3x3.

17 June 2005, page 861. A bunch of errors are in "Quaternion to Matrix" and "Matrix to Quaternion". This material was from "3D Game Engine Design", where the same errors occur, but I had posted the correction for that book in February 2001. If R = [rij], then equation (A.8) should have
  r00 = 1-2y^2-2z^2
  r01 = 2xy-2wz
  r02 = 2xz+2wy
  r10 = 2xy+2wz
  r11 = 1-2x^2-2z^2
  r12 = 2yz-2wx
  r20 = 2xz-2wy
  r21 = 2yz+2wx
  r22 = 1-2x^2-2y^2
The last sentence of the first paragraph of "Matrix to Quaternion" has
  y = (r20-r02)/(4w)
but it should be
  y = (r02-r20)/(4w)
In the second paragraph of that section. When r00 is the maximum, the book has
  w = (r12-r21)/(4x)
but it should be
  w = (r21-r12)/(4x)
When r11 is the maximum, the book has
  w = (r20-r02)/(4y)
but it should be
  w = (r02-r20)/(4y)
When r22 is the maximum, the book has
  w = (r01-r10)/(4z)
but it should be
  w = (r10-r01)/(4z)

16 May 2005, page 377. The equation at the bottom of the page for bar(t) has the wrong numerator. The equation should be "bar(t) = (b*d-a*e)/(a*c-b*b)". The pseudocode on page 380 does use the correct equation.

15 December 2004, pages 844 and 846. The second line from the bottom on page 844 should have "e2 = -98". The fifth line from the bottom of page 846 should have "(1.11596,1.14945,1.63258)".

14 December 2004, page 437. The displayed matrix equation is incorrect. Setting the gradient to zero leads to the system of equations
    a00*u + a01*v + a02*t + b0 = 0
    a01*u + a11*v + a12*t + b1 = 0
    a02*u + a12*v + a22*t + b2 = 0
The pseudocode on that page essentially solves this system using Cramer's rule (computes cofactors and determinant).

13 October 2004, page 534. The equations for Cramer's rule have errors. The inverse determinant should be
  invDet = 1/(a0*BC - b0*AC + c0*AB)
The Z value should be
  Z = (b0*AD + a0*DB - d0*AB)*invDet
I do not have a copy of Bowyer and Woodwark, the reference from which the equations are supposedly taken, so I do not know if these errors occur in that reference or were introduced in transcribing.

14 August 2004, page 92. Item "iv. Normal:" should be
  v_perp = v - Dot(u,v)*u/Dot(u,u)
The book has u as the first vector in the subtraction.

13 July 2004, page 367. In the case t0 > 0, the distance from a point Q to a ray P+t*D should be |Q-(P+t0*D)|. The book has |Q-(Q+t0*D)|, a typographical error.

14 June 2004, page 290. The last line of code is missing a left parenthesis. It should be
    v0.y = ((r2 * u.y) + rad) / denom;

14 June 2004, page 291. The first line of pseudocode with solution[0] has v.x when it should be v0.x:
    solution[0].setDir(-v0.y, v0.x);

19 May 2004, page 832. In the elimination step, the line of code
    a_{k,i} = a_{k,i} + m * a_{k,j}
should be
    a_{k,i} = a_{k,i} + m * a_{j,i}

18 May 2004, page 645. The upper bounds on the displayed inequalities are wrong. They should be
  P_{min,x} <= Q_x <= P_{max,x}
  P_{min,y} <= Q_y <= P_{max,y}
  P_{min,z} <= Q_z <= P_{max,z}
The pseudocode on that same page is inconsistent about its reference to the sphere center. I changed "sphere.c" to "sphere.center" in two places.

12 May 2004, page 366. This entry had been listed in the corrections organized by page number, but it was not duplicated in the section organized by date of change.

The pseudocode has a line
  t = Dot(l.direction, VectorSubtract(q, l.direction));
This should be
  t = Dot(l.direction, VectorSubtract(q, l.origin));

09 April 2004, page 187. The displayed equation for the recursion formulas has "B_{i,0}(t)" but should be "B_{i,1}(t)". The displayed equation after that one is the recursion, but should be for "2 <= j <= n" (the text has 1, not 2).

08 March 2004, page 328. The left image of Figure 9.2 has "r = d" and the right image has "r = -d". To be consistent with the paragraph on this page "More significantly...", these should be reversed. The origin in the left image has a negative signed distance from the plane (it is "behind" the plane). The origin in the right image has a positive signed distance from the plane (it is "in front of" the plane).

08 March 2004, page 490. The discussion in this section is about linear components (lines, rays, or line segments) and polygons. The text sometimes talks about lines, sometimes about rays. The pseudocode has a test "if (t < 0 ) { return false; }" which indicates it is about ray-polygon intersection, but is incorrectly named LinePolygonIntersection. Also, there is a line
  p = line.origin + t * r.d;
which should be
  p = line.origin + t * line.direction;
another indication that the author had a ray in mind. Near the end of the pseudocode is
  intersection = line.origin + t * line.direction;
However, the intersection p was calculated earlier, so just use
  intersection = p;

08 March 2004, page 493. The pseudocode has the declaration
  Point3D intersection
This is not needed.

08 March 2004, page 510. The solution for the quadratic formula is missing a minus sign on the "b" term. It should be
  t = (-b +/- sqrt(b^2-4*a*c))/(2*a)
The pseudocode has
  a = (tline.direction.x + tline.direction.y)^2;
This should be
  a = tline.direction.x^2 + tline.direction.y^2;

24 February 2004, page 635. The return statement in the pseudocode is incorrect. It appears to be designed for culling an AABB against the plane rather than indicating if the AABB is intersected by the plane. The return statement also indicates that the plane and AABB are considered not to be intersecting if the box is just touching the plane at a single vertex or single edge or single face. If the plane intersects the box at any point, I suggest replacing the return statement by
  return Dot(plane.normal,dMin)+plane.d <= 0 && Dot(plane.normal,dMax)+plane.d >= 0;
I much prefer the "test intersection" that I use for plane-OBB, which specializes easily when you have an AABB. No branching, just a few arithmetic calculations. Let N = (nx,ny,nz) be the plane normal and d be the plane constant. The plane is Dot(N,X)+d = 0. Let the AABB minima be x0, y0, and z0. Let the maxima be x1, y1, and z1. The test for intersection is
  |nx*(x1+x0)+ny*(y1+y0)+nz*(z1+z0)+2*d| <= |nx*(x1-x0)|+|ny*(y1-y0)|+|nz*(z1-z0)|

30 January 2004, page 59. The third displayed equation should be
  D^2 = (b + mx1 - y1)^2 + (b + mx2 - y2)^2 + (b + mx3 - y3)^2
The book has yn in the last term of the expression.

16 January 2004, page 769. The sentence that starts with "If < V_{i-1}, V_{i+2} > is a diagonal" should have V_{i+1}, not V_{i+2}.

16 January 2004, page 758-759. The input N is number of points, but the first mesh insert call returns a new triangle also named N. These are, of course, different things. For clarity, use Nj as the return value of the two "mesh.Insert" calls, one on page 758 and one on page 759. Use Nj as the argument to the two "stack.Push" calls, one on page 758 and one on page 759.

In the bottom block of pseudocode on page 759, the figure reference should be 13.46, not 13.45. After that line the pseudocode should be
  compute i0, i1, i2, i3 wwith T.v(i1) = A.v(i1) and T.v(i2) = A.v(i2);
In the "if" clause, two lines must be added to the end of the clause
  mesh.Remove(T);
  mesh.Remove(A);
After the two new triangles N0 and N1 are inserted into the mesh, the old triangles T and A must be removed.

14 January 2004, page 322. In the first paragraph of Section 8.16, the equation of the circle is listed as "(x-x_c)+(y-y_c)-r=0". It should be "(x-x_c)^2+(y-y_c)^2-r^2=0".

14 January 2004, page 385. The equation of the plane is listed "ax+by+cz+d=". It should be "ax+by+cz+d=0". (The zero is missing.)

14 January 2004, page 392. The last paragraph has an inequality "n_i * (P - V_{3-i} < 0". It should be "n_i * (P - V_{3-i}) < 0". (The right parenthesis is missing.)

14 January 2004, page 403. Equation (10.9) has "gradient q(X)" on the left-hand side. This should be "gradient q(x)". (The function argument is upper case but needs to be lower case.)

14 January 2004, page 409. Item 4 for the stopping criteria of Newton's method has
  ||(u_{i+1} - u_i S_u(u_i,v_i)) + (v_{i+1} - v_i S_u(u_i,v_i))|| <= epsilon_2
This should be
  ||(u_{i+1} - u_i) S_u(u_i,v_i) + (v_{i+1} - v_i) S_v(u_i,v_i)|| <= epsilon_2
(Two misplaced parentheses. The subscript on the second partial derivative was u, but needs to be v.)

14 January 2004, page 831. The "09 January 2003" entry mentions that the first displayed equation "(1) x1-3*x2+2*x3=6" should be "(1) x1-3*x2-2x3=6". The same correction must be made in the next occurrence of this equation near the bottom of the page (after "which gives us a new system").

14 January 2004, page 832. The pseudocode at the top of the page has a loop "(for k = 1 to n)". This should be "for (k = j+1 to n)". In that loop after the multiplier "m" is computed, you could set a_{k,j} = 0 to make it clear you are eliminating the column entries. If an implementation allows for the row-reduced matrix to be used in further calculations, the assignment of zero should occur in the code. But if the only goal is to solve the system, an implementation need not assign zero since a_{k,j} is never used again in the code.

The first line of pseudocode for back substitution has "x_n = c_n/a_{n,h}". This should be "x_n = c_n/a_{n,n}".

08 January 2004, page 305. The fourth sentence of the first paragraph in Section 8.9 has "r0 + 2*r" and "r0 - 2*r". These should be "r0 + 2*r1" and "r0 - 2*r1".

08 January 2004, page 319. Figure 8.32 has a vector d1 for the top line. This should be d0 since the lines are parallel and are using the same direction vector. The first and second displayed equations are identical. I think the first one was intended to be "L1(t) = P1 + t*d0". The third displayed equation should be "L1(t) = P0 + d*hat(d0)^perp + t*hat(d0)". Following this displayed equation is the statement "if L0 is normalized". The term "normalized line" is not a standard one. The intent is that d0 is a normalized vector, meaning it has length 1.

25 November 2003, page 377. The displayed equation for "e = -e1*(B-P)" should not have the leading minus sign. It should be "e = e1*(B-P)".

25 November 2003, page 380. The pseudocode at the bottom of the page is missing some conditions that lead to processing regions 0 and 5. The first part of the pseudocode should be
  if (s+t <= det)
  {
      if (s < 0)
      {
          if (t < 0)
          {
              region 4
          }
          else
          {
              region 3
          }
      }
      else if (t < 0)
      {
          region 5
      }
      else
      {
          region 0
      }
  }
The document from which that was taken Distance Between Point and Triangle in 3D has the correct information. The algorithm is implemented and available at the source code page at this website (see the Distance web page, files WmlDistVec3Tri3.{h,cpp}).

17 November 2003, page 838. The last set of equations with the d[i][j] should have alpha[2][2] replaced by alpha[2][0] and beta[2][2] replaced by beta[2][0].

All entries below this point are corrected in the second printing of the book.

21 October 2003, page 504. Section 11.3.2 starts out with a construction to find the points of intersection of a line and a sphere. The line is parameterized by P+t*D where D is a unit-length vector and t is a real number. The sphere has center C and radius R. A point X on the sphere satisfies |X-C|^2 = R^2. The intersection calculation is reduced to computing the roots of a quadratic polynomial q(t). The number of intersections is the number of real-valued roots. This is a find-intersection query.

For ray-sphere intersection, it is sufficient to compute the roots of q(t) and keep only those that are nonnegative. The argument is simply that the ray is P+t*D where t >= 0. The top of page 504 describes a method to inexpensively determine if there is an intersection. If the origin P of the ray is inside the sphere, then the ray must intersect the sphere. The point is inside when |P-C|^2 <= R^2. This is a test-intersection query. [If you actually need the point(s) of intersection, you still have to compute the roots of q(t). If P is outside the sphere, the ray might or might not intersect the sphere.]

The next paragraph starts with "In the case of a line segment, the same approach could be used, but instead checking that the root or roots are bounded between 0 and 1, inclusive". First, the statement is in error about the bound 1. The line segment is P+t*D for 0 <= t <= tmax for some positive value tmax that is not necessarily 1 (it is 1 only if the line segment has length 1 since |D| = 1). Ignoring this inaccuracy, in my opinion the implication of the statement is that a test-intersection query will be formulated. The remaining argument is based on a brief article at a link that is now stale (http://astronomy.swin.edu.au/~pbourke/geometry/sphereline/). This article, as well as its restatement in the book, strays from the argument, although the web page has a clear thesis statement: If the line containing the line segment does not intersect the sphere, then the line segment does not intersect the sphere .

Indeed, this is a test-intersection query that attempts to quickly determine there is no intersection. The ray-sphere test-intersection query is an attempt at showing there is an intersection. The essence of the argument given at the web page and in the book is:
  • Compute the point Q on the line that is closest to C.
  • If Q is on the line segment and Q is inside the sphere, then attempt to compute the intersection points using the quadratic formula
The claim is that this can lead to a more efficient intersection calculation. I disagree. The fact that Q is inside the sphere does not guarantee that the line segment intersects the sphere (unless you treat the sphere as a solid). You expend some cycles computing Q and determining it is inside the sphere, only to have to compute the quadratic formula, something you could have done from the start. Moreover, the argument does not tell you what to do if Q is not on the line segment.

Maybe the intent of the article/book was the following:
  • If one segment end point is inside the sphere and the other end point is outside the sphere, the line segment intersects the sphere. Compute the intersection (if desired) by using the quadratic formula.
  • If both segment end points are inside the sphere, then the line segment does not intersect the sphere.
  • If both segment end points are outside the sphere, the line segment intersects the sphere if and only if the distance L from C to the line segment satisfies L <= R. Compute the intersection (if desired) by using the quadratic formula.
At any rate, if you want to adhere to the thesis statement of the article, you compute the squared distance L^2 from the line to C. If L^2 > R^2, there is no intersection between the line segment and the sphere. The construction of L^2 uses only additions, multiplications, and a division (no sqrt calls). You can avoid the division. If L = a/b, then the test for no intersection is a^2 > (b*R)^2.

15 October 2003, page 834. In the paragraph "The row reduction..." is "The condition for existence of solutions is...". The equations should be corrected to
  0 = a0*b1 - a1*b0
    = (c0+c1*y)*b1 - a1*(d0+d1*y)
    = (b1*c0-a1*d0) + (b1*c1-a1*d1)*y
Later is "as long as b1*c1-b0*d1" which should be "as long as b1*c1-a1*d1". Finally, the last equation of the last displayed equation array should be
  = (c0*f0-e0*d0)+((c0*f1-e0*d1)+(c1*f0-e1*d0))*y + (c1*f1-e1*d1)*y^2

14 October 2003, page 125. The first sentence is: "This makes it completely obvious that the perp dot product reflects not only the angle between two vectors but the direction (that is, the sign) of the angle between them". I think what was intended by this is the following. Let u = (u1,u2) and v = (v1,v2). Then u-perp = (-u2,u1) and Dot(u-perp,v) = u1*v2-u2*v1. Write u and v as vectors in 3D, say u' = (u1,u2,0) and v' = (v1,v2,0). Then Cross(u',v') = (0,0,u1*v2-u2*v1) = (0,0,Dot(u-perp,v)). If you reverse the roles of u and v, then Dot(v-perp,u) = -Dot(u-perp,v) and Cross(v',u') = (0,0,Dot(v-perp,u)). The "sign of the angle" is a reference to measuring a counterclockwise (positive) or clockwise (negative) angle between two vectors when you cross the two using the right-hand rule.

14 October 2003, page 135. The revised proof mentioned in the July 25, 2003 entry is better served with a couple of additional steps:
  T(w) = T(Q-P)          by definition of w
       = T(Q) - T(P)     since T is an affine transformation
       = (Q+u) - (P+u)   applying the definition of T
       = Q - P           applying identities on page 80 and 81 to cancel the u terms
       = (P+w) - P       applying the definition of Q on page 134
       = (P-P)+w         by identity iii on page 81
       = zero_vector+w   by identity i on page 80
       = w               zero_vector is the additive identity for vector addition

14 October 2003, page 135. The clause "that is, the first three rows of the matrix..." should say "that is, the first n rows of the matrix...".

14 October 2003, page 138. The sentence before the displayed equation for T(v2) ends with "... transformed coordinates of basis vector v1:". This should end with "v2", not "v1".

14 October 2003, page 141. The last displayed equation has "T(v) = Dot(v,u-hat) + ...". This should be "T(v) = Dot(v,u-hat) u-hat + ...". The derivation of the formula uses identity (iii) on page 91, something that is not mentioned.

14 October 2003, page 152. The displayed equation that occurs after equation (4.32) ends with "u hat". It should end instead with "d hat perp".

14 October 2003, page 162. Equation (4.38) starts with "P = ". It should start with "T(P) = ".

14 October 2003, page 705. The second row of Table 13.1 has the four "Update" entries "i 0 m p". These entries should be "i 0 p m".

20 September 2003, page 24. The second equation should be "2*x1 + 4 = 12 + 17*x2 - 5*x2" in order that the second equation in the first displayed set of equations on page 25 be correct.

20 September 2003, page 51. The third displayed equation has "a(1,1)*x1 + x2*(...)". This should be "a(1,1)*x1 + a(1,2)*(...)".

20 September 2003, page 55. The second sentence of the second paragraph has "... be a set of basis vectors for a vector space vec(v)" where vec(v) denotes the letter v with an arrow over it. This should not be vec(v), but calligraphy(V) which is the symbol in the first paragraph "Given a Eucliden space calligraphy(V)...".

20 September 2003, pages 86, 90, 92. The notation for length of a vector v has been ||v|| (double bars). To be consistent with notation, the following changes should be made.

On page 86, Figure 3.27 has single bars on vector v and should be double bars.

On page 90, Equations (3.8), (3.9), and (3.10) have single bars and should be double bars.

On page 92, Figure 3.31 has single bars on vector v and should be double bars.

20 September 2003, page 94. The construction at the bottom of the page has some inconsistent notation. A couple of the "reasons" on the right-hand side need to be revised.
  • The "x" operator in equations 1 and 2 means "scalar multiplication".
  • The "dot" operator in equations 3 and 4 also means "scalar multiplication".
  • The "x" operator in equations 4, 5, 6 means "cross product".
  • The "dot" operator in equations 5, 6 means "vector dot product".
  • The "reason" for equation 3 should be "by definition of w_perp".
  • The "reason" for equation 5 is meaningless. Since Cross(u,v) and w_perp are parallel vectors that are in the same direction, the angle A between them is 0 radians. The dot product of the two vectors is
        Dot(Cross(u,v),w_perp)
            = Length(Cross(u,v))*Length(w_perp)*cos(A)
            = Length(Cross(u,v))*Length(w_perp)
    
    since cos(A) = 1.

15 September 2003, page 98. Figure 3.35 has a vector 3*v that should be 3*v_1 (the subscript 1 is missing).

15 September 2003, pages 128-129. On page 128, second paragraph from the bottom, the second affine frame is listed as F_G and has origin O_G. These should be F_B and O_B. Figure 4.6 has the origins swapped. What is typeset as O_A should be O_B, and vice versa.

07 September 2003, page 95. The equations at the bottom of page 94 use "psi" as the angle between vectors u and v. Figure 3.34 uses "phi". I changed the figure to use "psi".

04 July 2003, pages 303-304. The equations for the parallel lines at the bottom of the page are missing an "r" that should multiply the square roots.
  a0*x + b0*y + c0 +/- r*sqrt(a0*a0+b0*b0) = 0
  a1*x + b1*y + c1 +/- r*sqrt(a1*a1+b1*b1) = 0
The equations for x and y on page 304 are correct, but the pseudocode has sign errors. It should be
  center[0].x = -(L1.b*(L0.c+discrm0)-L0.b*(L1.c+discrm1))*invDenom;
  center[0].y = +(L1.a*(L0.c+discrm0)-L0.a*(L1.c+discrm1))*invDenom;
  center[1].x = -(L1.b*(L0.c+discrm0)-L0.b*(L1.c-discrm1))*invDenom;
  center[1].y = +(L1.a*(L0.c+discrm0)-L0.a*(L1.c-discrm1))*invDenom;
  center[2].x = -(L1.b*(L0.c-discrm0)-L0.b*(L1.c+discrm1))*invDenom;
  center[2].y = +(L1.a*(L0.c-discrm0)-L0.a*(L1.c+discrm1))*invDenom;
  center[3].x = -(L1.b*(L0.c-discrm0)-L0.b*(L1.c-discrm1))*invDenom;
  center[3].y = +(L1.a*(L0.c-discrm0)-L0.a*(L1.c-discrm1))*invDenom;

02 July 2003, page 304. The pseudocode at the bottom of the page has two occurences of "r" in the body. These should be "radius" to match the input parameter.

25 June 2003, pages 134-135. The second block of displayed equations on page 134 has vector w in two places. These should be vector u:
  T(origin) = origin + u
            = u + origin
The "proof" at the top of page 135 has problems. A formal proof is
  T(w) = T(Q-P)          by definition of w
       = T(Q) - T(P)     since T is an affine transformation
       = (P+w) - P       applying the definition of T
       = (P-P)+w         by identity iii on page 81
       = zero_vector+w   by identity i on page 80
       = w               zero_vector is the additive identity for vector addition

21 June 2003, page 21. The first displayed equation is a dot product whose second term on the right-hand side is a1*b2. This termshould be a2*b2.

20 June 2003, page 698. The pseudocode for PointInSubpolygon has a call to GetMiddleIndex that is missing its last parameter. The call should be
  mid = GetMiddleIndex(i0, i1, C.N);

29 May 2003, pages 535-536. On page 535, the pseudocode for intersection of triangle and plane has a parameter "p2" in the function signature; it should not be there. On page 536, the pseudocode has some test for points either all above the plane or all below the plane. These tests are incorrect. In the style of the remainder of the pseudocode, the tests are simply stated as
  if ( d1 > 0.0 && d2 > 0.0 && d3 > 0.0 )
  {
      // all points above plane
      return false;
  }
  if ( d1 < 0.0 && d2 < 0.0 && d3 < 0.0 )
  {
      // all points below plane
      return false;
  }

28 May 2003, pages 194-195. The displayed equations involving current minimum distance u have the inequalities in the wrong direction. On page 194 the statement
  |x(i) - a| <= u  and  |x(i+1) - a| <= u
should be
  |x(i) - a| >= u  and  |x(i+1) - a| >= u
A similar change for the y-equations on page 194:
  |y(i) - b| <= u  and  |y(i+1) - b| <= u
should be
  |y(i) - b| >= u  and  |y(i+1) - b| >= u
Page 195 has the squared versions that should be
  |x(i) - a|^2 >= u^2  and  |x(i+1) - a|^2 >= u^2
and
  |y(i) - b|^2 >= u^2  and  |y(i+1) - b|^2 >= u^2

26 May 2003, page 106. The end of the first paragraph of section 3.5.1 says "Of course, since the Q_i are representable in terms of A, we could rewrite R in terms of B". The last part should by "in terms of A".

25 May 2003, page 95. The first displayed formula for Vol(u,v,w), case -Dot(Cross(u,v),w) should have w_parallel parallel to -Cross(u,v). However, just as in the May 24th entry below, the use of "parallel" is incorrect terminology. The first case should state "when w_parallel is in the same direction as Cross(u,v)" and the second case should state "when w_parallel is in the opposite direction of Cross(u,v)".

24 May 2003, pages 90-91. The proof of the distributivity of the dot product over vector addition mentions two cases involving a vector u_parallel. I suspect what is intended is that you look at the projections of v, w, and v+w onto the vector u. If the angle between v+w and u is acute, then either the angle between v and u is acute or the angle between w and u is acute (possibly both angles are acute). If angle(v,u) > 0 and angle(w,u) > 0, then
  length(v_perp+w_perp) = length(v_perp) + length(w_perp)
This is case (a) in the book. If angle(v,u) > 0 and angle (w,u) < 0, then
  length(v_perp+w_perp) = length(v_perp) - length(w_perp)
This is case (b) in the book. Of course you must also handle the case when angle(v+w,u) < 0. Two similar cases occur. I also suspect that u_parallel is meant to be the component of u in the u-direction, which is just u itself. The statment u_parallel parallel to v_parallel must mean that angle(u_parallel,v_parallel) = 0. The statment u_parallel not parallel to v_parallel must mean angle (u_parallel,v_parallel) = 180 degrees. The use of the word "parallel" in this context is not good since regardless of the angle you have parallel vectors.

19 May 2003, page 505. The displayed equation for the quadratic coefficient c is missing a subtraction of the squared radius. It should be
  c = k*(Px-Cx)^2 + ell*(Py-Cy)^2 + m*(Pz-Cz)^2 - r^2
The pseudocode at the bottom of the page is also missing the r^2, and the scaling factors need to be members of ellipsoid e. It should be
  c = e.k*(ell.p.x - e.center.x)^2 + e.ell*(ell.p.y - e.center.y)^2 +
      e.m*(ell.p.z - e.center.z)^2 - e.r^2;
The pseudocode uses 'l' (lower-case ell) for variable names, making it difficult to distinguish between '1' (one) and 'l' (ell). The standard axis-aligned representation for an ellipsoid does not separate out the radius term. The standard form is
  ((x-Cx)/Lx)^2 + ((y-Cy)/Ly)^2 + ((z-Cz)/Lz)^2 = 1
where the center is (Cx,Cy,Cz) and Lx, Ly, and Lz are the lengths measured along the axes from the center to the extreme points of the ellipsoid in those directions. I prefer more verbose names for the member names. The line has an origin and a direction. The ellipsoid has a center and length values. Another pass at the pseudocode is given below.
  int LineEllipsoidIntersect (Ellipsoid E, Line3 L, float t[2])
  {
      Vector3 diff = L.origin - E.center;
      float invLx = 1/E.Lx;
      float invLy = 1/E.Ly;
      float invLz = 1/E.Lz;
      float tmp1x = L.direction.x * invLx;
      float tmp1y = L.direction.y * invLy;
      float tmp1z = L.direction.z * invLz;
      float tmp2x = diff.x * invLx;
      float tmp2y = diff.y * invLy;
      float tmp2z = diff.z * invLz;
      float a = tmp1x * tmp1x + tmp1y * tmp1y + tmp1z * tmp1z;
      float b = 2 * (tmp1x * tmp2x + tmp1y * tmp2y + tmp1z * tmp2z);
      float c = tmp2x * tmp2x + tmp2y * tmp2y + tmp2z * tmp2z - 1;
      float d = b * b - 4 * a * c;
      if ( d > 0 )
      {
          float invTwoA = 1/(2*a);
          float rootD = sqrt(d);
          t[0] = (-b + rootD)*invTwoA;
          t[1] = (-b - rootD)*invTwoA;
          return 2;
      }
      else if ( discrm == 0 )
      {
          t[0] = -b/(2*a);
          return 1;
      }
      return 0;  // no solutions
  }
The return value of the function is the number of solutions. If t[i] is valid (0 <= i < numberOfSolutions), then the corresponding point of intersection is: L.origin + t[i]*L.direction.

13 May 2003, page 487. The top of the page makes mention that "u+v <= 1" is necessary for the intersection, but the pseudocode does not make this check. The block of code starting with the comparison of v to 0 and 1 and ending with the assignment of t should be
  if ( v < 0.0 || v > 1.0 )
      return false;
  if ( u + v > 1.0 )
      return false;
  t = tmp * Dot(e2,q);

08 Apr 2003, page 410, 411, and 412. On page 410, the first displayed equation after (10.13) is
  v = Q0 - Q1
    = P0 + sc*d0 - P1 + tc*d1
should have parentheses
  v = Q0 - Q1
    = P0 + sc*d0 - (P1 + tc*d1)
However, the pair of displayed equations following these are correct.

The third line from the bottom defines "b = -d0*d1", but should be "b = d0*d1" to be consistent with the pseudocode.

The next to last line is missing a parenthesis and should be "f = (P0-P1)*(P0-P1)".

On page 411, with the corrected "b = d0*d1", the formula at the top of the page for tc is wrong. It should be
  tc = (a*e-b*d)/(a*c-b*b)
The formula is wrong even if you used the incorrect "b = -d0*d1". This change propagates to the next displayed equation:
  |L0(sc)-L1(tc)| = |(P0-P1)+ ((b*e-c*d)d0-(a*e-b*d)d1)/(a*c-b*b)|
The block of pseudocode for "det < epsilon" has "return d*s+f" which is wrong. It should be
"return -e*t+f". On page 412, the top line has the correct formula for t despite the errors in the text. However, the return value is wrong and should be
  return s*(a*s-b*t+2*d)+t*(-b*s+c*t-2*e)+f;

25 Mar 2003, page 366. The pseudocode has a line
  t = Dot(l.direction, VectorSubtract(q, l.direction));
This should be
  t = Dot(l.direction, VectorSubtract(q, l.origin));

15 Feb 2003, page 487. The pseudocode has a line near the bottom
  v = tmp * Dot(d, q);
This should be
  v = tmp * Dot(line.direction, q);

02 Feb 2003, page 628-629. The pseudocode on page 628, the case where "Ray parallel to planes" has a test "P.x < box.xMax". That should be "P.x > box.xMax". The pseudocode on page 629 for swapping t0 and t1 is incorrect. It should be "tmp = t1; t1 = t0; t0 = tmp;"

02 Feb 2003, page 502. The quadratic equation in the second line should say "at^2+bt+c=0" to be consistent with the footnote and the pseudocode. Equation (11.11) should be
    t = -d*(P-C) +- sqrt((d*(P-C))^2-((P-C)*(P-C)-r^2))
and is based on d being a unit-length vector (a = 1 in the quadratic equation). Note, though, that the pseudocode does not assume d is unit length. The website does have an implementation of this algorithm, the Intersection 3D page. (The expression in the line after equation (11.11) has an extraneous right parenthesis.)

09 Jan 2003, page 831. The fourth line from the top of the page is "(1) x1 - 3*x2 + 2*x3 = 6". It should be "(1) x1 - 3*x2 - 2*x3 = 6". The last line of the page lists the solution as (1,-3,-2). It should be (1,-3,2).

17 Dec 2002, page 67. The distributive rule (iii) should be (a+b)u = au + bu (the text has bv).

12 Dec 2002, pages 641 and 642. The last line of page 641 has a lower-case subscript 'a' on the unit vector u. That should be an upper-case 'A'. The last line of page 642 also has a lower-case subscript of 'a', but it should be an upper-case 'B'.

11 Oct 2002, page 77. The last sentence of the first paragraph states "Many values in the range..." but should be "Many values in the domain..."

Book Corrections Organized by Page Number

Section 2.2, page 15. The "Matrix" row of Table 2.1 has pseudocode "Matrix3x3 m, m;" that should be "Matrix3x3 m, n;".

Section 2.3.4, page 21. The first displayed equation is a dot product whose second term on the right-hand side is a1*b2. This termshould be a2*b2.

Section 2.4.1, page 24. The second equation should be "2*x1 + 4 = 12 + 17*x2 - 5*x2" in order that the second equation in the first displayed set of equations on page 25 be correct.

Section 2.6.3, page 43. Bullet 'v' in Section 2.6.3 is a repeat of bullet 'ii' and should be removed.

Section 2.6.5, page 44. The example combination of vectors in a displayed equation is
  2v1 + 3v3 +-1v3
but should be
  2v1 + 3v2 +-1v3
That is, the subscript on the middle term needed to be fixed.

Section 2.7.1, page 46. The sentence: "We can map real numbers to real numbers: ..." should say "We can map nonnegative real numbers fo nonnegative real numbers: ..."

Section 2.7.4, page 51. The third displayed equation has "a(1,1)*x1 + x2*(...)". This should be "a(1,1)*x1 + a(1,2)*(...)".

Section 2.8, page 53. The last sentence should say "... into Equations 2.4 and 2.5, ...".

Section 2.9.2, page 55. The second sentence of the second paragraph has "... be a set of basis vectors for a vector space vec(v)" where vec(v) denotes the letter v with an arrow over it. This should not be vec(v), but calligraphy(V) which is the symbol in the first paragraph "Given a Eucliden space calligraphy(V)...".

Section 2.10, page 59. The third displayed equation should be
  D^2 = (b + mx1 - y1)^2 + (b + mx2 - y2)^2 + (b + mx3 - y3)^2
The book has yn in the last term of the expression.

Section 3.1.5, page 67. The distributive rule (iii) should be (a+b)u = au + bu (the text has bv).

Section 3.2.6, page 77. The last sentence of the first paragraph states "Many values in the range..." but should be "Many values in the domain..."

Section 3.2.6, page 78. Figure 3.19 is confusing and has one typographical error. It is an attempt to describe the linear transformation T(x) = 2*x. In terms of basis vectors u and v: T(u) = 2*u and T(v) = 2*v. The vector s in the figure should be more clearly written:
  s = T(2*v)
    = 2*T(v)   // because T is linear
    = 2*(2*v)  // by definition of T
Call the top vector of the figure r, in which case
  r = T(3*u)
    = 3*T(u)   // because T is linear
    = 3*(2*u)  // by definition of T
Thus, the current figure should have "3T(u) = 2(3u)".

Section 3.2.6, page 79. Figure 3.20 is confusing and has two typographical errors. It is an attempt to describe the linear transformation that is defined on the basis vectors u and v: T(u) = 2*u and T(v) = 1.5*v. Given a vector x = a*u+b*v, the transformation is T(x) = T(a*u+b*v) = a*T(u)+b*T(v) = a*(2*u)+b*(1.5*v). The vector s in the figure is incorrect and should be
  s = T(2*v)
    = 2*T(v)   // because T is linear
    = 2*(1.5*v)  // by definition of T
Call the top vector of the figure r, in which case
  r = T(3*u)
    = 3*T(u)   // because T is linear
    = 3*(2*u)  // by definition of T
Thus, the current figure should have "3T(u) = 2(3u)".

Section 3.3, page 83. The sentence starting with "We've just done two forbidden operations on points--scaling by a vector and directly adding them." This should say "We've just done two forbidden operations on points--multiplying points by scalars and directly adding them."

Section 3.3, page 84. The first displayed equation has a term "alpha3*(P3-P2)". This term should be "alpha3*(P3-P1)".

Section 3.3.1, pages 86, 90, 92. The notation for length of a vector v has been ||v|| (double bars). To be consistent with notation, the following changes should be made.

On page 86, Figure 3.27 has single bars on vector v and should be double bars.

On page 90, Equations (3.8), (3.9), and (3.10) have single bars and should be double bars.

On page 92, Figure 3.31 has single bars on vector v and should be double bars.

Section 3.3.1, pages 90-91. The proof of the distributivity of the dot product over vector addition mentions two cases involving a vector u_parallel. I suspect what is intended is that you look at the projections of v, w, and v+w onto the vector u. If the angle between v+w and u is acute, then either the angle between v and u is acute or the angle between w and u is acute (possibly both angles are acute). If angle(v,u) > 0 and angle(w,u) > 0, then
  length(v_perp+w_perp) = length(v_perp) + length(w_perp)
This is case (a) in the book. If angle(v,u) > 0 and angle (w,u) < 0, then
  length(v_perp+w_perp) = length(v_perp) - length(w_perp)
This is case (b) in the book. Of course you must also handle the case when angle(v+w,u) < 0. Two similar cases occur. I also suspect that u_parallel is meant to be the component of u in the u-direction, which is just u itself. The statment u_parallel parallel to v_parallel must mean that angle(u_parallel,v_parallel) = 0. The statment u_parallel not parallel to v_parallel must mean angle (u_parallel,v_parallel) = 180 degrees. The use of the word "parallel" in this context is not good since regardless of the angle you have parallel vectors.

Section 3.3.1, page 92. Item "iv. Normal:" should be
  v_perp = v - Dot(u,v)*u/Dot(u,u)
The book has u as the first vector in the subtraction.

Section 3.3.2, page 94. The construction at the bottom of the page has some inconsistent notation. A couple of the "reasons" on the right-hand side need to be revised.
  • The "x" operator in equations 1 and 2 means "scalar multiplication".
  • The "dot" operator in equations 3 and 4 also means "scalar multiplication".
  • The "x" operator in equations 4, 5, 6 means "cross product".
  • The "dot" operator in equations 5, 6 means "vector dot product".
  • The "reason" for equation 3 should be "by definition of w_perp".
  • The "reason" for equation 5 is meaningless. Since Cross(u,v) and w_perp are parallel vectors that are in the same direction, the angle A between them is 0 radians. The dot product of the two vectors is
        Dot(Cross(u,v),w_perp)
            = Length(Cross(u,v))*Length(w_perp)*cos(A)
            = Length(Cross(u,v))*Length(w_perp)
    
    since cos(A) = 1.

Section 3.3.2, page 95. The equations at the bottom of page 94 use "psi" as the angle between vectors u and v. Figure 3.34 uses "phi". I changed the figure to use "psi".

Section 3.3.2, page 95. The first displayed formula for Vol(u,v,w), case -Dot(Cross(u,v),w) should have w_parallel parallel to -Cross(u,v). However, just as in the May 24th entry below, the use of "parallel" is incorrect terminology. The first case should state "when w_parallel is in the same direction as Cross(u,v)" and the second case should state "when w_parallel is in the opposite direction of Cross(u,v)".

Section 3.3.3, page 98. Figure 3.35 has a vector 3*v that should be 3*v_1 (the subscript 1 is missing).

Section 3.5.1, page 106. The end of the first paragraph of section 3.5.1 says "Of course, since the Q_i are representable in terms of A, we could rewrite R in terms of B". The last part should by "in terms of A".

Section 4.4.2, page 120. In the middle of the page, I do not understand why the statement "with, as before, Vector(u) retaining its usual matrix representation" is made. This seems to imply that Vector(u) will be represented as a row matrix as shown at the top of the page. When Cross(Vector(v),Vector(u)) is represented using a skew-symmetric matrix for Vector(v), Vector(u) must be written as a column matrix. Moreover, the skew-symmetric matrix Tilde(v) is the negative of what it should be. Listing the rows from top to bottom, it should be Tilde(v) = ((0,-v3,+v2),(+v3,0,-v1),(-v2,+v1,0)).

Section 4.4.4, page 125. The first sentence is: "This makes it completely obvious that the perp dot product reflects not only the angle between two vectors but the direction (that is, the sign) of the angle between them". I think what was intended by this is the following. Let u = (u1,u2) and v = (v1,v2). Then u-perp = (-u2,u1) and Dot(u-perp,v) = u1*v2-u2*v1. Write u and v as vectors in 3D, say u' = (u1,u2,0) and v' = (v1,v2,0). Then Cross(u',v') = (0,0,u1*v2-u2*v1) = (0,0,Dot(u-perp,v)). If you reverse the roles of u and v, then Dot(v-perp,u) = -Dot(u-perp,v) and Cross(v',u') = (0,0,Dot(v-perp,u)). The "sign of the angle" is a reference to measuring a counterclockwise (positive) or clockwise (negative) angle between two vectors when you cross the two using the right-hand rule.

Section 4.6, pages 128-129. On page 128, second paragraph from the bottom, the second affine frame is listed as F_G and has origin O_G. These should be F_B and O_B. Figure 4.6 has the origins swapped. What is typeset as O_A should be O_B, and vice versa.

Section 4.7.2, pages 134-135. The second block of displayed equations on page 134 has vector w in two places. These should be vector u:
  T(origin) = origin + u
            = u + origin
The "proof" at the top of page 135 has problems. A formal proof is
  T(w) = T(Q-P)          by definition of w
       = T(Q) - T(P)     since T is an affine transformation
       = (P+w) - P       applying the definition of T
       = (P-P)+w         by identity iii on page 81
       = zero_vector+w   by identity i on page 80
       = w               zero_vector is the additive identity for vector addition

Section 4.7.2, page 135. The revised proof mentioned in the previous entry is better served with a couple of additional steps:
  T(w) = T(Q-P)          by definition of w
       = T(Q) - T(P)     since T is an affine transformation
       = (Q+u) - (P+u)   applying the definition of T
       = Q - P           applying identities on page 80 and 81 to cancel the u terms
       = (P+w) - P       applying the definition of Q on page 134
       = (P-P)+w         by identity iii on page 81
       = zero_vector+w   by identity i on page 80
       = w               zero_vector is the additive identity for vector addition

Section 4.7.2, page 135. The clause "that is, the first three rows of the matrix..." should say "that is, the first n rows of the matrix...".

Section 4.7.3, page 138. The sentence before the displayed equation for T(v2) ends with "... transformed coordinates of basis vector v1:". This should end with "v2", not "v1".

Section 4.7.3, page 141. The last displayed equation has "T(v) = Dot(v,u-hat) + ...". This should be "T(v) = Dot(v,u-hat) u-hat + ...". The derivation of the formula uses identity (iii) on page 91, something that is not mentioned.

Section 4.7.5, page 152. The displayed equation that occurs after equation (4.32) ends with "u hat". It should end instead with "d hat perp".

Section 4.7.6, pages 157-158. For simplicity of notation here, define N = hat(n), N^T is the transpose of N, and V = hat(v). The last row, first column entry of the general shearing matrix should be -(tan(theta))*Dot(N,Q)*V. The first line on page 158 makes an appeal to applying the point transformation T to a vector N, which makes no sense mathematically. Even if it did, it is not clear why the author claims T causes N to change at all. Let P = P1 + h*N, where h = Dot(N,P-Q). The shearing is P" = P + t*V, where from trigonometry, t = h*tan(theta), so P" = P + (tan(theta))*Dot(N,P-Q)*V and P"-Q = [I + (tan(theta))*V*N^T](P-Q). Add Q to both sides and factor to P" = [I + (tan(theta))*V*N^T]*P - (tan(theta))*V*N^T*Q = [I + (tan(theta))*V*N^T]*P - (tan(theta))*Dot(N,Q)*V.

Section 4.8.2, page 162. Equation (4.38) starts with "P = ". It should start with "T(P) = ".

Section 5.6.2, page 187. The displayed equation for the recursion formulas has "B_{i,0}(t)" but should be "B_{i,1}(t)". The displayed equation after that one is the recursion, but should be for "2 <= j <= n" (the text has 1, not 2).

Section 6.2, pages 194-195. The displayed equations involving current minimum distance u have the inequalities in the wrong direction. On page 194 the statement
  |x(i) - a| <= u  and  |x(i+1) - a| <= u
should be
  |x(i) - a| >= u  and  |x(i+1) - a| >= u
A similar change for the y-equations on page 194:
  |y(i) - b| <= u  and  |y(i+1) - b| <= u
should be
  |y(i) - b| >= u  and  |y(i+1) - b| >= u
Page 195 has the squared versions that should be
  |x(i) - a|^2 >= u^2  and  |x(i+1) - a|^2 >= u^2
and
  |y(i) - b|^2 >= u^2  and  |y(i+1) - b|^2 >= u^2

Section 6.3.1, pages 208-210. The algorithm uses barycentric coordinates to determine closest edge, but this does not work for a triangle with an obtuse angle. Thus, the SquaredDistance pseudocode is incorrect. The online source code does not use this algorithm; it uses the previously discussed (correct) algorithm.

Section 6.3.3, page 213. In the next to last paragraph, the equation for a solid cone is given as Dot(hat(a),X-V) >= cos(theta) but should be Dot(hat(a),X-V) >= cos(theta)*Length(X-V).

Section 6.6.4, page 225. The denominator of Equation (6.17) should be dot(perp(d0),d1) rather than dot(perp(d1),d0); this has the effect of negating the right-hand side of the equation.

Section 6.6.6, page 228. The symbols in Equation (6.23) on page 228 are not typeset correctly. In that equation, let us say that the first row of the equation is region R0 and the last row of the equation is R8 (the parameter plane is partitioned into 9 regions). The corrections are
    R2:  0 < tilde(t0) < T0 and bar(t1) >= T1
    R4:  0 < tilde(t1) < T1 and bar(t0) >= T0
    R6:  hat(t0) >= T0 and tilde(t1) <= 0
    R7:  tilde(t0) <= 0 and hat(t1) >= T1
    R8:  tilde(t0) >= T0 and tilde(t1) >= T1

Section 6.6.6, page 229. Equation (6.24) has 5 cases (top to bottom). The final inequality in case 2 should be Dot(d0,delta+T0*d0) ≤ 0. The final inequality in case 4 should be Dot(d0,delta+T0*d0-T1*d1) ≤ 0.

Section 7.7.2, page 272. In the GetExtremeIndex function, the input direction D is listed as type Point but should be type Vector. The call to GetMiddleIndex has only two parameters but should have three parameters: mid = GetMiddleIndex(i0,i1,C.N);

Section 8.4, page 290. The last line of code is missing a left parenthesis. It should be
    v0.y = ((r2 * u.y) + rad) / denom;

Section 8.4, page 291. The first line of pseudocode with solution[0] has v.x when it should be v0.x:
    solution[0].setDir(-v0.y, v0.x);

Section 8.7, page 301-302. The pseudocode for CircleThroughPointTangentToLineGivenRadius has several typographical errors and code errors. On page 301, the "if (cPrime < 0)" block assigns the negative line direction to (a,b,c). The block also needs to negate cPrime. The declaration of tmp2 has a sum of squares, so it can never be less than -epsilon. It should be "tmp2 = r * r - tmp1 * tmp1". There are several occurrences of "r" in the pseudocode. These should be "radius". On page 302, the last two assignments of the center values should be "center[0] = tmpPt + tmp2 * D" and "center[1] = tmpPt - tmp2 * D" where "Point2D D = { b, -a };"

Section 8.8, pages 303-304. The equations for the parallel lines at the bottom of the page are missing an "r" that should multiply the square roots.
  a0*x + b0*y + c0 +/- r*sqrt(a0*a0+b0*b0) = 0
  a1*x + b1*y + c1 +/- r*sqrt(a1*a1+b1*b1) = 0
The equations for x and y on page 304 are correct, but the pseudocode has sign errors. It should be
  center[0].x = -(L1.b*(L0.c+discrm0)-L0.b*(L1.c+discrm1))*invDenom;
  center[0].y = +(L1.a*(L0.c+discrm0)-L0.a*(L1.c+discrm1))*invDenom;
  center[1].x = -(L1.b*(L0.c+discrm0)-L0.b*(L1.c-discrm1))*invDenom;
  center[1].y = +(L1.a*(L0.c+discrm0)-L0.a*(L1.c-discrm1))*invDenom;
  center[2].x = -(L1.b*(L0.c-discrm0)-L0.b*(L1.c+discrm1))*invDenom;
  center[2].y = +(L1.a*(L0.c-discrm0)-L0.a*(L1.c+discrm1))*invDenom;
  center[3].x = -(L1.b*(L0.c-discrm0)-L0.b*(L1.c-discrm1))*invDenom;
  center[3].y = +(L1.a*(L0.c-discrm0)-L0.a*(L1.c-discrm1))*invDenom;

Section 8.8, page 304. The pseudocode at the bottom of the page has two occurences of "r" in the body. These should be "radius" to match the input parameter.

Section 8.9, page 305. The fourth sentence of the first paragraph in Section 8.9 has "r0 + 2*r" and "r0 - 2*r". These should be "r0 + 2*r1" and "r0 - 2*r1".

Section 8.14, page 319. Figure 8.32 has a vector d1 for the top line. This should be d0 since the lines are parallel and are using the same direction vector. The first and second displayed equations are identical. I think the first one was intended to be "L1(t) = P1 + t*d0". The third displayed equation should be "L1(t) = P0 + d*hat(d0)^perp + t*hat(d0)". Following this displayed equation is the statement "if L0 is normalized". The term "normalized line" is not a standard one. The intent is that d0 is a normalized vector, meaning it has length 1.

Section 8.16, page 322. In the first paragraph of Section 8.16, the equation of the circle is listed as "(x-x_c)+(y-y_c)-r=0". It should be "(x-x_c)^2+(y-y_c)^2-r^2=0".

Section 9.2.1, page 328. The left image of Figure 9.2 has "r = d" and the right image has "r = -d". To be consistent with the paragraph on this page "More significantly...", these should be reversed. The origin in the left image has a negative signed distance from the plane (it is "behind" the plane). The origin in the right image has a positive signed distance from the plane (it is "in front of" the plane).

Section 9.5, page 356. The first equation for the torus is incorrect. It should be the following equation. I have factored it in a manner that leads to more efficient numerical computation:
  (x^2 + y^2 + z^2)^2 - 2*(r0^2 + r1^2)*(x^2 + y^2 + z^2) + 4*r0^2*z^2 + (r0^2 - r1^2)^2 = 0
My torus class (in Wm3Torus files) has the correct formula listed.

Section 10.2, page 366. The pseudocode has a line
  t = Dot(l.direction, VectorSubtract(q, l.direction));
This should be
  t = Dot(l.direction, VectorSubtract(q, l.origin));

Section 10.3.2, page 377. The displayed equation for "e = -e1*(B-P)" should not have the leading minus sign. It should be "e = e1*(B-P)".

Section 10.3.2, page 377. The equation at the bottom of the page for bar(t) has the wrong numerator. The equation should be "bar(t) = (b*d-a*e)/(a*c-b*b)". The pseudocode on page 380 does use the correct equation.

Section 10.2.1, page 367. In the case t0 > 0, the distance from a point Q to a ray P+t*D should be |Q-(P+t0*D)|. The book has |Q-(Q+t0*D)|, a typographical error.

Section 10.3, page 380. The pseudocode at the bottom of the page is missing some conditions that lead to processing regions 0 and 5. The first part of the pseudocode should be
  if (s+t <= det)
  {
      if (s < 0)
      {
          if (t < 0)
          {
              region 4
          }
          else
          {
              region 3
          }
      }
      else if (t < 0)
      {
          region 5
      }
      else
      {
          region 0
      }
  }
The document from which that was taken Distance Between Point and Triangle in 3D has the correct information. The algorithm is implemented and available at the source code page at this website (see the Distance web page, files WmlDistVec3Tri3.{h,cpp}).

Section 10.3.3, page 384. The first displayed equation has a Q that should be a P:
    Q' = P + s*e0 + t*e1

Section 10.3.4, page 385. The equation of the plane is listed "ax+by+cz+d=". It should be "ax+by+cz+d=0". (The zero is missing.)

Section 10.3.4, page 387. The top line of pseudocode has "Dot(p,n)-d", but should be "Dot(p,n)+d".

Section 10.4.1, page 392. The last paragraph has an inequality "n_i * (P - V_{3-i} < 0". It should be "n_i * (P - V_{3-i}) < 0". (The right parenthesis is missing.)

Section 10.5.2, page 403. Equation (10.9) has "gradient q(X)" on the left-hand side. This should be "gradient q(x)". (The function argument is upper case but needs to be lower case.)

Section 10.7, page 409. Item 4 for the stopping criteria of Newton's method has
  ||(u_{i+1} - u_i S_u(u_i,v_i)) + (v_{i+1} - v_i S_u(u_i,v_i))|| <= epsilon_2
This should be
  ||(u_{i+1} - u_i) S_u(u_i,v_i) + (v_{i+1} - v_i) S_v(u_i,v_i)|| <= epsilon_2
(Two misplaced parentheses. The subscript on the second partial derivative was u, but needs to be v.)

Section 10.8.1, page 410, 411, and 412. On page 410, the first displayed equation after (10.13) is
  v = Q0 - Q1
    = P0 + sc*d0 - P1 + tc*d1
should have parentheses
  v = Q0 - Q1
    = P0 + sc*d0 - (P1 + tc*d1)
However, the pair of displayed equations following these are correct.

The third line from the bottom defines "b = -d0*d1", but should be "b = d0*d1" to be consistent with the pseudocode.

The next to last line is missing a parenthesis and should be "f = (P0-P1)*(P0-P1)".

On page 411, with the corrected "b = d0*d1", the formula at the top of the page for tc is wrong. It should be
  tc = (a*e-b*d)/(a*c-b*b)
The formula is wrong even if you used the incorrect "b = -d0*d1". This change propagates to the next displayed equation:
  |L0(sc)-L1(tc)| = |(P0-P1)+ ((b*e-c*d)d0-(a*e-b*d)d1)/(a*c-b*b)|
The block of pseudocode for "det < epsilon" has "return d*s+f" which is wrong. It should be
"return -e*t+f". On page 412, the top line has the correct formula for t despite the errors in the text. However, the return value is wrong and should be
  return s*(a*s-b*t+2*d)+t*(-b*s+c*t-2*e)+f;

Section 10.8, pages 417-426. These pages cover distance calculations between linear components. The pseudocode has major problems. Apparently, one of these blocks of code was written and then cut-paste-modified for the other blocks. They all exhibit the same ills. The difference between closest points needs to be of the form
  v = (base0 + s*direction0) - (base1 + t*direction1)
That is, the book parentheses are incorrect. Worse is the handling of the distinction between parallel and nonparallel components. If "det" is ever zero, the pseudocode will do a division by it. Also, the epsilon test is problematic when the segment direction vectors have large length. If D0 and D1 are the line segment lengths, each computed as a difference of the line segment end points, and if A is the angle between the directions, it is the case that
  det = |Cross(D0,D1)|^2 = |D0|^2 * |D1|^2 * (sin(A))^2
If the angle A is very small, but the lengths of the directions are very large, "det" can be quite large, causing you to classify the segments are nonparallel when they are nearly parallel.

My own implementation uses a line segment represention of C+s*D, where C is the midpoint of the segment, D is a unit-length direction, and |s| <= e where e is half the length of the segment. This leads to a test
  det = (sin(A))^2
so that the epsilon test is a relative error measurement, not an absolute error measurement. Section 10.8.3 of the book describes how I think of distance measurements for linear components, and I much prefer this approach based on actually writing code and seeing how it performs on real data.

Section 10.8.2, page 419. The pseudocode for c should be
  c = Dot(ray.direction, ray.direction);

Section 10.9.1, page 437. The displayed matrix equation is incorrect. Setting the gradient to zero leads to the system of equations
    a00*u + a01*v + a02*t + b0 = 0
    a01*u + a11*v + a12*t + b1 = 0
    a02*u + a12*v + a22*t + b2 = 0
The pseudocode on that page essentially solves this system using Cramer's rule (computes cofactors and determinant).

Section 10.8.3, page 428. The paragraph about the case of region 2 has "the gradient at the corner (1,1) cannot point inside the unit square". This should say "outside" rather than "inside".

Section 10.9.2, pages 442-446. Section 10.9.2 for computing the distance between a line and a rectangle in 3D is filled with errors--what an atrocious description. As a replacement, please use Distance from Line to Rectangle in 3D, which describes a robust and efficient algorithm for computing the distance.
  • On page 442, Figure 10.42 is misleading because the line is drawn as if it is neither parallel nor perpendicular to the rectangle, yet the closest point Q is shown to be on the rectangle. This is not possible for a non-parallel/non-perpendicular line.
  • On page 442, the displayed |R(u,v) - L(t)| is missing a power 2..
  • On page 442, in the second displayed equation involving Q(u,v,t) it would have been helpful to use Δ for V-P in order to simplify reading the equations. Of course, things are worse in that the equation is wrong. It should be
    (e0 • e0)u2 + (e1 • e1)v2 + (d • d)t2 + 2(e0 • e1)u v - 2(e0 • d)u t - 2(e1 • d)v t
     + 2(e0 • Δ)u + 2(e1 • Δ)v - 2(d • Δ)t + |Δ|2
    The definitions for the quadratic coefficients on page 443 must be adjusted accordingly.
    a00 = (e0 • e0), a01 = 2(e0 • e1), a11 = (e1 • e1), a02 = -2(e0 • d), a12 = -2(e1 • d),
    a22 = (d • d), b0 = 2(e0 • Δ), b1 = 2(e1 • Δ), b2 = -2(d • Δ), c = |Δ|2 
    The book's a12 is listed twice; the first occurrence is wrong. The book's b2 is also wrong. The book also failed to define c in that list.
  • On page 443, "...then the solution (bar(u),bar(v),bar(t))..." does not say what problem is being solved. It is the solution to the linear system of equations mentioned later. Perhaps "...then the minimizer (bar(u),bar(v),bar(t))..." is better said.
  • On page 443, there is an expression |T(u,v,t)-L(t)| that should be |R(u,v)-L(t)|.
  • On page 443, "Otherwise, the minimum of ∇ Q = (bar(u),bar(v),bar(t))..." is incorrect terminology. The minimum we seek is for Q. For a convex quadratic with unconstrained domain, the global minimum occurs when ∇ Q is (0,0,0), which leads to a minimizer (bar(u),bar(v),bar(t)). In the constrained domain, the minimum can occur on a face of the convex domain, but any gradient is computed in the domain of that face. Better is "Otherwise, the minimizer (bar(u),bar(v),bar(t)) of Q will be on a face...".
  • Speaking of that linear system at the bottom of page 443, it is incorrectly written. The variables are written on the right-hand side and the b-vector on the left. These are swapped, the b-vector needs to be negated and the matrix of coefficients should be written to show the symmetry.
  • On page 444, a01, a02, and a12> each needs an additional factor of 2. On page 445, b0 and b1 each needs an additional factor of, and b2 needs an additional factor of -2.
  • On page 445, The block of pseudocode that flips signs when the determinant is negative does not make sense. Theoretically, the determinant is nonnegative, and it can be computed with an algorithm different from that of the pseudocode to guarantee that it is nonnegative even with floating-point rounding errors.
  • On page 445, the test for "(near) parallelism" is one way to try to avoid floating-point rounding errors when computing the distance. But as the PDF mentioned previously points out, subtractive cancellation can cause gross errors in the result. I recommend using the constrained conjugate gradient algorithm of that PDF.
  • On page 445, there are terms "rectangle[i]" for i in {0,1,2,3}. Previously in the pseudocode, the notation was "rectangle.V[i]".
  • On page 446, the pseudocode first computes u and v but then the conditional tests involves s and t. The s and t are meant to be u and v, respectively.

Section 10.9.4, page 454. The pseudocode after the "Closest point on the box" comment should be
  qPrime.x = box.extents.x;
  qPrime.y = line.origin.y;
  qPrime.z = line.origin.z;

Section 11.1, page 487. The pseudocode has a line near the bottom
  v = tmp * Dot(d, q);
This should be
  v = tmp * Dot(line.direction, q);

Section 11.1, page 487. The top of the page makes mention that "u+v <= 1" is necessary for the intersection, but the pseudocode does not make this check. The block of code starting with the comparison of v to 0 and 1 and ending with the assignment of t should be
  if ( v < 0.0 || v > 1.0 )
      return false;
  if ( u + v > 1.0 )
      return false;
  t = tmp * Dot(e2,q);

Section 11.1.3, page 490. The discussion in this section is about linear components (lines, rays, or line segments) and polygons. The text sometimes talks about lines, sometimes about rays. The pseudocode has a test "if (t < 0 ) { return false; }" which indicates it is about ray-polygon intersection, but is incorrectly named LinePolygonIntersection. Also, there is a line
  p = line.origin + t * r.d;
which should be
  p = line.origin + t * line.direction;
another indication that the author had a ray in mind. Near the end of the pseudocode is
  intersection = line.origin + t * line.direction;
However, the intersection p was calculated earlier, so just use
  intersection = p;

Section 11.1.4, page 493. The pseudocode has the declaration
  Point3D intersection
This is not needed.

Section 11.3, page 502. The quadratic equation in the second line should say "at^2+bt+c=0" to be consistent with the footnote and the pseudocode. Equation (11.11) should be
    t = -d*(P-C) +- sqrt((d*(P-C))^2-((P-C)*(P-C)-r^2))
and is based on d being a unit-length vector (a = 1 in the quadratic equation). Note, though, that the pseudocode does not assume d is unit length. The website does have an implementation of this algorithm, the Intersection page. (The expression in the line after equation (11.11) has an extraneous right parenthesis.)

Section 11.3.2, page 504. Section 11.3.2 starts out with a construction to find the points of intersection of a line and a sphere. The line is parameterized by P+t*D where D is a unit-length vector and t is a real number. The sphere has center C and radius R. A point X on the sphere satisfies |X-C|^2 = R^2. The intersection calculation is reduced to computing the roots of a quadratic polynomial q(t). The number of intersections is the number of real-valued roots. This is a find-intersection query.

For ray-sphere intersection, it is sufficient to compute the roots of q(t) and keep only those that are nonnegative. The argument is simply that the ray is P+t*D where t >= 0. The top of page 504 describes a method to inexpensively determine if there is an intersection. If the origin P of the ray is inside the sphere, then the ray must intersect the sphere. The point is inside when |P-C|^2 <= R^2. This is a test-intersection query. [If you actually need the point(s) of intersection, you still have to compute the roots of q(t). If P is outside the sphere, the ray might or might not intersect the sphere.]

The next paragraph starts with "In the case of a line segment, the same approach could be used, but instead checking that the root or roots are bounded between 0 and 1, inclusive". First, the statement is in error about the bound 1. The line segment is P+t*D for 0 <= t <= tmax for some positive value tmax that is not necessarily 1 (it is 1 only if the line segment has length 1 since |D| = 1). Ignoring this inaccuracy, in my opinion the implication of the statement is that a test-intersection query will be formulated. The remaining argument is based on a brief article at a link that is now stale (http://astronomy.swin.edu.au/~pbourke/geometry/sphereline/). This article, as well as its restatement in the book, strays from the argument, although the web page has a clear thesis statement: If the line containing the line segment does not intersect the sphere, then the line segment does not intersect the sphere .

Indeed, this is a test-intersection query that attempts to quickly determine there is no intersection. The ray-sphere test-intersection query is an attempt at showing there is an intersection. The essence of the argument given at the web page and in the book is:
  • Compute the point Q on the line that is closest to C.
  • If Q is on the line segment and Q is inside the sphere, then attempt to compute the intersection points using the quadratic formula
The claim is that this can lead to a more efficient intersection calculation. I disagree. The fact that Q is inside the sphere does not guarantee that the line segment intersects the sphere (unless you treat the sphere as a solid). You expend some cycles computing Q and determining it is inside the sphere, only to have to compute the quadratic formula, something you could have done from the start. Moreover, the argument does not tell you what to do if Q is not on the line segment.

Maybe the intent of the article/book was the following:
  • If one segment end point is inside the sphere and the other end point is outside the sphere, the line segment intersects the sphere. Compute the intersection (if desired) by using the quadratic formula.
  • If both segment end points are inside the sphere, then the line segment does not intersect the sphere.
  • If both segment end points are outside the sphere, the line segment intersects the sphere if and only if the distance L from C to the line segment satisfies L <= R. Compute the intersection (if desired) by using the quadratic formula.
At any rate, if you want to adhere to the thesis statement of the article, you compute the squared distance L^2 from the line to C. If L^2 > R^2, there is no intersection between the line segment and the sphere. The construction of L^2 uses only additions, multiplications, and a division (no sqrt calls). You can avoid the division. If L = a/b, then the test for no intersection is a^2 > (b*R)^2.

Section 11.3.3, page 505. The displayed equation for the quadratic coefficient c is missing a subtraction of the squared radius. It should be
  c = k*(Px-Cx)^2 + ell*(Py-Cy)^2 + m*(Pz-Cz)^2 - r^2
The pseudocode at the bottom of the page is also missing the r^2, and the scaling factors need to be members of ellipsoid e. It should be
  c = e.k*(ell.p.x - e.center.x)^2 + e.ell*(ell.p.y - e.center.y)^2 +
      e.m*(ell.p.z - e.center.z)^2 - e.r^2;
The pseudocode uses 'l' (lower-case ell) for variable names, making it difficult to distinguish between '1' (one) and 'l' (ell). The standard axis-aligned representation for an ellipsoid does not separate out the radius term. The standard form is
  ((x-Cx)/Lx)^2 + ((y-Cy)/Ly)^2 + ((z-Cz)/Lz)^2 = 1
where the center is (Cx,Cy,Cz) and Lx, Ly, and Lz are the lengths measured along the axes from the center to the extreme points of the ellipsoid in those directions. I prefer more verbose names for the member names. The line has an origin and a direction. The ellipsoid has a center and length values. Another pass at the pseudocode is given below.
  int LineEllipsoidIntersect (Ellipsoid E, Line3 L, float t[2])
  {
      Vector3 diff = L.origin - E.center;
      float invLx = 1/E.Lx;
      float invLy = 1/E.Ly;
      float invLz = 1/E.Lz;
      float tmp1x = L.direction.x * invLx;
      float tmp1y = L.direction.y * invLy;
      float tmp1z = L.direction.z * invLz;
      float tmp2x = diff.x * invLx;
      float tmp2y = diff.y * invLy;
      float tmp2z = diff.z * invLz;
      float a = tmp1x * tmp1x + tmp1y * tmp1y + tmp1z * tmp1z;
      float b = 2 * (tmp1x * tmp2x + tmp1y * tmp2y + tmp1z * tmp2z);
      float c = tmp2x * tmp2x + tmp2y * tmp2y + tmp2z * tmp2z - 1;
      float d = b * b - 4 * a * c;
      if ( d > 0 )
      {
          float invTwoA = 1/(2*a);
          float rootD = sqrt(d);
          t[0] = (-b + rootD)*invTwoA;
          t[1] = (-b - rootD)*invTwoA;
          return 2;
      }
      else if ( discrm == 0 )
      {
          t[0] = -b/(2*a);
          return 1;
      }
      return 0;  // no solutions
  }
The return value of the function is the number of solutions. If t[i] is valid (0 <= i < numberOfSolutions), then the corresponding point of intersection is: L.origin + t[i]*L.direction.

Section 11.3.4, page 510. The solution for the quadratic formula is missing a minus sign on the "b" term. It should be
  t = (-b +/- sqrt(b^2-4*a*c))/(2*a)
The pseudocode has
  a = (tline.direction.x + tline.direction.y)^2;
This should be
  a = tline.direction.x^2 + tline.direction.y^2;

Section 11.3.4, page 511. Lines 4 and 5 are missing right parentheses. The pseudocode should be
  t[0] = (-b + sqrt(discrm))/(2*a);
  t[1] = (-b - sqrt(discrm))/(2*a);

Section 11.3.5, page 518. The page has 7 references to array elements "[3]". These should be "[2]".

Section 11.5.1, page 531. The pseudocode for b at the bottom of the page should be
  b = (s1 * n1n2dot - s2 * n1normsqr) / (n1n2dot^2 - n1normsqr * n2normsqr);
That is, the first occurrence of n2normsqr in the book should have been n1normsqr.

Section 11.5.2, page 534. The equations for Cramer's rule have errors. The inverse determinant should be
  invDet = 1/(a0*BC - b0*AC + c0*AB)
The Z value should be
  Z = (b0*AD + a0*DB - d0*AB)*invDet
I do not have a copy of Bowyer and Woodwark, the reference from which the equations are supposedly taken, so I do not know if these errors occur in that reference or were introduced in transcribing.

Section 11.5.3, pages 535-536. On page 535, the pseudocode for intersection of triangle and plane has a parameter "p2" in the function signature; it should not be there. On page 536, the pseudocode has some test for points either all above the plane or all below the plane. These tests are incorrect. In the style of the remainder of the pseudocode, the tests are simply stated as
  if ( d1 > 0.0 && d2 > 0.0 && d3 > 0.0 )
  {
      // all points above plane
      return false;
  }
  if ( d1 < 0.0 && d2 < 0.0 && d3 < 0.0 )
  {
      // all points below plane
      return false;
  }

Section 11.7.1, page 548. In matrix M, entry 6 of row 2 says "2(ux vy + vx vy)". This should be "2(ux vy + vx uy)".

Section 11.7.2, page 548. Equation (11.20) has "r" when it should be "r^2" (r squared).

Section 11.7.2, page 550. To be consistent with the other pseudocode, "if (b < r)" should be written as "if (b < sphere.radius)".

Section 11.7.2, page 550. The last line of pseudocode on this page should be
    float radius = sqrt(sphere.radius^2 - b^2);
The pseudocode also rejects the "just touching" case as "not an intersection"; that is, "b < sphere.radius" should be "b ≤ sphere.radius". In fact, since the pseudocode allows a small amount of error ("radius < epsilon"), it might even be good to support an error term in the first test: "if (b ≤ sphere.radius + epsilon)".

Section 11.7.2, page 550. The pseudocode computes b as an absolute value, but the absolute value should only be used in the comparison to the sphere radius on the line after computing b. The pseudocode on this page has had other corrections, so the revision is listed next. I do not like the book's formulation--the revision is more like how I write the actual code. This is the mathematically correct formulation, without concern for "epsilon" arithmetic to deal with numerical round-off errors. The plane normal is assumed to be a unit-length vector. The Circle3 data structure is for a circle in three dimensions.
    bool PlaneSphereIntersection (Plane plane, Sphere sphere, Circle3& circle)
    {
        float b = dotProd(plane.normal, sphere.center - plane.pointOnPlane);
        float absB = fabs(b);
        circle.center = sphere.center - b * plane.normal;
        circle.normal = plane.normal;
        if (absB <= sphere.radius)
        {
            // The sphere intersects the plane in a circle.  The circle is
            // degenerate when fabs(b) is equal to sphere.radius, in which
            // case the circle radius is zero.
            circle.radius = sqrt(sphere.radius * sphere.radius - b * b);
            return true;
        }
        return false;
    }
The circle parameter is an output of the function. It is valid only when the function's return value is true. The plane's normal vector is copied to the circle's normal so that you can actually reconstruct the circle. If the circle center is K, the normal is N, and the radius is R, then choose vectors U and V so that {U,V,N} is an orthonormal set (vectors are unit length and mutually perpendicular). A circle parameterization is
  X(t) = K + R*(cos(t)*U + sin(t)*V),  0 <= t < 2*pi
A numerically robust way to compute U and V is
    if (fabs(N.x) >= fabs(N.y))
    {
        // N.x or N.z is the largest magnitude component, swap them.
        invLength = 1/sqrt(N.x * N.x + N.z * N.z);
        U = Vector(-N.z * invLength, 0, N.x * invLength);
        V = Vector(N.y * U.z, N.z * U.x - N.x * U.z, -N.y * U.x);  // = Cross(N,U)
    }
    else
    {
        // N.y or N.z is the largest magnitude component, swap them.
        invLength = 1/sqrt(N.y * N.y + N.z * N.z);
        U = Vector(0, N.z * invLength, -N.y * invLength);
        V = Vector(N.y * U.z - N.z * U.y, -N.x * U.z, N.x * U.y);  // = Cross(N,U)
    }
If you want to allow for numerical round-off error, say, by requiring the intersection to be a single point when b is within epsilon of the radius, use
    bool PlaneSphereIntersection (Plane plane, Sphere sphere, Circle3& circle)
    {
        float b = dotProd(plane.normal, sphere.center - plane.pointOnPlane);
        float absB = fabs(b);
        circle.center = sphere.center - b * plane.normal;
        circle.normal = plane.normal;
        if (absB <= sphere.radius - epsilon)
        {
            // The sphere intersects the plane in a circle.  The circle is
            // degenerate when fabs(b) is equal to sphere.radius, in which
            // case the circle radius is zero.
            circle.radius = sqrt(sphere.radius * sphere.radius - b * b);
            return true;
        }
        else if (absB <= sphere.radius + epsilon)
        {
            // The sphere is _nearly_ tangent to the plane.  Choose the
            // intersection to be the point of tangency.
            circle.radius = 0;
            return true;
        }
        return false;
    }
The magnitude of the nonnegative number epsilon is up to the programmer. Observe that if you choose epsilon to be zero, you effectively have the pseudocode for the mathematically correct case.

Section 11.7.3, page 562. The pseudocode has a few issues. The first line of code that computes d should be a signed distance calculation. The first comparison of cos(theta) should have the abs about the cosine term, not about the epsilon term. If an epsilon comparison is made for when cos(theta) is nearly zero, to be consistent an epsilon comparison should also be made for when cos(theta) is nearly one (the circle case). Here is my own rewrite of the pseudocode without epsilon tests. The following data structures are assumed.
    struct Plane    { Point origin; Vector normal; };
    struct Cylinder { Point origin; Vector direction; float radius; };
    struct Line     { Point origin; Vector direction; };
    struct Circle   { Point center; Vector normal; float radius; };
    struct Ellipse  { Point center; Vector normal, major, minor; float majorLength, minorLength; };
    
    enum Type
    {
        EMPTY_SET,
        ONE_LINE,
        TWO_LINES,
        CIRCLE,
        ELLIPSE
    };
    struct Intersection
    {
        Type type;
        Line line;          // valid when type = ONE_LINE
        Line line1, line2;  // valid when type = TWO_LINES
        Circle circle;      // valid when type = CIRCLE
        Ellipse ellipse;    // valid when type = ELLIPSE
    };
The plane normal, cylinder direction, line direction, circle normal, ellipse normal, ellipse major direction, and ellipse minor direction are all required to be unit-length vectors.

The output is intr. The returned bool value is true when there is an intersection, in which case intr is valid (of the type specified in its Type member). The returned value is false when there is no intersection, in which case the Type member has value EMPTY_SET.
    bool PlaneCylinderIntersection (Plane plane, Cylinder cylinder, Intersection& intr)
    {
        float d = Dot(plane.normal, cylinder.origin - plane.origin);
        Point projectedCylinderOrigin = cylinder.origin - d * plane.normal;
        float cosTheta = Dot(cylinder.direction, plane.normal);
        if (abs(cosTheta) > 0)  // cylinder axis intersects plane in a unique point
        {
            if (abs(cosTheta) < 1)  // intersection is an ellipse
            {
                intr.type = ELLIPSE;
                intr.ellipse.normal = plane.normal;
                intr.ellipse.center = projectedCylinderOrigin - (d/cosTheta) * cylinder.direction;
                intr.ellipse.major = cylinder.direction - cosTheta * plane.normal;
                intr.ellipse.minor = Cross(plane.normal, intr.ellipse.major);
                intr.ellipse.majorLength = cylinder.radius / abs(cosTheta);
                intr.ellipse.minorLength = cylinder.radius;
                Normalize(intr.ellipse.major);
                Normalize(intr.ellipse.minor);
                return true;
            }
            else  // abs(cosTheta) == 1,  intersection is a circle
            {
                intr.type = CIRCLE;
                intr.circle.normal = plane.normal;
                intr.circle.center = projectedCylinderOrigin;
                intr.circle.radius = cylinder.radius;
                return true;
            }
        }
        else  // abs(cosTheta) == 0,  cylinder is parallel to plane
        {
            if (abs(d) < cylinder.radius)  // cylinder intersects plane in two lines
            {
                Vector offset = Cross(cylinder.direction, plane.normal);
                float e = sqrt(cylinder.radius * cylinder.radius - d * d);
                intr.type = TWO_LINES;
                intr.line1.origin = projectedCylinderOrigin - e * offset;
                intr.line1.direction = cylinder.direction;
                intr.line2.origin = projectedCylinderOrigin + e * offset;
                intr.line2.direction = cylinder.direction;
                return true;
            }
            else if (abs(d) == cylinder.radius)  // cylinder intersects plane in one line
            {
                intr.type = ONE_LINE;
                intr.line.origin = projectedCylinderOrigin;
                intr.line.direction = cylinder.direction;
                return true;
            }
            else // abs(d) > cylinder.radius,  cylinder does not intersect plane
            {
                intr.type = EMPTY_SET;
                return false;
            }
        }
    }

Section 11.7.4, pages 564-569. The first paragraph of Intersection Detection mentions a base point B and refers you to Figure 11.14, which shows the base point labeled as C. Later in Intersection Detection, there is a section Finite Cone that refers to the base point using name B. This point is as shown in Figure 11.14 and is not the same as the apex A.

On page 567, a vector hat(w) is defined. The author used the hat notation to denote unit-length vectors. However, the vector Cross(hat(d),Cross(hat(n),hat(d))) is not necessarily unit length. For example, if hat(n) is nearly parallel to hat(d), Cross(hat(n),hat(d)) is nearly the zero vector, in which case hat(w) is nearly the zero vector. The correct equation is
    hat(w) = Cross(hat(d),Cross(hat(n),hat(d)))/Length(Cross(hat(n),hat(d)))
On pages 567, 568, and 569, the author has occurrences of arrow(w), but these should be hat(w).

The variable a is the distance between points B and Ia, so the last equation on page 567 should be
    a = Length(Ia - B)
That is, the subtraction of h is incorrect. Just to be clear, the variable b is the distance between points B and Ic, and the variable c is the distance between points Ia and Ic. On pages 568 and 569, any occurrence of Length(Ia-B)-h should be replaced by Length(Ia-B).

Observe that with the correct hat(w) equation, it must be that Dot(hat(n),hat(w)) = Length(Cross(hat(n),hat(d)), so it is not necessary to explicitly compute hat(w) in the source code. The equation for b2 is
  b2 = a2 [1/Length(Cross(hat(n),hat(d)))2 - 1]
However, notice that this equation has numerical problems when hat(n) and hat(d) are nearly parallel. The length of the cross product is nearly zero, so computing the reciprocal length is a problem. Instead, you want to code up the test as
  v = Cross(hat(n),hat(d));
  lensqr = Dot(v,v);
  a*a*(1 - lensqr) ≤ lensqr*r*r;
Here is an alternative for the intersection test of plane and finite cone. Let the plane have unit-length normal N and contain point P. The plane equation is Dot(N,X-P) = 0. The idea is to compute the interval of projection of the finite cone onto the line containing P with direction N. This inteval has endpoints P + t0*N and P + t1*N, where t0 < t1. If the interval [t0,t1] contains 0, then the cone and plane intersect.

The projection interval endpoints are determined by one of the following.
  1. C and two points on the circle that terminates the cone wall. The two points project to the same interval endpoint. This happens when N and A are parallel.

  2. Two points on the circle.

  3. C and one point on the circle.
The cone has vertex C, unit-length axis direction D, height h, and angle A measured from the axis to the cone wall. The circle that terminates the cone wall is parameterized by
    X = C + h*D + h*tan(A)*(cos(B)*U + sin(B)*V)
for B in [0,2*pi). The vectors U and V are unit length and {U,V,D} is an orthonormal set (vectors are unit length and mutually perpendicular). The projection of X onto the normal line is P + t(B)*N, where
   t(B) = Dot(N,X-P) = Dot(N,C-P) + h*Dot(N,D) + h*tan(A)*(cos(B)*Dot(N,U) + sin(B)*Dot(N,V))
I have indicated that t is a function of B (different circle points project to different normal line points). The projection of C onto the normal line is Dot(N,C-P).

The minimum value of t(B) is t0 and the maximum value of t(B) is t1. Define d = Dot(N,C-P) and
 r = sqrt(Dot(N,U)^2 + Dot(N,V)^2) = sqrt(1 - Dot(N,D)^2)
t0 is the minimum, and t0 is the maximum, of the set
    {d, d + h*Dot(N,D) - h*tan(A)*r, d + h*Dot(N,D) + h*tan(A)*r}
The pseudocode is
    bool ConeIntersectsPlane (Cone cone, Plane plane)
    {
        Vector CmP = cone.C - plane.P;
        float d = Dot(plane.N, CmP);
        float NdotD = Dot(plane.N, cone.D);
        float r = sqrt(1 - NdotD*NdotD);
        float tmp0 = d + cone.h*NdotD;
        float tmp1 = cone.h*cone.tanAngle*r;
        float vm = tmp0 - tmp1;
        float vp = tmp0 + tmp1;
        float t0 = min(d, min(vm, vp));
        float t1 = max(d, max(vm, vp));
        return t0 <= 0 && 0 <= t1;
    }
This algorithm does not require you to explicitly determine which of the book-mentioned cases 1, 2, 3a, or 3b applies.

Section 11.7.4, page 571. The end of the last sentence of the first paragraph of the subsection Nondegenerate Plane-Cone Intersections should be: cos(theta) = Dot(hat(a),hat(n)).

Section 11.7.4, page 575. The style of the pseudocode differs from other pages in this section. In the first block starting with ``float d = '', remove the 'float' types (to be consistent). In the section about the parabola, consistent notation would be
    parabola.axis = Normalize(cone.axis - cosTheta * plane.normal);
    parabola.perpendicularAxis = Cross(parabola.axis, plane.normal);
    parabola.vertex = cone.vertex + d* plane.normal + e*parabola.axis;
    parabola.focalLength = d/2 * cosTheta;
Also on this page, the first paragraph has a condition on the dot product of the plane normal hat(n) and the cone axis hat(a) in order that the two vectors be parallel. The condition should be: |Dot(hat(n),hat(a)| = 1 (for the true mathematical classification) or |Dot(hat(n),hat(a))| >= 1-epsilon (if you want to worry about numerical round-off errors).

Section 11.7.4, page 580. The pseudocode has a line that computes the minor axis. It should be
    minorAxis = Cross(plane.normal, majorAxis);
because nowhere is there defined ``ellipse.u''.

Section 11.11.1, page 612. The last line of pseudocode is
    if (value < min) min = value; else max = value;
This should be
    if (value < min) min = value; else if (value > max) max = value;

Section 11.12, page 628-629. The pseudocode on page 628, the case where "Ray parallel to planes" has a test "P.x < box.xMax". That should be "P.x > box.xMax". The pseudocode on page 629 for swapping t0 and t1 is incorrect. It should be "tmp = t1; t1 = t0; t0 = tmp;"

Section 11.12, pages 628-629. The pseudocode for ray-AABB intersection has a loop that illustrates clipping only against the planes x=xmin and x=xmax. This makes it unclear how the the full set of tests fit in with the rest of the loop logic. The following pseudocode is the replacement. The Point and Vector data structures are assumed to have operator[] functions to access the coordinates. Index 0 corresponds to x (P[0] and P.x are equivalent), index 1 corresponds to y, and index 2 corresponds to z. It is my opinion that this function should also return the interval of intersection (first point of intersection and last point of intersection). This change shows up in the last input to the function. It is guaranteed that tIntersect[0] <= tIntersect[1]. If the ray just grazes the AABB, we will return the same value in both parameters.
  boolean RayIntersectAABB (Point P, Vector d, AABB box, float tIntersect[2])
  {
      // The first part of this code tests for line-AABB intersection.
      tNear = -INFINITY;
      tFar = INFINITY;
      for (i = 0; i < 3; i++)
      {
          // Check if the line is parallel to the min[i] and max[i] planes.
          if (abs(d[i]) < epsilon)
          {
              // The line is parallel to the min[i] and max[i] planes.
              if (P[i] < box.min[i] || P[i] > box.max[i])
              {
                  // The line is outside the slab defined by min[i] and max[i].
                  return false;
              }
              // else:  The line is inside the slab defined by min[i] and max[i].
              // Just continue testing with the next slab.
              continue;
          }
          // The line is not parallel to the slab planes, so compute the
          // parameters of intersection.  The line segment of intersection is
          // P+t*d, where t0 <= t <= t1.
          t0 = (box.min[i] - P[i])/d[i];
          t1 = (box.max[i] - P[i])/d[i];
          if (t0 > t1)
          {
              tmp = t1;
              t1 = t0;
              t0 = tmp;
          }
          
          // Compare with current values.  The current slab may have increased
          // tNear and/or decreased tFar.
          if (t0 > tNear)
          {
              tNear = t0;
          }
          if (t1 < tFar)
          {
              tFar = t1;
          }
          
          // Check if the line misses the AABB entirely.
          if (tNear > tFar)
          {
              return false;
          }
      }
      
      // At this place in the code, the line intersects the AABB in the line
      // segment P+t*d, tNear <= t <= tFar.  If all you care about is line-AABB
      // intersection, then use the following three lines of code and omit the
      // remaining pseudocode.
      //
      // tIntersect[0] = tNear;
      // tIntersect[1] = tFar;
      // return true;
      // The ray is (P+t*d, t >= 0).  We need to compute the intersection of
      // the interval tNear <= t <= tFar with the interval t >= 0.
      if (tFar < 0)
      {
          return false;
      }
      if (tNear < 0)
      {
          tNear = 0;
      }
      tIntersect[0] = tNear;
      tIntersect[1] = tFar;
      return true;
  }

Section 11.12.3, page 633. The pseudocode for the ray-OBB intersection is absurd. The first line of code is missing two right parentheses after "box.axis[i]". The inequality "-r + box.halfLength[i] > 0" should be "-r + box.halfLength[i] < 0". The rest of the pseudocode is worthless. The value s can very well be zero, leading to a division by zero when computing t0 and t1. If you find worthless pseudocode in the book, first let me know and, second, look at the code at my website. I wrote the real code (not the pseudocode), so if you find errors in the real code, please let me know.

Section 11.12.4, page 635. The return statement in the pseudocode is incorrect. It appears to be designed for culling an AABB against the plane rather than indicating if the AABB is intersected by the plane. The return statement also indicates that the plane and AABB are considered not to be intersecting if the box is just touching the plane at a single vertex or single edge or single face. If the plane intersects the box at any point, I suggest replacing the return statement by
  return Dot(plane.normal,dMin)+plane.d <= 0 && Dot(plane.normal,dMax)+plane.d >= 0;
I much prefer the "test intersection" that I use for plane-OBB, which specializes easily when you have an AABB. No branching, just a few arithmetic calculations. Let N = (nx,ny,nz) be the plane normal and d be the plane constant. The plane is Dot(N,X)+d = 0. Let the AABB minima be x0, y0, and z0. Let the maxima be x1, y1, and z1. The test for intersection is
  |nx*(x1+x0)+ny*(y1+y0)+nz*(z1+z0)+2*d| <= |nx*(x1-x0)|+|ny*(y1-y0)|+|nx*(z1-z0)|

Section 11.12.7, pages 641 and 642. The last line of page 641 has a lower-case subscript 'a' on the unit vector u. That should be an upper-case 'A'. The last line of page 642 also has a lower-case subscript of 'a', but it should be an upper-case 'B'.

Section 11.12.8, page 645. The upper bounds on the displayed inequalities are wrong. They should be
  P_{min,x} <= Q_x <= P_{max,x}
  P_{min,y} <= Q_y <= P_{max,y}
  P_{min,z} <= Q_z <= P_{max,z}
The pseudocode on that same page is inconsistent about its reference to the sphere center. I changed "sphere.c" to "sphere.center" in two places.

Section 11.12.10, page 660. The equation for c0 is missing a power 2. It should be
    c0 = (Dot(P,P) - (R^2+r^2))^2 - 4*R^2*(r^2 - P_z^2)

Section 11.12, pages 661-662. Yet another badly written topic.
  • Throughout the discussion, the torus point X is written as X = [Xx   Xy   Xz]. This is awkward notation when the torus point is simply (x,y,z).
  • On page 661, "cutting with the XZ plane" should be "cutting with the XY plane".
  • On page 661, the displayed equation for u has two occurrences of "arccos(θ)" that should be "arccos(Xx/ru)".
  • On page 661, "...if we want the parametric origin to be on the 'inside' of the torus...". In more mathematical terms: Given a circular cross section of radius r for the torus, we want the circle point closest to the origin to correspond to parameter v = 0.
  • On page 662, it is not clear in Figure 11.77 which angle is labeled with φ. Turns out that the right triangle has φ corresponding to the circle center. This angle starts at 0 for the circle point at distance R-r from the origin. It increase to π as you walk the circle counterclockwise to the circle point at distance R+r from the origin.
  • On page 662, the length rv is the same as length ru. Perhaps the distance from the origin to the projected point should have been named something simple like d.
  • On page 662, the displayed equation has "sin(φ) = Xz/rv" but should be "sin(φ) = Xz/r".
  • On page 662, the displayed equation specifies parameter u but should be parameter v. There are two occurrences of "arccos(φ)" that should be "arccos(-(rv-R)/r)".
  • The paragraphs starting with "The partial derivatives..." are mathematical gibberish. After a discussion of computing functions u(x,y,z) and v(x,y,z), one might assume that the partial derivatives of u and v with respect to x, y and z were to be computed. This is not the case. The notation ∂u and ∂v is not standard. The vector (X'-O) lives in 3D, so what does (X'-O) mean? What does it mean "the partial derivatives are computed in the local space of the torus"? The intent of the paragraphs is to motivate computing tangent vectors to the torus, presumably as a parametric surface P(u,v) =(x(u,v),y(u,v),z(u,v)). Linearly independent tangent vectors are ∂P/∂u and ∂P/∂v. Computing them requires knowing how x, y and z depend on u and v, but the previous discussion shows only how to compute u and v as functions of x, y and z.
A simpler presentation is the following.

A torus with origin (0,0,0), outer radius r0 and inner radius r1 (with r0 >= r1) is defined implicitly as follows. The point P0 = (x,y,z) is on the torus. Its projection onto the xy-plane is P1 = (x,y,0). The circular cross section of the torus that contains the projection has radius r0 and center P2 = r0*(x,y,0)/sqrt(x2+y2). The points triangle <P0,P1,P2> is a right triangle with right angle at P1. The hypotenuse <P0,P2> has length r1, leg <P1,P2> has length z and leg <P0,P1> has length |r0 - sqrt(x^2+y^2)|. The Pythagorean theorem says z2 + |r0 - sqrt(x2+y2)|2 = r12. This can be algebraically manipulated to (x2 + y2 + z2 + r02 - r12)2 - 4 * r02 * (x2 + y2) = 0 A parametric form is x = (r0 + r1 * cos(v)) * cos(u), y = (r0 + r1 * cos(v)) * sin(u), z = r1 * sin(v) for u in [0,2*pi) and v in [0,2*pi). The position is P = (x,y,z) and the partial derivatives are ∂P/∂u = (r0 + r1 * cos(v)) * (-sin(u), cos(u), 0) and ∂P/∂v = r1(-cos(u) * sin(v), -sin(u) * sin(v), cos(v)). It is simple to show that the partial derivatives are perpendicular. To compute (u,v) from (x,y,z), u = atan2(y,x) and v = atan2(z,sqrt(x2+y2)-r0).

Section 13.1.4, page 686. The pseudocode has two sets named CoPos and CoNeg. These should be CoSame and CoDiff, respectively.

Section 13.3.2, page 698. The pseudocode for PointInSubpolygon has a call to GetMiddleIndex that is missing its last parameter. The call should be
  mid = GetMiddleIndex(i0, i1, C.N);

Section 13.3.3, page 705. The second row of Table 13.1 has the four "Update" entries "i 0 m p". These entries should be "i 0 p m".

Section 13.8.1, page 758-759. The input N is number of points, but the first mesh insert call returns a new triangle also named N. These are, of course, different things. For clarity, use Nj as the return value of the two "mesh.Insert" calls, one on page 758 and one on page 759. Use Nj as the argument to the two "stack.Push" calls, one on page 758 and one on page 759.

In the bottom block of pseudocode on page 759, the figure reference should be 13.46, not 13.45. After that line the pseudocode should be
  compute i0, i1, i2, i3 wwith T.v(i1) = A.v(i1) and T.v(i2) = A.v(i2);
In the "if" clause, two lines must be added to the end of the clause
  mesh.Remove(T);
  mesh.Remove(A);
After the two new triangles N0 and N1 are inserted into the mesh, the old triangles T and A must be removed.

Section 13.9.1, pages 767-768. The visibility graph does not contain an arc (i,i), where i is an index of a polygon vertex. The text should also state that adjacent polygon vertices are not considered to be visible to each other. Therefore, the visibility graph does not contain arcs (i,(i+1) modulo n), where n is the number of polygon vertices and 0 <= i <= n-1. On page 768, the text states that the adjacency matrix for a convex polygon has 3 zero diagonals, the rest of the entries are ones. However, the entries at (row,column) of (0,n-1) and (n-1,0) are also zero since these correspond to adjacent polygon vertices.

Section 13.9.1, page 769. The sentence that starts with "If < V_{i-1}, V_{i+2} > is a diagonal" should have V_{i+1}, not V_{i+2}.

Section 13.9.1, page 771. In the SegmentInCode listing, the case for the "vertex is convex" has the inequalities reversed. The line of code should be
                      return Kross(diff, edgeR) < 0 and Kross(diff, edgeL) > 0
                  

Section 13.9.3, page 777. The first sentence of the last paragraph has a subject-verb mismatch. It should be "An edge sk has ...". In the same paragraph, the end of the second line has "s1.y" but should be "s1.y0".

Appendix A.1, page 831. The fourth line from the top of the page is "(1) x1 - 3*x2 + 2*x3 = 6". It should be "(1) x1 - 3*x2 - 2*x3 = 6". The last line of the page lists the solution as (1,-3,-2). It should be (1,-3,2).

Appendix A.1, page 831. The "09 January 2003" entry mentions that the first displayed equation "(1) x1-3*x2+2*x3=6" should be "(1) x1-3*x2-2x3=6". The same correction must be made in the next occurrence of this equation near the bottom of the page (after "which gives us a new system").

Appendix A.1, page 832. The pseudocode at the top of the page has a loop "(for k = 1 to n)". This should be "for (k = j+1 to n)". In that loop after the multiplier "m" is computed, you could set a_{k,j} = 0 to make it clear you are eliminating the column entries. If an implementation allows for the row-reduced matrix to be used in further calculations, the assignment of zero should occur in the code. But if the only goal is to solve the system, an implementation need not assign zero since a_{k,j} is never used again in the code.

The first line of pseudocode for back substitution has "x_n = c_n/a_{n,h}". This should be "x_n = c_n/a_{n,n}".

Appendix A.1, page 832. In the elimination step, the line of code
    a_{k,i} = a_{k,i} + m * a_{k,j}
should be
    a_{k,i} = a_{k,i} + m * a_{j,i}

Appendix A.2, page 834. In the paragraph "The row reduction..." is "The condition for existence of solutions is...". The equations should be corrected to
  0 = a0*b1 - a1*b0
    = (c0+c1*y)*b1 - a1*(d0+d1*y)
    = (b1*c0-a1*d0) + (b1*c1-a1*d1)*y
Later is "as long as b1*c1-b0*d1" which should be "as long as b1*c1-a1*d1". Finally, the last equation of the last displayed equation array should be
  = (c0*f0-e0*d0)+((c0*f1-e0*d1)+(c1*f0-e1*d0))*y + (c1*f1-e1*d1)*y^2

Appendix A.2, page 838. The last set of equations with the d[i][j] should have alpha[2][2] replaced by alpha[2][0] and beta[2][2] replaced by beta[2][0].

Appendix A.2, pages 844 and 846. The second line from the bottom on page 844 should have "e2 = -98". The fifth line from the bottom of page 846 should have "(1.11596,1.14945,1.63258)".

Appendix A.4.3, page 861. A bunch of errors are in "Quaternion to Matrix" and "Matrix to Quaternion". This material was from "3D Game Engine Design", where the same errors occur, but I had posted the correction for that book in February 2001. If R = [rij], then equation (A.8) should have
  r00 = 1-2y^2-2z^2
  r01 = 2xy-2wz
  r02 = 2xz+2wy
  r10 = 2xy+2wz
  r11 = 1-2x^2-2z^2
  r12 = 2yz-2wx
  r20 = 2xz-2wy
  r21 = 2yz+2wx
  r22 = 1-2x^2-2y^2
The last sentence of the first paragraph of "Matrix to Quaternion" has
  y = (r20-r02)/(4w)
but it should be
  y = (r02-r20)/(4w)
In the second paragraph of that section. When r00 is the maximum, the book has
  w = (r12-r21)/(4x)
but it should be
  w = (r21-r12)/(4x)
When r11 is the maximum, the book has
  w = (r20-r02)/(4y)
but it should be
  w = (r02-r20)/(4y)
When r22 is the maximum, the book has
  w = (r01-r10)/(4z)
but it should be
  w = (r10-r01)/(4z)

Appendix A.5.1, page 872. The first displayed equation is
  max{|a_0|, 1+|a_1|, ..., 1+|a_{n-1}|} = 1 + max{|a_0|, ..., |a_{n-1}|}
but the last equaltity should be an inequality
  max{|a_0|, 1+|a_1|, ..., 1+|a_{n-1}|} ≤ 1 + max{|a_0|, ..., |a_{n-1}|}

Appendix A.7.7, page 888. The last displayed equation is
    E(hat(c)) = (sum_{i=0}^n hat(c) * vector(v))^2 = hat(c)^T M hat(c)
which is the square of a sum, but it should be a sum of squares,
    E(hat(c)) = sum_{i=0}^n (hat(c) * vector(v))^2 = hat(c)^T M hat(c)

Appendix A.7.8, page 889. Line 5 of the first paragraph of section A.7.8 has
    M = sum_{i=0}^2 vector(v) vector(v)^T
The upper limit should be n,
    M = sum_{i=0}^n vector(v) vector(v)^T

Appendix B.3.3, page 941. All the product formulas are incorrect. They should be
                      sin(a)*sin(b) = (cos(a-b) - cos(a+b))/2
                      sin(a)*cos(b) = (sin(a+b) + sin(a-b))/2
                      cos(a)*cos(b) = (cos(a-b) + cos(a+b))/2
                  
The last formula in the book appears to be an attempt to cover all the cases, but the first three equations are all the cases.