SkDQuadLineIntersection.cpp 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470
  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/SkPathOpsCurve.h"
  9. #include "src/pathops/SkPathOpsLine.h"
  10. #include "src/pathops/SkPathOpsQuad.h"
  11. /*
  12. Find the interection of a line and quadratic by solving for valid t values.
  13. From http://stackoverflow.com/questions/1853637/how-to-find-the-mathematical-function-defining-a-bezier-curve
  14. "A Bezier curve is a parametric function. A quadratic Bezier curve (i.e. three
  15. control points) can be expressed as: F(t) = A(1 - t)^2 + B(1 - t)t + Ct^2 where
  16. A, B and C are points and t goes from zero to one.
  17. This will give you two equations:
  18. x = a(1 - t)^2 + b(1 - t)t + ct^2
  19. y = d(1 - t)^2 + e(1 - t)t + ft^2
  20. If you add for instance the line equation (y = kx + m) to that, you'll end up
  21. with three equations and three unknowns (x, y and t)."
  22. Similar to above, the quadratic is represented as
  23. x = a(1-t)^2 + 2b(1-t)t + ct^2
  24. y = d(1-t)^2 + 2e(1-t)t + ft^2
  25. and the line as
  26. y = g*x + h
  27. Using Mathematica, solve for the values of t where the quadratic intersects the
  28. line:
  29. (in) t1 = Resultant[a*(1 - t)^2 + 2*b*(1 - t)*t + c*t^2 - x,
  30. d*(1 - t)^2 + 2*e*(1 - t)*t + f*t^2 - g*x - h, x]
  31. (out) -d + h + 2 d t - 2 e t - d t^2 + 2 e t^2 - f t^2 +
  32. g (a - 2 a t + 2 b t + a t^2 - 2 b t^2 + c t^2)
  33. (in) Solve[t1 == 0, t]
  34. (out) {
  35. {t -> (-2 d + 2 e + 2 a g - 2 b g -
  36. Sqrt[(2 d - 2 e - 2 a g + 2 b g)^2 -
  37. 4 (-d + 2 e - f + a g - 2 b g + c g) (-d + a g + h)]) /
  38. (2 (-d + 2 e - f + a g - 2 b g + c g))
  39. },
  40. {t -> (-2 d + 2 e + 2 a g - 2 b g +
  41. Sqrt[(2 d - 2 e - 2 a g + 2 b g)^2 -
  42. 4 (-d + 2 e - f + a g - 2 b g + c g) (-d + a g + h)]) /
  43. (2 (-d + 2 e - f + a g - 2 b g + c g))
  44. }
  45. }
  46. Using the results above (when the line tends towards horizontal)
  47. A = (-(d - 2*e + f) + g*(a - 2*b + c) )
  48. B = 2*( (d - e ) - g*(a - b ) )
  49. C = (-(d ) + g*(a ) + h )
  50. If g goes to infinity, we can rewrite the line in terms of x.
  51. x = g'*y + h'
  52. And solve accordingly in Mathematica:
  53. (in) t2 = Resultant[a*(1 - t)^2 + 2*b*(1 - t)*t + c*t^2 - g'*y - h',
  54. d*(1 - t)^2 + 2*e*(1 - t)*t + f*t^2 - y, y]
  55. (out) a - h' - 2 a t + 2 b t + a t^2 - 2 b t^2 + c t^2 -
  56. g' (d - 2 d t + 2 e t + d t^2 - 2 e t^2 + f t^2)
  57. (in) Solve[t2 == 0, t]
  58. (out) {
  59. {t -> (2 a - 2 b - 2 d g' + 2 e g' -
  60. Sqrt[(-2 a + 2 b + 2 d g' - 2 e g')^2 -
  61. 4 (a - 2 b + c - d g' + 2 e g' - f g') (a - d g' - h')]) /
  62. (2 (a - 2 b + c - d g' + 2 e g' - f g'))
  63. },
  64. {t -> (2 a - 2 b - 2 d g' + 2 e g' +
  65. Sqrt[(-2 a + 2 b + 2 d g' - 2 e g')^2 -
  66. 4 (a - 2 b + c - d g' + 2 e g' - f g') (a - d g' - h')])/
  67. (2 (a - 2 b + c - d g' + 2 e g' - f g'))
  68. }
  69. }
  70. Thus, if the slope of the line tends towards vertical, we use:
  71. A = ( (a - 2*b + c) - g'*(d - 2*e + f) )
  72. B = 2*(-(a - b ) + g'*(d - e ) )
  73. C = ( (a ) - g'*(d ) - h' )
  74. */
  75. class LineQuadraticIntersections {
  76. public:
  77. enum PinTPoint {
  78. kPointUninitialized,
  79. kPointInitialized
  80. };
  81. LineQuadraticIntersections(const SkDQuad& q, const SkDLine& l, SkIntersections* i)
  82. : fQuad(q)
  83. , fLine(&l)
  84. , fIntersections(i)
  85. , fAllowNear(true) {
  86. i->setMax(5); // allow short partial coincidence plus discrete intersections
  87. }
  88. LineQuadraticIntersections(const SkDQuad& q)
  89. : fQuad(q)
  90. SkDEBUGPARAMS(fLine(nullptr))
  91. SkDEBUGPARAMS(fIntersections(nullptr))
  92. SkDEBUGPARAMS(fAllowNear(false)) {
  93. }
  94. void allowNear(bool allow) {
  95. fAllowNear = allow;
  96. }
  97. void checkCoincident() {
  98. int last = fIntersections->used() - 1;
  99. for (int index = 0; index < last; ) {
  100. double quadMidT = ((*fIntersections)[0][index] + (*fIntersections)[0][index + 1]) / 2;
  101. SkDPoint quadMidPt = fQuad.ptAtT(quadMidT);
  102. double t = fLine->nearPoint(quadMidPt, nullptr);
  103. if (t < 0) {
  104. ++index;
  105. continue;
  106. }
  107. if (fIntersections->isCoincident(index)) {
  108. fIntersections->removeOne(index);
  109. --last;
  110. } else if (fIntersections->isCoincident(index + 1)) {
  111. fIntersections->removeOne(index + 1);
  112. --last;
  113. } else {
  114. fIntersections->setCoincident(index++);
  115. }
  116. fIntersections->setCoincident(index);
  117. }
  118. }
  119. int intersectRay(double roots[2]) {
  120. /*
  121. solve by rotating line+quad so line is horizontal, then finding the roots
  122. set up matrix to rotate quad to x-axis
  123. |cos(a) -sin(a)|
  124. |sin(a) cos(a)|
  125. note that cos(a) = A(djacent) / Hypoteneuse
  126. sin(a) = O(pposite) / Hypoteneuse
  127. since we are computing Ts, we can ignore hypoteneuse, the scale factor:
  128. | A -O |
  129. | O A |
  130. A = line[1].fX - line[0].fX (adjacent side of the right triangle)
  131. O = line[1].fY - line[0].fY (opposite side of the right triangle)
  132. for each of the three points (e.g. n = 0 to 2)
  133. quad[n].fY' = (quad[n].fY - line[0].fY) * A - (quad[n].fX - line[0].fX) * O
  134. */
  135. double adj = (*fLine)[1].fX - (*fLine)[0].fX;
  136. double opp = (*fLine)[1].fY - (*fLine)[0].fY;
  137. double r[3];
  138. for (int n = 0; n < 3; ++n) {
  139. r[n] = (fQuad[n].fY - (*fLine)[0].fY) * adj - (fQuad[n].fX - (*fLine)[0].fX) * opp;
  140. }
  141. double A = r[2];
  142. double B = r[1];
  143. double C = r[0];
  144. A += C - 2 * B; // A = a - 2*b + c
  145. B -= C; // B = -(b - c)
  146. return SkDQuad::RootsValidT(A, 2 * B, C, roots);
  147. }
  148. int intersect() {
  149. addExactEndPoints();
  150. if (fAllowNear) {
  151. addNearEndPoints();
  152. }
  153. double rootVals[2];
  154. int roots = intersectRay(rootVals);
  155. for (int index = 0; index < roots; ++index) {
  156. double quadT = rootVals[index];
  157. double lineT = findLineT(quadT);
  158. SkDPoint pt;
  159. if (pinTs(&quadT, &lineT, &pt, kPointUninitialized) && uniqueAnswer(quadT, pt)) {
  160. fIntersections->insert(quadT, lineT, pt);
  161. }
  162. }
  163. checkCoincident();
  164. return fIntersections->used();
  165. }
  166. int horizontalIntersect(double axisIntercept, double roots[2]) {
  167. double D = fQuad[2].fY; // f
  168. double E = fQuad[1].fY; // e
  169. double F = fQuad[0].fY; // d
  170. D += F - 2 * E; // D = d - 2*e + f
  171. E -= F; // E = -(d - e)
  172. F -= axisIntercept;
  173. return SkDQuad::RootsValidT(D, 2 * E, F, roots);
  174. }
  175. int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
  176. addExactHorizontalEndPoints(left, right, axisIntercept);
  177. if (fAllowNear) {
  178. addNearHorizontalEndPoints(left, right, axisIntercept);
  179. }
  180. double rootVals[2];
  181. int roots = horizontalIntersect(axisIntercept, rootVals);
  182. for (int index = 0; index < roots; ++index) {
  183. double quadT = rootVals[index];
  184. SkDPoint pt = fQuad.ptAtT(quadT);
  185. double lineT = (pt.fX - left) / (right - left);
  186. if (pinTs(&quadT, &lineT, &pt, kPointInitialized) && uniqueAnswer(quadT, pt)) {
  187. fIntersections->insert(quadT, lineT, pt);
  188. }
  189. }
  190. if (flipped) {
  191. fIntersections->flip();
  192. }
  193. checkCoincident();
  194. return fIntersections->used();
  195. }
  196. bool uniqueAnswer(double quadT, const SkDPoint& pt) {
  197. for (int inner = 0; inner < fIntersections->used(); ++inner) {
  198. if (fIntersections->pt(inner) != pt) {
  199. continue;
  200. }
  201. double existingQuadT = (*fIntersections)[0][inner];
  202. if (quadT == existingQuadT) {
  203. return false;
  204. }
  205. // check if midway on quad is also same point. If so, discard this
  206. double quadMidT = (existingQuadT + quadT) / 2;
  207. SkDPoint quadMidPt = fQuad.ptAtT(quadMidT);
  208. if (quadMidPt.approximatelyEqual(pt)) {
  209. return false;
  210. }
  211. }
  212. #if ONE_OFF_DEBUG
  213. SkDPoint qPt = fQuad.ptAtT(quadT);
  214. SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
  215. qPt.fX, qPt.fY);
  216. #endif
  217. return true;
  218. }
  219. int verticalIntersect(double axisIntercept, double roots[2]) {
  220. double D = fQuad[2].fX; // f
  221. double E = fQuad[1].fX; // e
  222. double F = fQuad[0].fX; // d
  223. D += F - 2 * E; // D = d - 2*e + f
  224. E -= F; // E = -(d - e)
  225. F -= axisIntercept;
  226. return SkDQuad::RootsValidT(D, 2 * E, F, roots);
  227. }
  228. int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
  229. addExactVerticalEndPoints(top, bottom, axisIntercept);
  230. if (fAllowNear) {
  231. addNearVerticalEndPoints(top, bottom, axisIntercept);
  232. }
  233. double rootVals[2];
  234. int roots = verticalIntersect(axisIntercept, rootVals);
  235. for (int index = 0; index < roots; ++index) {
  236. double quadT = rootVals[index];
  237. SkDPoint pt = fQuad.ptAtT(quadT);
  238. double lineT = (pt.fY - top) / (bottom - top);
  239. if (pinTs(&quadT, &lineT, &pt, kPointInitialized) && uniqueAnswer(quadT, pt)) {
  240. fIntersections->insert(quadT, lineT, pt);
  241. }
  242. }
  243. if (flipped) {
  244. fIntersections->flip();
  245. }
  246. checkCoincident();
  247. return fIntersections->used();
  248. }
  249. protected:
  250. // add endpoints first to get zero and one t values exactly
  251. void addExactEndPoints() {
  252. for (int qIndex = 0; qIndex < 3; qIndex += 2) {
  253. double lineT = fLine->exactPoint(fQuad[qIndex]);
  254. if (lineT < 0) {
  255. continue;
  256. }
  257. double quadT = (double) (qIndex >> 1);
  258. fIntersections->insert(quadT, lineT, fQuad[qIndex]);
  259. }
  260. }
  261. void addNearEndPoints() {
  262. for (int qIndex = 0; qIndex < 3; qIndex += 2) {
  263. double quadT = (double) (qIndex >> 1);
  264. if (fIntersections->hasT(quadT)) {
  265. continue;
  266. }
  267. double lineT = fLine->nearPoint(fQuad[qIndex], nullptr);
  268. if (lineT < 0) {
  269. continue;
  270. }
  271. fIntersections->insert(quadT, lineT, fQuad[qIndex]);
  272. }
  273. this->addLineNearEndPoints();
  274. }
  275. void addLineNearEndPoints() {
  276. for (int lIndex = 0; lIndex < 2; ++lIndex) {
  277. double lineT = (double) lIndex;
  278. if (fIntersections->hasOppT(lineT)) {
  279. continue;
  280. }
  281. double quadT = ((SkDCurve*) &fQuad)->nearPoint(SkPath::kQuad_Verb,
  282. (*fLine)[lIndex], (*fLine)[!lIndex]);
  283. if (quadT < 0) {
  284. continue;
  285. }
  286. fIntersections->insert(quadT, lineT, (*fLine)[lIndex]);
  287. }
  288. }
  289. void addExactHorizontalEndPoints(double left, double right, double y) {
  290. for (int qIndex = 0; qIndex < 3; qIndex += 2) {
  291. double lineT = SkDLine::ExactPointH(fQuad[qIndex], left, right, y);
  292. if (lineT < 0) {
  293. continue;
  294. }
  295. double quadT = (double) (qIndex >> 1);
  296. fIntersections->insert(quadT, lineT, fQuad[qIndex]);
  297. }
  298. }
  299. void addNearHorizontalEndPoints(double left, double right, double y) {
  300. for (int qIndex = 0; qIndex < 3; qIndex += 2) {
  301. double quadT = (double) (qIndex >> 1);
  302. if (fIntersections->hasT(quadT)) {
  303. continue;
  304. }
  305. double lineT = SkDLine::NearPointH(fQuad[qIndex], left, right, y);
  306. if (lineT < 0) {
  307. continue;
  308. }
  309. fIntersections->insert(quadT, lineT, fQuad[qIndex]);
  310. }
  311. this->addLineNearEndPoints();
  312. }
  313. void addExactVerticalEndPoints(double top, double bottom, double x) {
  314. for (int qIndex = 0; qIndex < 3; qIndex += 2) {
  315. double lineT = SkDLine::ExactPointV(fQuad[qIndex], top, bottom, x);
  316. if (lineT < 0) {
  317. continue;
  318. }
  319. double quadT = (double) (qIndex >> 1);
  320. fIntersections->insert(quadT, lineT, fQuad[qIndex]);
  321. }
  322. }
  323. void addNearVerticalEndPoints(double top, double bottom, double x) {
  324. for (int qIndex = 0; qIndex < 3; qIndex += 2) {
  325. double quadT = (double) (qIndex >> 1);
  326. if (fIntersections->hasT(quadT)) {
  327. continue;
  328. }
  329. double lineT = SkDLine::NearPointV(fQuad[qIndex], top, bottom, x);
  330. if (lineT < 0) {
  331. continue;
  332. }
  333. fIntersections->insert(quadT, lineT, fQuad[qIndex]);
  334. }
  335. this->addLineNearEndPoints();
  336. }
  337. double findLineT(double t) {
  338. SkDPoint xy = fQuad.ptAtT(t);
  339. double dx = (*fLine)[1].fX - (*fLine)[0].fX;
  340. double dy = (*fLine)[1].fY - (*fLine)[0].fY;
  341. if (fabs(dx) > fabs(dy)) {
  342. return (xy.fX - (*fLine)[0].fX) / dx;
  343. }
  344. return (xy.fY - (*fLine)[0].fY) / dy;
  345. }
  346. bool pinTs(double* quadT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
  347. if (!approximately_one_or_less_double(*lineT)) {
  348. return false;
  349. }
  350. if (!approximately_zero_or_more_double(*lineT)) {
  351. return false;
  352. }
  353. double qT = *quadT = SkPinT(*quadT);
  354. double lT = *lineT = SkPinT(*lineT);
  355. if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && qT != 0 && qT != 1)) {
  356. *pt = (*fLine).ptAtT(lT);
  357. } else if (ptSet == kPointUninitialized) {
  358. *pt = fQuad.ptAtT(qT);
  359. }
  360. SkPoint gridPt = pt->asSkPoint();
  361. if (SkDPoint::ApproximatelyEqual(gridPt, (*fLine)[0].asSkPoint())) {
  362. *pt = (*fLine)[0];
  363. *lineT = 0;
  364. } else if (SkDPoint::ApproximatelyEqual(gridPt, (*fLine)[1].asSkPoint())) {
  365. *pt = (*fLine)[1];
  366. *lineT = 1;
  367. }
  368. if (fIntersections->used() > 0 && approximately_equal((*fIntersections)[1][0], *lineT)) {
  369. return false;
  370. }
  371. if (gridPt == fQuad[0].asSkPoint()) {
  372. *pt = fQuad[0];
  373. *quadT = 0;
  374. } else if (gridPt == fQuad[2].asSkPoint()) {
  375. *pt = fQuad[2];
  376. *quadT = 1;
  377. }
  378. return true;
  379. }
  380. private:
  381. const SkDQuad& fQuad;
  382. const SkDLine* fLine;
  383. SkIntersections* fIntersections;
  384. bool fAllowNear;
  385. };
  386. int SkIntersections::horizontal(const SkDQuad& quad, double left, double right, double y,
  387. bool flipped) {
  388. SkDLine line = {{{ left, y }, { right, y }}};
  389. LineQuadraticIntersections q(quad, line, this);
  390. return q.horizontalIntersect(y, left, right, flipped);
  391. }
  392. int SkIntersections::vertical(const SkDQuad& quad, double top, double bottom, double x,
  393. bool flipped) {
  394. SkDLine line = {{{ x, top }, { x, bottom }}};
  395. LineQuadraticIntersections q(quad, line, this);
  396. return q.verticalIntersect(x, top, bottom, flipped);
  397. }
  398. int SkIntersections::intersect(const SkDQuad& quad, const SkDLine& line) {
  399. LineQuadraticIntersections q(quad, line, this);
  400. q.allowNear(fAllowNear);
  401. return q.intersect();
  402. }
  403. int SkIntersections::intersectRay(const SkDQuad& quad, const SkDLine& line) {
  404. LineQuadraticIntersections q(quad, line, this);
  405. fUsed = q.intersectRay(fT[0]);
  406. for (int index = 0; index < fUsed; ++index) {
  407. fPt[index] = quad.ptAtT(fT[0][index]);
  408. }
  409. return fUsed;
  410. }
  411. int SkIntersections::HorizontalIntercept(const SkDQuad& quad, SkScalar y, double* roots) {
  412. LineQuadraticIntersections q(quad);
  413. return q.horizontalIntersect(y, roots);
  414. }
  415. int SkIntersections::VerticalIntercept(const SkDQuad& quad, SkScalar x, double* roots) {
  416. LineQuadraticIntersections q(quad);
  417. return q.verticalIntersect(x, roots);
  418. }
  419. // SkDQuad accessors to Intersection utilities
  420. int SkDQuad::horizontalIntersect(double yIntercept, double roots[2]) const {
  421. return SkIntersections::HorizontalIntercept(*this, yIntercept, roots);
  422. }
  423. int SkDQuad::verticalIntersect(double xIntercept, double roots[2]) const {
  424. return SkIntersections::VerticalIntercept(*this, xIntercept, roots);
  425. }