SkGeometry.cpp | Explore in Territory

/*
 * Copyright 2006 The Android Open Source Project
 *
 * Use of this source code is governed by a BSD-style license that can be
 * found in the LICENSE file.
 */

#include "src/core/SkGeometry.h"

#include "include/core/SkMatrix.h"
#include "include/core/SkPoint3.h"
#include "include/core/SkRect.h"
#include "include/core/SkScalar.h"
#include "include/private/base/SkDebug.h"
#include "include/private/base/SkFloatingPoint.h"
#include "include/private/base/SkTPin.h"
#include "include/private/base/SkTo.h"
#include "src/base/SkBezierCurves.h"
#include "src/base/SkCubics.h"
#include "src/base/SkUtils.h"
#include "src/base/SkVx.h"
#include "src/core/SkPointPriv.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <cstdint>

namespace {

float2;
float4;

SkVector to_vector(const float2& x) { … }

////////////////////////////////////////////////////////////////////////

int is_not_monotonic(SkScalar a, SkScalar b, SkScalar c) { … }

////////////////////////////////////////////////////////////////////////

int valid_unit_divide(SkScalar numer, SkScalar denom, SkScalar* ratio) { … }

// Just returns its argument, but makes it easy to set a break-point to know when
// SkFindUnitQuadRoots is going to return 0 (an error).
int return_check_zero(int value) { … }

} // namespace

/** From Numerical Recipes in C.

    Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C])
    x1 = Q / A
    x2 = C / Q
*/
int SkFindUnitQuadRoots(SkScalar A, SkScalar B, SkScalar C, SkScalar roots[2]) { … }

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

void SkEvalQuadAt(const SkPoint src[3], SkScalar t, SkPoint* pt, SkVector* tangent) { … }

SkPoint SkEvalQuadAt(const SkPoint src[3], SkScalar t) { … }

SkVector SkEvalQuadTangentAt(const SkPoint src[3], SkScalar t) { … }

static inline float2 interp(const float2& v0,
                            const float2& v1,
                            const float2& t) { … }

void SkChopQuadAt(const SkPoint src[3], SkPoint dst[5], SkScalar t) { … }

void SkChopQuadAtHalf(const SkPoint src[3], SkPoint dst[5]) { … }

float SkMeasureAngleBetweenVectors(SkVector a, SkVector b) { … }

SkVector SkFindBisector(SkVector a, SkVector b) { … }

float SkFindQuadMidTangent(const SkPoint src[3]) { … }

/** Quad'(t) = At + B, where
    A = 2(a - 2b + c)
    B = 2(b - a)
    Solve for t, only if it fits between 0 < t < 1
*/
int SkFindQuadExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar tValue[1]) { … }

static inline void flatten_double_quad_extrema(SkScalar coords[14]) { … }

/*  Returns 0 for 1 quad, and 1 for two quads, either way the answer is
 stored in dst[]. Guarantees that the 1/2 quads will be monotonic.
 */
int SkChopQuadAtYExtrema(const SkPoint src[3], SkPoint dst[5]) { … }

/*  Returns 0 for 1 quad, and 1 for two quads, either way the answer is
    stored in dst[]. Guarantees that the 1/2 quads will be monotonic.
 */
int SkChopQuadAtXExtrema(const SkPoint src[3], SkPoint dst[5]) { … }

//  F(t)    = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2
//  F'(t)   = 2 (b - a) + 2 (a - 2b + c) t
//  F''(t)  = 2 (a - 2b + c)
//
//  A = 2 (b - a)
//  B = 2 (a - 2b + c)
//
//  Maximum curvature for a quadratic means solving
//  Fx' Fx'' + Fy' Fy'' = 0
//
//  t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2)
//
SkScalar SkFindQuadMaxCurvature(const SkPoint src[3]) { … }

int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5]) { … }

void SkConvertQuadToCubic(const SkPoint src[3], SkPoint dst[4]) { … }

//////////////////////////////////////////////////////////////////////////////
///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS /////
//////////////////////////////////////////////////////////////////////////////

