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 "CurveIntersection.h"
      8 #include "Extrema.h"
      9 #include "IntersectionUtilities.h"
     10 #include "LineParameters.h"
     11 
     12 static double interp_cubic_coords(const double* src, double t)
     13 {
     14     double ab = interp(src[0], src[2], t);
     15     double bc = interp(src[2], src[4], t);
     16     double cd = interp(src[4], src[6], t);
     17     double abc = interp(ab, bc, t);
     18     double bcd = interp(bc, cd, t);
     19     return interp(abc, bcd, t);
     20 }
     21 
     22 static int coincident_line(const Cubic& cubic, Cubic& reduction) {
     23     reduction[0] = reduction[1] = cubic[0];
     24     return 1;
     25 }
     26 
     27 static int vertical_line(const Cubic& cubic, ReduceOrder_Styles reduceStyle, Cubic& reduction) {
     28     double tValues[2];
     29     reduction[0] = cubic[0];
     30     reduction[1] = cubic[3];
     31     if (reduceStyle == kReduceOrder_TreatAsFill) {
     32         return 2;
     33     }
     34     int smaller = reduction[1].y > reduction[0].y;
     35     int larger = smaller ^ 1;
     36     int roots = findExtrema(cubic[0].y, cubic[1].y, cubic[2].y, cubic[3].y, tValues);
     37     for (int index = 0; index < roots; ++index) {
     38         double yExtrema = interp_cubic_coords(&cubic[0].y, tValues[index]);
     39         if (reduction[smaller].y > yExtrema) {
     40             reduction[smaller].y = yExtrema;
     41             continue;
     42         }
     43         if (reduction[larger].y < yExtrema) {
     44             reduction[larger].y = yExtrema;
     45         }
     46     }
     47     return 2;
     48 }
     49 
     50 static int horizontal_line(const Cubic& cubic, ReduceOrder_Styles reduceStyle, Cubic& reduction) {
     51     double tValues[2];
     52     reduction[0] = cubic[0];
     53     reduction[1] = cubic[3];
     54     if (reduceStyle == kReduceOrder_TreatAsFill) {
     55         return 2;
     56     }
     57     int smaller = reduction[1].x > reduction[0].x;
     58     int larger = smaller ^ 1;
     59     int roots = findExtrema(cubic[0].x, cubic[1].x, cubic[2].x, cubic[3].x, tValues);
     60     for (int index = 0; index < roots; ++index) {
     61         double xExtrema = interp_cubic_coords(&cubic[0].x, tValues[index]);
     62         if (reduction[smaller].x > xExtrema) {
     63             reduction[smaller].x = xExtrema;
     64             continue;
     65         }
     66         if (reduction[larger].x < xExtrema) {
     67             reduction[larger].x = xExtrema;
     68         }
     69     }
     70     return 2;
     71 }
     72 
     73 // check to see if it is a quadratic or a line
     74 static int check_quadratic(const Cubic& cubic, Cubic& reduction) {
     75     double dx10 = cubic[1].x - cubic[0].x;
     76     double dx23 = cubic[2].x - cubic[3].x;
     77     double midX = cubic[0].x + dx10 * 3 / 2;
     78     if (!AlmostEqualUlps(midX - cubic[3].x, dx23 * 3 / 2)) {
     79         return 0;
     80     }
     81     double dy10 = cubic[1].y - cubic[0].y;
     82     double dy23 = cubic[2].y - cubic[3].y;
     83     double midY = cubic[0].y + dy10 * 3 / 2;
     84     if (!AlmostEqualUlps(midY - cubic[3].y, dy23 * 3 / 2)) {
     85         return 0;
     86     }
     87     reduction[0] = cubic[0];
     88     reduction[1].x = midX;
     89     reduction[1].y = midY;
     90     reduction[2] = cubic[3];
     91     return 3;
     92 }
     93 
     94 static int check_linear(const Cubic& cubic, ReduceOrder_Styles reduceStyle,
     95         int minX, int maxX, int minY, int maxY, Cubic& reduction) {
     96     int startIndex = 0;
     97     int endIndex = 3;
     98     while (cubic[startIndex].approximatelyEqual(cubic[endIndex])) {
     99         --endIndex;
    100         if (endIndex == 0) {
    101             printf("%s shouldn't get here if all four points are about equal\n", __FUNCTION__);
    102             SkASSERT(0);
    103         }
    104     }
    105     if (!isLinear(cubic, startIndex, endIndex)) {
    106         return 0;
    107     }
    108     // four are colinear: return line formed by outside
    109     reduction[0] = cubic[0];
    110     reduction[1] = cubic[3];
    111     if (reduceStyle == kReduceOrder_TreatAsFill) {
    112         return 2;
    113     }
    114     int sameSide1;
    115     int sameSide2;
    116     bool useX = cubic[maxX].x - cubic[minX].x >= cubic[maxY].y - cubic[minY].y;
    117     if (useX) {
    118         sameSide1 = sign(cubic[0].x - cubic[1].x) + sign(cubic[3].x - cubic[1].x);
    119         sameSide2 = sign(cubic[0].x - cubic[2].x) + sign(cubic[3].x - cubic[2].x);
    120     } else {
    121         sameSide1 = sign(cubic[0].y - cubic[1].y) + sign(cubic[3].y - cubic[1].y);
    122         sameSide2 = sign(cubic[0].y - cubic[2].y) + sign(cubic[3].y - cubic[2].y);
    123     }
    124     if (sameSide1 == sameSide2 && (sameSide1 & 3) != 2) {
    125         return 2;
    126     }
    127     double tValues[2];
    128     int roots;
    129     if (useX) {
    130         roots = findExtrema(cubic[0].x, cubic[1].x, cubic[2].x, cubic[3].x, tValues);
    131     } else {
    132         roots = findExtrema(cubic[0].y, cubic[1].y, cubic[2].y, cubic[3].y, tValues);
    133     }
    134     for (int index = 0; index < roots; ++index) {
    135         _Point extrema;
    136         extrema.x = interp_cubic_coords(&cubic[0].x, tValues[index]);
    137         extrema.y = interp_cubic_coords(&cubic[0].y, tValues[index]);
    138         // sameSide > 0 means mid is smaller than either [0] or [3], so replace smaller
    139         int replace;
    140         if (useX) {
    141             if (extrema.x < cubic[0].x ^ extrema.x < cubic[3].x) {
    142                 continue;
    143             }
    144             replace = (extrema.x < cubic[0].x | extrema.x < cubic[3].x)
    145                     ^ (cubic[0].x < cubic[3].x);
    146         } else {
    147             if (extrema.y < cubic[0].y ^ extrema.y < cubic[3].y) {
    148                 continue;
    149             }
    150             replace = (extrema.y < cubic[0].y | extrema.y < cubic[3].y)
    151                     ^ (cubic[0].y < cubic[3].y);
    152         }
    153         reduction[replace] = extrema;
    154     }
    155     return 2;
    156 }
    157 
    158 bool isLinear(const Cubic& cubic, int startIndex, int endIndex) {
    159     LineParameters lineParameters;
    160     lineParameters.cubicEndPoints(cubic, startIndex, endIndex);
    161     // FIXME: maybe it's possible to avoid this and compare non-normalized
    162     lineParameters.normalize();
    163     double distance = lineParameters.controlPtDistance(cubic, 1);
    164     if (!approximately_zero(distance)) {
    165         return false;
    166     }
    167     distance = lineParameters.controlPtDistance(cubic, 2);
    168     return approximately_zero(distance);
    169 }
    170 
    171 /* food for thought:
    172 http://objectmix.com/graphics/132906-fast-precision-driven-cubic-quadratic-piecewise-degree-reduction-algos-2-a.html
    173 
    174 Given points c1, c2, c3 and c4 of a cubic Bezier, the points of the
    175 corresponding quadratic Bezier are (given in convex combinations of
    176 points):
    177 
    178 q1 = (11/13)c1 + (3/13)c2 -(3/13)c3 + (2/13)c4
    179 q2 = -c1 + (3/2)c2 + (3/2)c3 - c4
    180 q3 = (2/13)c1 - (3/13)c2 + (3/13)c3 + (11/13)c4
    181 
    182 Of course, this curve does not interpolate the end-points, but it would
    183 be interesting to see the behaviour of such a curve in an applet.
    184 
    185 --
    186 Kalle Rutanen
    187 http://kaba.hilvi.org
    188 
    189 */
    190 
    191 // reduce to a quadratic or smaller
    192 // look for identical points
    193 // look for all four points in a line
    194     // note that three points in a line doesn't simplify a cubic
    195 // look for approximation with single quadratic
    196     // save approximation with multiple quadratics for later
    197 int reduceOrder(const Cubic& cubic, Cubic& reduction, ReduceOrder_Quadratics allowQuadratics,
    198         ReduceOrder_Styles reduceStyle) {
    199     int index, minX, maxX, minY, maxY;
    200     int minXSet, minYSet;
    201     minX = maxX = minY = maxY = 0;
    202     minXSet = minYSet = 0;
    203     for (index = 1; index < 4; ++index) {
    204         if (cubic[minX].x > cubic[index].x) {
    205             minX = index;
    206         }
    207         if (cubic[minY].y > cubic[index].y) {
    208             minY = index;
    209         }
    210         if (cubic[maxX].x < cubic[index].x) {
    211             maxX = index;
    212         }
    213         if (cubic[maxY].y < cubic[index].y) {
    214             maxY = index;
    215         }
    216     }
    217     for (index = 0; index < 4; ++index) {
    218         double cx = cubic[index].x;
    219         double cy = cubic[index].y;
    220         double denom = SkTMax(fabs(cx), SkTMax(fabs(cy),
    221                 SkTMax(fabs(cubic[minX].x), fabs(cubic[minY].y))));
    222         if (denom == 0) {
    223             minXSet |= 1 << index;
    224             minYSet |= 1 << index;
    225             continue;
    226         }
    227         double inv = 1 / denom;
    228         if (approximately_equal_half(cx * inv, cubic[minX].x * inv)) {
    229             minXSet |= 1 << index;
    230         }
    231         if (approximately_equal_half(cy * inv, cubic[minY].y * inv)) {
    232             minYSet |= 1 << index;
    233         }
    234     }
    235     if (minXSet == 0xF) { // test for vertical line
    236         if (minYSet == 0xF) { // return 1 if all four are coincident
    237             return coincident_line(cubic, reduction);
    238         }
    239         return vertical_line(cubic, reduceStyle, reduction);
    240     }
    241     if (minYSet == 0xF) { // test for horizontal line
    242         return horizontal_line(cubic, reduceStyle, reduction);
    243     }
    244     int result = check_linear(cubic, reduceStyle, minX, maxX, minY, maxY, reduction);
    245     if (result) {
    246         return result;
    247     }
    248     if (allowQuadratics == kReduceOrder_QuadraticsAllowed
    249             && (result = check_quadratic(cubic, reduction))) {
    250         return result;
    251     }
    252     memcpy(reduction, cubic, sizeof(Cubic));
    253     return 4;
    254 }
    255