SkDCubicLineIntersection.cpp 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455
  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 "src/pathops/SkIntersections.h"
  8. #include "src/pathops/SkPathOpsCubic.h"
  9. #include "src/pathops/SkPathOpsCurve.h"
  10. #include "src/pathops/SkPathOpsLine.h"
  11. /*
  12. Find the interection of a line and cubic by solving for valid t values.
  13. Analogous to line-quadratic intersection, solve line-cubic intersection by
  14. representing the cubic as:
  15. x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3
  16. y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3
  17. and the line as:
  18. y = i*x + j (if the line is more horizontal)
  19. or:
  20. x = i*y + j (if the line is more vertical)
  21. Then using Mathematica, solve for the values of t where the cubic intersects the
  22. line:
  23. (in) Resultant[
  24. a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x,
  25. e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x]
  26. (out) -e + j +
  27. 3 e t - 3 f t -
  28. 3 e t^2 + 6 f t^2 - 3 g t^2 +
  29. e t^3 - 3 f t^3 + 3 g t^3 - h t^3 +
  30. i ( a -
  31. 3 a t + 3 b t +
  32. 3 a t^2 - 6 b t^2 + 3 c t^2 -
  33. a t^3 + 3 b t^3 - 3 c t^3 + d t^3 )
  34. if i goes to infinity, we can rewrite the line in terms of x. Mathematica:
  35. (in) Resultant[
  36. a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j,
  37. e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
  38. (out) a - j -
  39. 3 a t + 3 b t +
  40. 3 a t^2 - 6 b t^2 + 3 c t^2 -
  41. a t^3 + 3 b t^3 - 3 c t^3 + d t^3 -
  42. i ( e -
  43. 3 e t + 3 f t +
  44. 3 e t^2 - 6 f t^2 + 3 g t^2 -
  45. e t^3 + 3 f t^3 - 3 g t^3 + h t^3 )
  46. Solving this with Mathematica produces an expression with hundreds of terms;
  47. instead, use Numeric Solutions recipe to solve the cubic.
  48. The near-horizontal case, in terms of: Ax^3 + Bx^2 + Cx + D == 0
  49. A = (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d) )
  50. B = 3*(-( e - 2*f + g ) + i*( a - 2*b + c ) )
  51. C = 3*(-(-e + f ) + i*(-a + b ) )
  52. D = (-( e ) + i*( a ) + j )
  53. The near-vertical case, in terms of: Ax^3 + Bx^2 + Cx + D == 0
  54. A = ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h) )
  55. B = 3*( ( a - 2*b + c ) - i*( e - 2*f + g ) )
  56. C = 3*( (-a + b ) - i*(-e + f ) )
  57. D = ( ( a ) - i*( e ) - j )
  58. For horizontal lines:
  59. (in) Resultant[
  60. a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j,
  61. e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
  62. (out) e - j -
  63. 3 e t + 3 f t +
  64. 3 e t^2 - 6 f t^2 + 3 g t^2 -
  65. e t^3 + 3 f t^3 - 3 g t^3 + h t^3
  66. */
  67. class LineCubicIntersections {
  68. public:
  69. enum PinTPoint {
  70. kPointUninitialized,
  71. kPointInitialized
  72. };
  73. LineCubicIntersections(const SkDCubic& c, const SkDLine& l, SkIntersections* i)
  74. : fCubic(c)
  75. , fLine(l)
  76. , fIntersections(i)
  77. , fAllowNear(true) {
  78. i->setMax(4);
  79. }
  80. void allowNear(bool allow) {
  81. fAllowNear = allow;
  82. }
  83. void checkCoincident() {
  84. int last = fIntersections->used() - 1;
  85. for (int index = 0; index < last; ) {
  86. double cubicMidT = ((*fIntersections)[0][index] + (*fIntersections)[0][index + 1]) / 2;
  87. SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
  88. double t = fLine.nearPoint(cubicMidPt, nullptr);
  89. if (t < 0) {
  90. ++index;
  91. continue;
  92. }
  93. if (fIntersections->isCoincident(index)) {
  94. fIntersections->removeOne(index);
  95. --last;
  96. } else if (fIntersections->isCoincident(index + 1)) {
  97. fIntersections->removeOne(index + 1);
  98. --last;
  99. } else {
  100. fIntersections->setCoincident(index++);
  101. }
  102. fIntersections->setCoincident(index);
  103. }
  104. }
  105. // see parallel routine in line quadratic intersections
  106. int intersectRay(double roots[3]) {
  107. double adj = fLine[1].fX - fLine[0].fX;
  108. double opp = fLine[1].fY - fLine[0].fY;
  109. SkDCubic c;
  110. SkDEBUGCODE(c.fDebugGlobalState = fIntersections->globalState());
  111. for (int n = 0; n < 4; ++n) {
  112. c[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp;
  113. }
  114. double A, B, C, D;
  115. SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D);
  116. int count = SkDCubic::RootsValidT(A, B, C, D, roots);
  117. for (int index = 0; index < count; ++index) {
  118. SkDPoint calcPt = c.ptAtT(roots[index]);
  119. if (!approximately_zero(calcPt.fX)) {
  120. for (int n = 0; n < 4; ++n) {
  121. c[n].fY = (fCubic[n].fY - fLine[0].fY) * opp
  122. + (fCubic[n].fX - fLine[0].fX) * adj;
  123. }
  124. double extremeTs[6];
  125. int extrema = SkDCubic::FindExtrema(&c[0].fX, extremeTs);
  126. count = c.searchRoots(extremeTs, extrema, 0, SkDCubic::kXAxis, roots);
  127. break;
  128. }
  129. }
  130. return count;
  131. }
  132. int intersect() {
  133. addExactEndPoints();
  134. if (fAllowNear) {
  135. addNearEndPoints();
  136. }
  137. double rootVals[3];
  138. int roots = intersectRay(rootVals);
  139. for (int index = 0; index < roots; ++index) {
  140. double cubicT = rootVals[index];
  141. double lineT = findLineT(cubicT);
  142. SkDPoint pt;
  143. if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized) && uniqueAnswer(cubicT, pt)) {
  144. fIntersections->insert(cubicT, lineT, pt);
  145. }
  146. }
  147. checkCoincident();
  148. return fIntersections->used();
  149. }
  150. static int HorizontalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) {
  151. double A, B, C, D;
  152. SkDCubic::Coefficients(&c[0].fY, &A, &B, &C, &D);
  153. D -= axisIntercept;
  154. int count = SkDCubic::RootsValidT(A, B, C, D, roots);
  155. for (int index = 0; index < count; ++index) {
  156. SkDPoint calcPt = c.ptAtT(roots[index]);
  157. if (!approximately_equal(calcPt.fY, axisIntercept)) {
  158. double extremeTs[6];
  159. int extrema = SkDCubic::FindExtrema(&c[0].fY, extremeTs);
  160. count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kYAxis, roots);
  161. break;
  162. }
  163. }
  164. return count;
  165. }
  166. int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
  167. addExactHorizontalEndPoints(left, right, axisIntercept);
  168. if (fAllowNear) {
  169. addNearHorizontalEndPoints(left, right, axisIntercept);
  170. }
  171. double roots[3];
  172. int count = HorizontalIntersect(fCubic, axisIntercept, roots);
  173. for (int index = 0; index < count; ++index) {
  174. double cubicT = roots[index];
  175. SkDPoint pt = { fCubic.ptAtT(cubicT).fX, axisIntercept };
  176. double lineT = (pt.fX - left) / (right - left);
  177. if (pinTs(&cubicT, &lineT, &pt, kPointInitialized) && uniqueAnswer(cubicT, pt)) {
  178. fIntersections->insert(cubicT, lineT, pt);
  179. }
  180. }
  181. if (flipped) {
  182. fIntersections->flip();
  183. }
  184. checkCoincident();
  185. return fIntersections->used();
  186. }
  187. bool uniqueAnswer(double cubicT, const SkDPoint& pt) {
  188. for (int inner = 0; inner < fIntersections->used(); ++inner) {
  189. if (fIntersections->pt(inner) != pt) {
  190. continue;
  191. }
  192. double existingCubicT = (*fIntersections)[0][inner];
  193. if (cubicT == existingCubicT) {
  194. return false;
  195. }
  196. // check if midway on cubic is also same point. If so, discard this
  197. double cubicMidT = (existingCubicT + cubicT) / 2;
  198. SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
  199. if (cubicMidPt.approximatelyEqual(pt)) {
  200. return false;
  201. }
  202. }
  203. #if ONE_OFF_DEBUG
  204. SkDPoint cPt = fCubic.ptAtT(cubicT);
  205. SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
  206. cPt.fX, cPt.fY);
  207. #endif
  208. return true;
  209. }
  210. static int VerticalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) {
  211. double A, B, C, D;
  212. SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D);
  213. D -= axisIntercept;
  214. int count = SkDCubic::RootsValidT(A, B, C, D, roots);
  215. for (int index = 0; index < count; ++index) {
  216. SkDPoint calcPt = c.ptAtT(roots[index]);
  217. if (!approximately_equal(calcPt.fX, axisIntercept)) {
  218. double extremeTs[6];
  219. int extrema = SkDCubic::FindExtrema(&c[0].fX, extremeTs);
  220. count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kXAxis, roots);
  221. break;
  222. }
  223. }
  224. return count;
  225. }
  226. int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
  227. addExactVerticalEndPoints(top, bottom, axisIntercept);
  228. if (fAllowNear) {
  229. addNearVerticalEndPoints(top, bottom, axisIntercept);
  230. }
  231. double roots[3];
  232. int count = VerticalIntersect(fCubic, axisIntercept, roots);
  233. for (int index = 0; index < count; ++index) {
  234. double cubicT = roots[index];
  235. SkDPoint pt = { axisIntercept, fCubic.ptAtT(cubicT).fY };
  236. double lineT = (pt.fY - top) / (bottom - top);
  237. if (pinTs(&cubicT, &lineT, &pt, kPointInitialized) && uniqueAnswer(cubicT, pt)) {
  238. fIntersections->insert(cubicT, lineT, pt);
  239. }
  240. }
  241. if (flipped) {
  242. fIntersections->flip();
  243. }
  244. checkCoincident();
  245. return fIntersections->used();
  246. }
  247. protected:
  248. void addExactEndPoints() {
  249. for (int cIndex = 0; cIndex < 4; cIndex += 3) {
  250. double lineT = fLine.exactPoint(fCubic[cIndex]);
  251. if (lineT < 0) {
  252. continue;
  253. }
  254. double cubicT = (double) (cIndex >> 1);
  255. fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
  256. }
  257. }
  258. /* Note that this does not look for endpoints of the line that are near the cubic.
  259. These points are found later when check ends looks for missing points */
  260. void addNearEndPoints() {
  261. for (int cIndex = 0; cIndex < 4; cIndex += 3) {
  262. double cubicT = (double) (cIndex >> 1);
  263. if (fIntersections->hasT(cubicT)) {
  264. continue;
  265. }
  266. double lineT = fLine.nearPoint(fCubic[cIndex], nullptr);
  267. if (lineT < 0) {
  268. continue;
  269. }
  270. fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
  271. }
  272. this->addLineNearEndPoints();
  273. }
  274. void addLineNearEndPoints() {
  275. for (int lIndex = 0; lIndex < 2; ++lIndex) {
  276. double lineT = (double) lIndex;
  277. if (fIntersections->hasOppT(lineT)) {
  278. continue;
  279. }
  280. double cubicT = ((SkDCurve*) &fCubic)->nearPoint(SkPath::kCubic_Verb,
  281. fLine[lIndex], fLine[!lIndex]);
  282. if (cubicT < 0) {
  283. continue;
  284. }
  285. fIntersections->insert(cubicT, lineT, fLine[lIndex]);
  286. }
  287. }
  288. void addExactHorizontalEndPoints(double left, double right, double y) {
  289. for (int cIndex = 0; cIndex < 4; cIndex += 3) {
  290. double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y);
  291. if (lineT < 0) {
  292. continue;
  293. }
  294. double cubicT = (double) (cIndex >> 1);
  295. fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
  296. }
  297. }
  298. void addNearHorizontalEndPoints(double left, double right, double y) {
  299. for (int cIndex = 0; cIndex < 4; cIndex += 3) {
  300. double cubicT = (double) (cIndex >> 1);
  301. if (fIntersections->hasT(cubicT)) {
  302. continue;
  303. }
  304. double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y);
  305. if (lineT < 0) {
  306. continue;
  307. }
  308. fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
  309. }
  310. this->addLineNearEndPoints();
  311. }
  312. void addExactVerticalEndPoints(double top, double bottom, double x) {
  313. for (int cIndex = 0; cIndex < 4; cIndex += 3) {
  314. double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x);
  315. if (lineT < 0) {
  316. continue;
  317. }
  318. double cubicT = (double) (cIndex >> 1);
  319. fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
  320. }
  321. }
  322. void addNearVerticalEndPoints(double top, double bottom, double x) {
  323. for (int cIndex = 0; cIndex < 4; cIndex += 3) {
  324. double cubicT = (double) (cIndex >> 1);
  325. if (fIntersections->hasT(cubicT)) {
  326. continue;
  327. }
  328. double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x);
  329. if (lineT < 0) {
  330. continue;
  331. }
  332. fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
  333. }
  334. this->addLineNearEndPoints();
  335. }
  336. double findLineT(double t) {
  337. SkDPoint xy = fCubic.ptAtT(t);
  338. double dx = fLine[1].fX - fLine[0].fX;
  339. double dy = fLine[1].fY - fLine[0].fY;
  340. if (fabs(dx) > fabs(dy)) {
  341. return (xy.fX - fLine[0].fX) / dx;
  342. }
  343. return (xy.fY - fLine[0].fY) / dy;
  344. }
  345. bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
  346. if (!approximately_one_or_less(*lineT)) {
  347. return false;
  348. }
  349. if (!approximately_zero_or_more(*lineT)) {
  350. return false;
  351. }
  352. double cT = *cubicT = SkPinT(*cubicT);
  353. double lT = *lineT = SkPinT(*lineT);
  354. SkDPoint lPt = fLine.ptAtT(lT);
  355. SkDPoint cPt = fCubic.ptAtT(cT);
  356. if (!lPt.roughlyEqual(cPt)) {
  357. return false;
  358. }
  359. // FIXME: if points are roughly equal but not approximately equal, need to do
  360. // a binary search like quad/quad intersection to find more precise t values
  361. if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) {
  362. *pt = lPt;
  363. } else if (ptSet == kPointUninitialized) {
  364. *pt = cPt;
  365. }
  366. SkPoint gridPt = pt->asSkPoint();
  367. if (gridPt == fLine[0].asSkPoint()) {
  368. *lineT = 0;
  369. } else if (gridPt == fLine[1].asSkPoint()) {
  370. *lineT = 1;
  371. }
  372. if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) {
  373. *cubicT = 0;
  374. } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) {
  375. *cubicT = 1;
  376. }
  377. return true;
  378. }
  379. private:
  380. const SkDCubic& fCubic;
  381. const SkDLine& fLine;
  382. SkIntersections* fIntersections;
  383. bool fAllowNear;
  384. };
  385. int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y,
  386. bool flipped) {
  387. SkDLine line = {{{ left, y }, { right, y }}};
  388. LineCubicIntersections c(cubic, line, this);
  389. return c.horizontalIntersect(y, left, right, flipped);
  390. }
  391. int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x,
  392. bool flipped) {
  393. SkDLine line = {{{ x, top }, { x, bottom }}};
  394. LineCubicIntersections c(cubic, line, this);
  395. return c.verticalIntersect(x, top, bottom, flipped);
  396. }
  397. int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) {
  398. LineCubicIntersections c(cubic, line, this);
  399. c.allowNear(fAllowNear);
  400. return c.intersect();
  401. }
  402. int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) {
  403. LineCubicIntersections c(cubic, line, this);
  404. fUsed = c.intersectRay(fT[0]);
  405. for (int index = 0; index < fUsed; ++index) {
  406. fPt[index] = cubic.ptAtT(fT[0][index]);
  407. }
  408. return fUsed;
  409. }
  410. // SkDCubic accessors to Intersection utilities
  411. int SkDCubic::horizontalIntersect(double yIntercept, double roots[3]) const {
  412. return LineCubicIntersections::HorizontalIntersect(*this, yIntercept, roots);
  413. }
  414. int SkDCubic::verticalIntersect(double xIntercept, double roots[3]) const {
  415. return LineCubicIntersections::VerticalIntersect(*this, xIntercept, roots);
  416. }