static SkVector eval_cubic_derivative(const SkPoint src[4], SkScalar t) { … }

static SkVector eval_cubic_2ndDerivative(const SkPoint src[4], SkScalar t) { … }

void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc,
                   SkVector* tangent, SkVector* curvature) { … }

/** Cubic'(t) = At^2 + Bt + C, where
    A = 3(-a + 3(b - c) + d)
    B = 6(a - 2b + c)
    C = 3(b - a)
    Solve for t, keeping only those that fit betwee 0 < t < 1
*/
int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d,
                       SkScalar tValues[2]) { … }

// This does not return b when t==1, but it otherwise seems to get better precision than
// "a*(1 - t) + b*t" for things like chopping cubics on exact cusp points.
// The responsibility falls on the caller to check that t != 1 before calling.
template<int N, typename T>
inline static skvx::Vec<N,T> unchecked_mix(const skvx::Vec<N,T>& a, const skvx::Vec<N,T>& b,
                                           const skvx::Vec<N,T>& t) { … }

void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t) { … }

void SkChopCubicAt(const SkPoint src[4], SkPoint dst[10], float t0, float t1) { … }

void SkChopCubicAt(const SkPoint src[4], SkPoint dst[],
                   const SkScalar tValues[], int tCount) { … }

void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7]) { … }

float SkMeasureNonInflectCubicRotation(const SkPoint pts[4]) { … }

static skvx::float4 fma(const skvx::float4& f, float m, const skvx::float4& a) { … }

// Finds the root nearest 0.5. Returns 0.5 if the roots are undefined or outside 0..1.
static float solve_quadratic_equation_for_midtangent(float a, float b, float c, float discr) { … }

static float solve_quadratic_equation_for_midtangent(float a, float b, float c) { … }

float SkFindCubicMidTangent(const SkPoint src[4]) { … }

static void flatten_double_cubic_extrema(SkScalar coords[14]) { … }

/** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that
    the resulting beziers are monotonic in Y. This is called by the scan
    converter.  Depending on what is returned, dst[] is treated as follows:
    0   dst[0..3] is the original cubic
    1   dst[0..3] and dst[3..6] are the two new cubics
    2   dst[0..3], dst[3..6], dst[6..9] are the three new cubics
    If dst == null, it is ignored and only the count is returned.
*/
int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) { … }

int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) { … }

/** http://www.faculty.idc.ac.il/arik/quality/appendixA.html

    Inflection means that curvature is zero.
    Curvature is [F' x F''] / [F'^3]
    So we solve F'x X F''y - F'y X F''y == 0
    After some canceling of the cubic term, we get
    A = b - a
    B = c - 2b + a
    C = d - 3c + 3b - a
    (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0
*/
int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[2]) { … }

int SkChopCubicAtInflections(const SkPoint src[4], SkPoint dst[10]) { … }

// Assumes the third component of points is 1.
// Calcs p0 . (p1 x p2)
static double calc_dot_cross_cubic(const SkPoint& p0, const SkPoint& p1, const SkPoint& p2) { … }

// Returns a positive power of 2 that, when multiplied by n, and excepting the two edge cases listed
// below, shifts the exponent of n to yield a magnitude somewhere inside [1..2).
// Returns 2^1023 if abs(n) < 2^-1022 (including 0).
// Returns NaN if n is Inf or NaN.
inline static double previous_inverse_pow2(double n) { … }

inline static void write_cubic_inflection_roots(double t0, double s0, double t1, double s1,
                                                double* t, double* s) { … }

SkCubicType SkClassifyCubic(const SkPoint P[4], double t[2], double s[2], double d[4]) { … }

template <typename T> void bubble_sort(T array[], int count) { … }

/**
 *  Given an array and count, remove all pair-wise duplicates from the array,
 *  keeping the existing sorting, and return the new count
 */
static int collaps_duplicates(SkScalar array[], int count) { … }

#ifdef SK_DEBUG

#define TEST_COLLAPS_ENTRY(array) …

static void test_collaps_duplicates() { … }
#endif

static SkScalar SkScalarCubeRoot(SkScalar x) { … }

