michael@0: /* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
michael@0: /* This Source Code Form is subject to the terms of the Mozilla Public
michael@0:  * License, v. 2.0. If a copy of the MPL was not distributed with this
michael@0:  * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
michael@0: 
michael@0: #include "nsSMILKeySpline.h"
michael@0: #include <stdint.h>
michael@0: #include <math.h>
michael@0: 
michael@0: #define NEWTON_ITERATIONS          4
michael@0: #define NEWTON_MIN_SLOPE           0.02
michael@0: #define SUBDIVISION_PRECISION      0.0000001
michael@0: #define SUBDIVISION_MAX_ITERATIONS 10
michael@0: 
michael@0: const double nsSMILKeySpline::kSampleStepSize =
michael@0:                                         1.0 / double(kSplineTableSize - 1);
michael@0: 
michael@0: void
michael@0: nsSMILKeySpline::Init(double aX1,
michael@0:                       double aY1,
michael@0:                       double aX2,
michael@0:                       double aY2)
michael@0: {
michael@0:   mX1 = aX1;
michael@0:   mY1 = aY1;
michael@0:   mX2 = aX2;
michael@0:   mY2 = aY2;
michael@0: 
michael@0:   if (mX1 != mY1 || mX2 != mY2)
michael@0:     CalcSampleValues();
michael@0: }
michael@0: 
michael@0: double
michael@0: nsSMILKeySpline::GetSplineValue(double aX) const
michael@0: {
michael@0:   if (mX1 == mY1 && mX2 == mY2)
michael@0:     return aX;
michael@0: 
michael@0:   return CalcBezier(GetTForX(aX), mY1, mY2);
michael@0: }
michael@0: 
michael@0: void
michael@0: nsSMILKeySpline::GetSplineDerivativeValues(double aX, double& aDX, double& aDY) const
michael@0: {
michael@0:   double t = GetTForX(aX);
michael@0:   aDX = GetSlope(t, mX1, mX2);
michael@0:   aDY = GetSlope(t, mY1, mY2);
michael@0: }
michael@0: 
michael@0: void
michael@0: nsSMILKeySpline::CalcSampleValues()
michael@0: {
michael@0:   for (uint32_t i = 0; i < kSplineTableSize; ++i) {
michael@0:     mSampleValues[i] = CalcBezier(double(i) * kSampleStepSize, mX1, mX2);
michael@0:   }
michael@0: }
michael@0: 
michael@0: /*static*/ double
michael@0: nsSMILKeySpline::CalcBezier(double aT,
michael@0:                             double aA1,
michael@0:                             double aA2)
michael@0: {
michael@0:   // use Horner's scheme to evaluate the Bezier polynomial
michael@0:   return ((A(aA1, aA2)*aT + B(aA1, aA2))*aT + C(aA1))*aT;
michael@0: }
michael@0: 
michael@0: /*static*/ double
michael@0: nsSMILKeySpline::GetSlope(double aT,
michael@0:                           double aA1,
michael@0:                           double aA2)
michael@0: {
michael@0:   return 3.0 * A(aA1, aA2)*aT*aT + 2.0 * B(aA1, aA2) * aT + C(aA1);
michael@0: }
michael@0: 
michael@0: double
michael@0: nsSMILKeySpline::GetTForX(double aX) const
michael@0: {
michael@0:   // Find interval where t lies
michael@0:   double intervalStart = 0.0;
michael@0:   const double* currentSample = &mSampleValues[1];
michael@0:   const double* const lastSample = &mSampleValues[kSplineTableSize - 1];
michael@0:   for (; currentSample != lastSample && *currentSample <= aX;
michael@0:         ++currentSample) {
michael@0:     intervalStart += kSampleStepSize;
michael@0:   }
michael@0:   --currentSample; // t now lies between *currentSample and *currentSample+1
michael@0: 
michael@0:   // Interpolate to provide an initial guess for t
michael@0:   double dist = (aX - *currentSample) /
michael@0:                 (*(currentSample+1) - *currentSample);
michael@0:   double guessForT = intervalStart + dist * kSampleStepSize;
michael@0: 
michael@0:   // Check the slope to see what strategy to use. If the slope is too small
michael@0:   // Newton-Raphson iteration won't converge on a root so we use bisection
michael@0:   // instead.
michael@0:   double initialSlope = GetSlope(guessForT, mX1, mX2);
michael@0:   if (initialSlope >= NEWTON_MIN_SLOPE) {
michael@0:     return NewtonRaphsonIterate(aX, guessForT);
michael@0:   } else if (initialSlope == 0.0) {
michael@0:     return guessForT;
michael@0:   } else {
michael@0:     return BinarySubdivide(aX, intervalStart, intervalStart + kSampleStepSize);
michael@0:   }
michael@0: }
michael@0: 
michael@0: double
michael@0: nsSMILKeySpline::NewtonRaphsonIterate(double aX, double aGuessT) const
michael@0: {
michael@0:   // Refine guess with Newton-Raphson iteration
michael@0:   for (uint32_t i = 0; i < NEWTON_ITERATIONS; ++i) {
michael@0:     // We're trying to find where f(t) = aX,
michael@0:     // so we're actually looking for a root for: CalcBezier(t) - aX
michael@0:     double currentX = CalcBezier(aGuessT, mX1, mX2) - aX;
michael@0:     double currentSlope = GetSlope(aGuessT, mX1, mX2);
michael@0: 
michael@0:     if (currentSlope == 0.0)
michael@0:       return aGuessT;
michael@0: 
michael@0:     aGuessT -= currentX / currentSlope;
michael@0:   }
michael@0: 
michael@0:   return aGuessT;
michael@0: }
michael@0: 
michael@0: double
michael@0: nsSMILKeySpline::BinarySubdivide(double aX, double aA, double aB) const
michael@0: {
michael@0:   double currentX;
michael@0:   double currentT;
michael@0:   uint32_t i = 0;
michael@0: 
michael@0:   do
michael@0:   {
michael@0:     currentT = aA + (aB - aA) / 2.0;
michael@0:     currentX = CalcBezier(currentT, mX1, mX2) - aX;
michael@0: 
michael@0:     if (currentX > 0.0) {
michael@0:       aB = currentT;
michael@0:     } else {
michael@0:       aA = currentT;
michael@0:     }
michael@0:   } while (fabs(currentX) > SUBDIVISION_PRECISION
michael@0:            && ++i < SUBDIVISION_MAX_ITERATIONS);
michael@0: 
michael@0:   return currentT;
michael@0: }