Home | History | Annotate | Download | only in Intersection 
      1 /*
      2  * Copyright 2012 Google Inc.
      3  *
      4  * Use of this source code is governed by a BSD-style license that can be
      5  * found in the LICENSE file.
      6  */
      7 #include "DataTypes.h"
      8 #include "Extrema.h"
      9 
     10 static int validUnitDivide(double numer, double denom, double* ratio)
     11 {
     12     if (numer < 0) {
     13         numer = -numer;
     14         denom = -denom;
     15     }
     16     if (denom == 0 || numer == 0 || numer >= denom)
     17         return 0;
     18     double r = numer / denom;
     19     if (r == 0) { // catch underflow if numer <<<< denom
     20         return 0;
     21     }
     22     *ratio = r;
     23     return 1;
     24 }
     25 
     26 /** From Numerical Recipes in C.
     27 
     28     Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C])
     29     x1 = Q / A
     30     x2 = C / Q
     31 */
     32 static int findUnitQuadRoots(double A, double B, double C, double roots[2])
     33 {
     34     if (A == 0)
     35         return validUnitDivide(-C, B, roots);
     36 
     37     double* r = roots;
     38 
     39     double R = B*B - 4*A*C;
     40     if (R < 0) {  // complex roots
     41         return 0;
     42     }
     43     R = sqrt(R);
     44 
     45     double Q = (B < 0) ? -(B-R)/2 : -(B+R)/2;
     46     r += validUnitDivide(Q, A, r);
     47     r += validUnitDivide(C, Q, r);
     48     if (r - roots == 2 && AlmostEqualUlps(roots[0], roots[1])) { // nearly-equal?
     49         r -= 1; // skip the double root
     50     }
     51     return (int)(r - roots);
     52 }
     53 
     54 /** Cubic'(t) = At^2 + Bt + C, where
     55     A = 3(-a + 3(b - c) + d)
     56     B = 6(a - 2b + c)
     57     C = 3(b - a)
     58     Solve for t, keeping only those that fit between 0 < t < 1
     59 */
     60 int findExtrema(double a, double b, double c, double d, double tValues[2])
     61 {
     62     // we divide A,B,C by 3 to simplify
     63     double A = d - a + 3*(b - c);
     64     double B = 2*(a - b - b + c);
     65     double C = b - a;
     66 
     67     return findUnitQuadRoots(A, B, C, tValues);
     68 }
     69 
     70 /** Quad'(t) = At + B, where
     71     A = 2(a - 2b + c)
     72     B = 2(b - a)
     73     Solve for t, only if it fits between 0 < t < 1
     74 */
     75 int findExtrema(double a, double b, double c, double tValue[1])
     76 {
     77     /*  At + B == 0
     78         t = -B / A
     79     */
     80     return validUnitDivide(a - b, a - b - b + c, tValue);
     81 }
     82