/*  Solve coeff(t) == 0, returning the number of roots that
    lie withing 0 < t < 1.
    coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3]

    Eliminates repeated roots (so that all tValues are distinct, and are always
    in increasing order.
*/
static int solve_cubic_poly(const SkScalar coeff[4], SkScalar tValues[3]) { … }

/*  Looking for F' dot F'' == 0

    A = b - a
    B = c - 2b + a
    C = d - 3c + 3b - a

    F' = 3Ct^2 + 6Bt + 3A
    F'' = 6Ct + 6B

    F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
*/
static void formulate_F1DotF2(const SkScalar src[], SkScalar coeff[4]) { … }

/*  Looking for F' dot F'' == 0

    A = b - a
    B = c - 2b + a
    C = d - 3c + 3b - a

    F' = 3Ct^2 + 6Bt + 3A
    F'' = 6Ct + 6B

    F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
*/
int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3]) { … }

int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13],
                              SkScalar tValues[3]) { … }

// Returns a constant proportional to the dimensions of the cubic.
// Constant found through experimentation -- maybe there's a better way....
static SkScalar calc_cubic_precision(const SkPoint src[4]) { … }

// Returns true if both points src[testIndex], src[testIndex+1] are in the same half plane defined
// by the line segment src[lineIndex], src[lineIndex+1].
static bool on_same_side(const SkPoint src[4], int testIndex, int lineIndex) { … }

// Return location (in t) of cubic cusp, if there is one.
// Note that classify cubic code does not reliably return all cusp'd cubics, so
// it is not called here.
SkScalar SkFindCubicCusp(const SkPoint src[4]) { … }

static bool close_enough_to_zero(double x) { … }

static bool first_axis_intersection(const double coefficients[8], bool yDirection,
                                    double axisIntercept, double* solution) { … }

bool SkChopMonoCubicAtY(const SkPoint src[4], SkScalar y, SkPoint dst[7]) { … }

bool SkChopMonoCubicAtX(const SkPoint src[4], SkScalar x, SkPoint dst[7]) { … }

///////////////////////////////////////////////////////////////////////////////
//
// NURB representation for conics.  Helpful explanations at:
//
// http://citeseerx.ist.psu.edu/viewdoc/
//   download?doi=10.1.1.44.5740&rep=rep1&type=ps
// and
// http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/NURBS/RB-conics.html
//
// F = (A (1 - t)^2 + C t^2 + 2 B (1 - t) t w)
//     ------------------------------------------
//         ((1 - t)^2 + t^2 + 2 (1 - t) t w)
//
//   = {t^2 (P0 + P2 - 2 P1 w), t (-2 P0 + 2 P1 w), P0}
//     ------------------------------------------------
//             {t^2 (2 - 2 w), t (-2 + 2 w), 1}
//

// F' = 2 (C t (1 + t (-1 + w)) - A (-1 + t) (t (-1 + w) - w) + B (1 - 2 t) w)
//
//  t^2 : (2 P0 - 2 P2 - 2 P0 w + 2 P2 w)
//  t^1 : (-2 P0 + 2 P2 + 4 P0 w - 4 P1 w)
//  t^0 : -2 P0 w + 2 P1 w
//
//  We disregard magnitude, so we can freely ignore the denominator of F', and
//  divide the numerator by 2
//
//    coeff[0] for t^2
//    coeff[1] for t^1
//    coeff[2] for t^0
//
static void conic_deriv_coeff(const SkScalar src[],
                              SkScalar w,
                              SkScalar coeff[3]) { … }

static bool conic_find_extrema(const SkScalar src[], SkScalar w, SkScalar* t) { … }

// We only interpolate one dimension at a time (the first, at +0, +3, +6).
static void p3d_interp(const SkScalar src[7], SkScalar dst[7], SkScalar t) { … }

static void ratquad_mapTo3D(const SkPoint src[3], SkScalar w, SkPoint3 dst[3]) { … }

static SkPoint project_down(const SkPoint3& src) { … }

// return false if infinity or NaN is generated; caller must check
bool SkConic::chopAt(SkScalar t, SkConic dst[2]) const { … }

