CubicReduceOrder.cpp revision 47d73daa7a971e7eee5822def7922f7d43b2dc47
1a5728872c7702ddd09537c95bc3cbd20e1f2fb09Daniel Dunbar/* 2d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor * Copyright 2012 Google Inc. 3d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor * 4d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor * Use of this source code is governed by a BSD-style license that can be 5d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor * found in the LICENSE file. 6d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor */ 7d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor#include "CurveIntersection.h" 8d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor#include "Extrema.h" 9d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor#include "IntersectionUtilities.h" 10d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor#include "LineParameters.h" 11d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor 12d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregorstatic double interp_cubic_coords(const double* src, double t) 13d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor{ 14d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor double ab = interp(src[0], src[2], t); 15d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor double bc = interp(src[2], src[4], t); 16d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor double cd = interp(src[4], src[6], t); 17d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor double abc = interp(ab, bc, t); 18d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor double bcd = interp(bc, cd, t); 19d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor return interp(abc, bcd, t); 20d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor} 21d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregor 22d60e105e6d1624da647ef7dd35a9cf6fad1b763eDouglas Gregorstatic int coincident_line(const Cubic& cubic, Cubic& reduction) { 23d1102433214bd33b5bef5b493944292a1e82c2fbDouglas Gregor reduction[0] = reduction[1] = cubic[0]; 24d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor return 1; 25d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor} 26d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor 27d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregorstatic int vertical_line(const Cubic& cubic, ReduceOrder_Styles reduceStyle, Cubic& reduction) { 28d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor double tValues[2]; 29d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor reduction[0] = cubic[0]; 30d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor reduction[1] = cubic[3]; 31d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor if (reduceStyle == kReduceOrder_TreatAsFill) { 32d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor return 2; 33d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor } 34d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor int smaller = reduction[1].y > reduction[0].y; 35357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor int larger = smaller ^ 1; 36357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor int roots = findExtrema(cubic[0].y, cubic[1].y, cubic[2].y, cubic[3].y, tValues); 37357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor for (int index = 0; index < roots; ++index) { 38357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor double yExtrema = interp_cubic_coords(&cubic[0].y, tValues[index]); 39357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor if (reduction[smaller].y > yExtrema) { 40357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor reduction[smaller].y = yExtrema; 41357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor continue; 42357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor } 43d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor if (reduction[larger].y < yExtrema) { 445ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor reduction[larger].y = yExtrema; 455ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor } 465ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor } 475ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor return 2; 485ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor} 495ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor 505ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregorstatic int horizontal_line(const Cubic& cubic, ReduceOrder_Styles reduceStyle, Cubic& reduction) { 515ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor double tValues[2]; 525ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor reduction[0] = cubic[0]; 5318ef5e28a9a2677f8b1dce1fb2638d66e0a1621fDouglas Gregor reduction[1] = cubic[3]; 545ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor if (reduceStyle == kReduceOrder_TreatAsFill) { 555ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor return 2; 565ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor } 578320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor int smaller = reduction[1].x > reduction[0].x; 588320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor int larger = smaller ^ 1; 598320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor int roots = findExtrema(cubic[0].x, cubic[1].x, cubic[2].x, cubic[3].x, tValues); 608320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor for (int index = 0; index < roots; ++index) { 618320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor double xExtrema = interp_cubic_coords(&cubic[0].x, tValues[index]); 628320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor if (reduction[smaller].x > xExtrema) { 638320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor reduction[smaller].x = xExtrema; 648320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor continue; 658320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor } 663cd4d1ece34d36317ec5352855d86e256cb6aa27Douglas Gregor if (reduction[larger].x < xExtrema) { 673cd4d1ece34d36317ec5352855d86e256cb6aa27Douglas Gregor reduction[larger].x = xExtrema; 68d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor } 69d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor } 708320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor return 2; 718320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor} 723cd4d1ece34d36317ec5352855d86e256cb6aa27Douglas Gregor 733cd4d1ece34d36317ec5352855d86e256cb6aa27Douglas Gregor// check to see if it is a quadratic or a line 743cd4d1ece34d36317ec5352855d86e256cb6aa27Douglas Gregorstatic int check_quadratic(const Cubic& cubic, Cubic& reduction) { 753cd4d1ece34d36317ec5352855d86e256cb6aa27Douglas Gregor double dx10 = cubic[1].x - cubic[0].x; 768320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor double dx23 = cubic[2].x - cubic[3].x; 778320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor double midX = cubic[0].x + dx10 * 3 / 2; 7818ef5e28a9a2677f8b1dce1fb2638d66e0a1621fDouglas Gregor if (!AlmostEqualUlps(midX - cubic[3].x, dx23 * 3 / 2)) { 798320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor return 0; 808320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor } 815ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor double dy10 = cubic[1].y - cubic[0].y; 82d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor double dy23 = cubic[2].y - cubic[3].y; 83d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor double midY = cubic[0].y + dy10 * 3 / 2; 84d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor if (!AlmostEqualUlps(midY - cubic[3].y, dy23 * 3 / 2)) { 85d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor return 0; 86d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor } 87357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor reduction[0] = cubic[0]; 88357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor reduction[1].x = midX; 89357bbd022c1d340c8e255aea7a684ddb34bc76e5Douglas Gregor reduction[1].y = midY; 905ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor reduction[2] = cubic[3]; 915ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor return 3; 925ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor} 93aa0373107968aa7a26bf63f4a2673b8325b800afDouglas Gregor 945ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregorstatic int check_linear(const Cubic& cubic, ReduceOrder_Styles reduceStyle, 955ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor int minX, int maxX, int minY, int maxY, Cubic& reduction) { 965ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor int startIndex = 0; 975ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor int endIndex = 3; 985ec178ff237d77e7acf5e747c19bf4f2de77a779Douglas Gregor while (cubic[startIndex].approximatelyEqual(cubic[endIndex])) { 998320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor --endIndex; 1008320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor if (endIndex == 0) { 1018320ad400439def135a75aa1dd0cafbcc1c9e5e0Douglas Gregor printf("%s shouldn't get here if all four points are about equal\n", __FUNCTION__); 1023cd4d1ece34d36317ec5352855d86e256cb6aa27Douglas Gregor SkASSERT(0); 1033cd4d1ece34d36317ec5352855d86e256cb6aa27Douglas Gregor } 1043cd4d1ece34d36317ec5352855d86e256cb6aa27Douglas Gregor } 105d6350aefb1396b299e199c7f1fe51bb40c12e75eDouglas Gregor if (!isLinear(cubic, startIndex, endIndex)) { 10611a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor return 0; 10711a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor } 10811a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor // four are colinear: return line formed by outside 10911a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor reduction[0] = cubic[0]; 11011a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor reduction[1] = cubic[3]; 11111a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor if (reduceStyle == kReduceOrder_TreatAsFill) { 11211a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor return 2; 11311a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor } 11411a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor int sameSide1; 11511a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor int sameSide2; 11611a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor bool useX = cubic[maxX].x - cubic[minX].x >= cubic[maxY].y - cubic[minY].y; 11711a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor if (useX) { 11811a82401eaf3396dfb6466fa691a606204c8ddb0Douglas Gregor sameSide1 = sign(cubic[0].x - cubic[1].x) + sign(cubic[3].x - cubic[1].x); 119699ee52b2f383b62865013c3575510b520055811Douglas Gregor sameSide2 = sign(cubic[0].x - cubic[2].x) + sign(cubic[3].x - cubic[2].x); 120699ee52b2f383b62865013c3575510b520055811Douglas Gregor } else { 121699ee52b2f383b62865013c3575510b520055811Douglas Gregor sameSide1 = sign(cubic[0].y - cubic[1].y) + sign(cubic[3].y - cubic[1].y); 122699ee52b2f383b62865013c3575510b520055811Douglas Gregor sameSide2 = sign(cubic[0].y - cubic[2].y) + sign(cubic[3].y - cubic[2].y); 123699ee52b2f383b62865013c3575510b520055811Douglas Gregor } 124699ee52b2f383b62865013c3575510b520055811Douglas Gregor if (sameSide1 == sameSide2 && (sameSide1 & 3) != 2) { 125699ee52b2f383b62865013c3575510b520055811Douglas Gregor return 2; 126699ee52b2f383b62865013c3575510b520055811Douglas Gregor } 127699ee52b2f383b62865013c3575510b520055811Douglas Gregor double tValues[2]; 128699ee52b2f383b62865013c3575510b520055811Douglas Gregor int roots; 129699ee52b2f383b62865013c3575510b520055811Douglas Gregor if (useX) { 130699ee52b2f383b62865013c3575510b520055811Douglas Gregor roots = findExtrema(cubic[0].x, cubic[1].x, cubic[2].x, cubic[3].x, tValues); 131699ee52b2f383b62865013c3575510b520055811Douglas Gregor } else { 132699ee52b2f383b62865013c3575510b520055811Douglas Gregor roots = findExtrema(cubic[0].y, cubic[1].y, cubic[2].y, cubic[3].y, tValues); 133699ee52b2f383b62865013c3575510b520055811Douglas Gregor } 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 158bool 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: 172http://objectmix.com/graphics/132906-fast-precision-driven-cubic-quadratic-piecewise-degree-reduction-algos-2-a.html 173 174Given points c1, c2, c3 and c4 of a cubic Bezier, the points of the 175corresponding quadratic Bezier are (given in convex combinations of 176points): 177 178q1 = (11/13)c1 + (3/13)c2 -(3/13)c3 + (2/13)c4 179q2 = -c1 + (3/2)c2 + (3/2)c3 - c4 180q3 = (2/13)c1 - (3/13)c2 + (3/13)c3 + (11/13)c4 181 182Of course, this curve does not interpolate the end-points, but it would 183be interesting to see the behaviour of such a curve in an applet. 184 185-- 186Kalle Rutanen 187http://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 197int 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