// Copyright 2014 The Chromium Authors. All rights reserved.
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file.

#include "ui/gfx/geometry/cubic_bezier.h"

#include <algorithm>
#include <cmath>

#include "base/logging.h"

namespace gfx {

static const double kBezierEpsilon = 1e-7;

CubicBezier::CubicBezier(double p1x, double p1y, double p2x, double p2y) {
  InitCoefficients(p1x, p1y, p2x, p2y);
  InitGradients(p1x, p1y, p2x, p2y);
  InitRange(p1y, p2y);
}

CubicBezier::CubicBezier(const CubicBezier& other) = default;

void CubicBezier::InitCoefficients(double p1x,
                                   double p1y,
                                   double p2x,
                                   double p2y) {
  // Calculate the polynomial coefficients, implicit first and last control
  // points are (0,0) and (1,1).
  cx_ = 3.0 * p1x;
  bx_ = 3.0 * (p2x - p1x) - cx_;
  ax_ = 1.0 - cx_ - bx_;

  cy_ = 3.0 * p1y;
  by_ = 3.0 * (p2y - p1y) - cy_;
  ay_ = 1.0 - cy_ - by_;
}

void CubicBezier::InitGradients(double p1x,
                                double p1y,
                                double p2x,
                                double p2y) {
  // End-point gradients are used to calculate timing function results
  // outside the range [0, 1].
  //
  // There are three possibilities for the gradient at each end:
  // (1) the closest control point is not horizontally coincident with regard to
  //     (0, 0) or (1, 1). In this case the line between the end point and
  //     the control point is tangent to the bezier at the end point.
  // (2) the closest control point is coincident with the end point. In
  //     this case the line between the end point and the far control
  //     point is tangent to the bezier at the end point.
  // (3) the closest control point is horizontally coincident with the end
  //     point, but vertically distinct. In this case the gradient at the
  //     end point is Infinite. However, this causes issues when
  //     interpolating. As a result, we break down to a simple case of
  //     0 gradient under these conditions.

  if (p1x > 0)
    start_gradient_ = p1y / p1x;
  else if (!p1y && p2x > 0)
    start_gradient_ = p2y / p2x;
  else
    start_gradient_ = 0;

  if (p2x < 1)
    end_gradient_ = (p2y - 1) / (p2x - 1);
  else if (p2x == 1 && p1x < 1)
    end_gradient_ = (p1y - 1) / (p1x - 1);
  else
    end_gradient_ = 0;
}

// This works by taking taking the derivative of the cubic bezier, on the y
// axis. We can then solve for where the derivative is zero to find the min
// and max distance along the line. We the have to solve those in terms of time
// rather than distance on the x-axis
void CubicBezier::InitRange(double p1y, double p2y) {
  range_min_ = 0;
  range_max_ = 1;
  if (0 <= p1y && p1y < 1 && 0 <= p2y && p2y <= 1)
    return;

  const double epsilon = kBezierEpsilon;

  // Represent the function's derivative in the form at^2 + bt + c
  // as in sampleCurveDerivativeY.
  // (Technically this is (dy/dt)*(1/3), which is suitable for finding zeros
  // but does not actually give the slope of the curve.)
  const double a = 3.0 * ay_;
  const double b = 2.0 * by_;
  const double c = cy_;

  // Check if the derivative is constant.
  if (std::abs(a) < epsilon && std::abs(b) < epsilon)
    return;

  // Zeros of the function's derivative.
  double t1 = 0;
  double t2 = 0;

  if (std::abs(a) < epsilon) {
    // The function's derivative is linear.
    t1 = -c / b;
  } else {
    // The function's derivative is a quadratic. We find the zeros of this
    // quadratic using the quadratic formula.
    double discriminant = b * b - 4 * a * c;
    if (discriminant < 0)
      return;
    double discriminant_sqrt = sqrt(discriminant);
    t1 = (-b + discriminant_sqrt) / (2 * a);
    t2 = (-b - discriminant_sqrt) / (2 * a);
  }

  double sol1 = 0;
  double sol2 = 0;

  // If the solution is in the range [0,1] then we include it, otherwise we
  // ignore it.

  // An interesting fact about these beziers is that they are only
  // actually evaluated in [0,1]. After that we take the tangent at that point
  // and linearly project it out.
  if (0 < t1 && t1 < 1)
    sol1 = SampleCurveY(t1);

  if (0 < t2 && t2 < 1)
    sol2 = SampleCurveY(t2);

  range_min_ = std::min(std::min(range_min_, sol1), sol2);
  range_max_ = std::max(std::max(range_max_, sol1), sol2);
}

double CubicBezier::GetDefaultEpsilon() {
  return kBezierEpsilon;
}

double CubicBezier::SolveCurveX(double x, double epsilon) const {
  DCHECK_GE(x, 0.0);
  DCHECK_LE(x, 1.0);

  double t0;
  double t1;
  double t2;
  double x2;
  double d2;
  int i;

  // First try a few iterations of Newton's method -- normally very fast.
  for (t2 = x, i = 0; i < 8; i++) {
    x2 = SampleCurveX(t2) - x;
    if (fabs(x2) < epsilon)
      return t2;
    d2 = SampleCurveDerivativeX(t2);
    if (fabs(d2) < 1e-6)
      break;
    t2 = t2 - x2 / d2;
  }

  // Fall back to the bisection method for reliability.
  t0 = 0.0;
  t1 = 1.0;
  t2 = x;

  while (t0 < t1) {
    x2 = SampleCurveX(t2);
    if (fabs(x2 - x) < epsilon)
      return t2;
    if (x > x2)
      t0 = t2;
    else
      t1 = t2;
    t2 = (t1 - t0) * .5 + t0;
  }

  // Failure.
  return t2;
}

double CubicBezier::Solve(double x) const {
  return SolveWithEpsilon(x, kBezierEpsilon);
}

double CubicBezier::SlopeWithEpsilon(double x, double epsilon) const {
  x = std::min(std::max(x, 0.0), 1.0);
  double t = SolveCurveX(x, epsilon);
  double dx = SampleCurveDerivativeX(t);
  double dy = SampleCurveDerivativeY(t);
  return dy / dx;
}

double CubicBezier::Slope(double x) const {
  return SlopeWithEpsilon(x, kBezierEpsilon);
}

double CubicBezier::GetX1() const {
  return cx_ / 3.0;
}

double CubicBezier::GetY1() const {
  return cy_ / 3.0;
}

double CubicBezier::GetX2() const {
  return (bx_ + cx_) / 3.0 + GetX1();
}

double CubicBezier::GetY2() const {
  return (by_ + cy_) / 3.0 + GetY1();
}

}  // namespace gfx