void SkConic::chopAt(SkScalar t1, SkScalar t2, SkConic* dst) const { … }

SkPoint SkConic::evalAt(SkScalar t) const { … }

SkVector SkConic::evalTangentAt(SkScalar t) const { … }

void SkConic::evalAt(SkScalar t, SkPoint* pt, SkVector* tangent) const { … }

static SkScalar subdivide_w_value(SkScalar w) { … }

#if defined(SK_SUPPORT_LEGACY_CONIC_CHOP)
void SkConic::chop(SkConic * SK_RESTRICT dst) const { … }
#else
void SkConic::chop(SkConic * SK_RESTRICT dst) const {

    // Observe that scale will always be smaller than 1 because fW > 0.
    const float scale = SkScalarInvert(SK_Scalar1 + fW);

    // The subdivided control points below are the sums of the following three terms. Because the
    // terms are multiplied by something <1, and the resulting control points lie within the
    // control points of the original then the terms and the sums below will not overflow. Note
    // that fW * scale approaches 1 as fW becomes very large.
    float2 t0 = from_point(fPts[0]) * scale;
    float2 t1 = from_point(fPts[1]) * (fW * scale);
    float2 t2 = from_point(fPts[2]) * scale;

    // Calculate the subdivided control points
    const SkPoint p1 = to_point(t0 + t1);
    const SkPoint p3 = to_point(t1 + t2);

    // p2 = (t0 + 2*t1 + t2) / 2. Divide the terms by 2 before the sum to keep the sum for p2
    // from overflowing.
    const SkPoint p2 = to_point(0.5f * t0 + t1 + 0.5f * t2);

    SkASSERT(p1.isFinite() && p2.isFinite() && p3.isFinite());

    dst[0].fPts[0] = fPts[0];
    dst[0].fPts[1] = p1;
    dst[0].fPts[2] = p2;
    dst[1].fPts[0] = p2;
    dst[1].fPts[1] = p3;
    dst[1].fPts[2] = fPts[2];

    // Update w.
    dst[0].fW = dst[1].fW = subdivide_w_value(fW);
}
#endif  // SK_SUPPORT_LEGACY_CONIC_CHOP

/*
 *  "High order approximation of conic sections by quadratic splines"
 *      by Michael Floater, 1993
 */
#define AS_QUAD_ERROR_SETUP …

void SkConic::computeAsQuadError(SkVector* err) const { … }

bool SkConic::asQuadTol(SkScalar tol) const { … }

// Limit the number of suggested quads to approximate a conic
#define kMaxConicToQuadPOW2 …

int SkConic::computeQuadPOW2(SkScalar tol) const { … }

// This was originally developed and tested for pathops: see SkOpTypes.h
// returns true if (a <= b <= c) || (a >= b >= c)
static bool between(SkScalar a, SkScalar b, SkScalar c) { … }

static SkPoint* subdivide(const SkConic& src, SkPoint pts[], int level) { … }

int SkConic::chopIntoQuadsPOW2(SkPoint pts[], int pow2) const { … }

float SkConic::findMidTangent() const { … }

bool SkConic::findXExtrema(SkScalar* t) const { … }

bool SkConic::findYExtrema(SkScalar* t) const { … }

bool SkConic::chopAtXExtrema(SkConic dst[2]) const { … }

bool SkConic::chopAtYExtrema(SkConic dst[2]) const { … }

void SkConic::computeTightBounds(SkRect* bounds) const { … }

void SkConic::computeFastBounds(SkRect* bounds) const { … }

#if 0  // unimplemented
bool SkConic::findMaxCurvature(SkScalar* t) const {
    // TODO: Implement me
    return false;
}
#endif

SkScalar SkConic::TransformW(const SkPoint pts[3], SkScalar w, const SkMatrix& matrix) { … }

int SkConic::BuildUnitArc(const SkVector& uStart, const SkVector& uStop, SkRotationDirection dir,
                          const SkMatrix* userMatrix, SkConic dst[kMaxConicsForArc]) { … }
chromium/third_party/skia/src/core/SkGeometry.cpp
