experimental/Intersection/QuadraticImplicit.cpp - Issue 867213004: remove prototype pathops code

Side by Side Diff: experimental/Intersection/QuadraticImplicit.cpp

Issue 867213004: remove prototype pathops code (Closed) Base URL: https://skia.googlesource.com/skia.git@master

Patch Set: Created 5 years, 10 months ago

Use n/p to move between diff chunks; N/P to move between comments. Draft comments are only viewable by you.

Jump to:

OLD	NEW
	(Empty)
1 // Another approach is to start with the implicit form of one curve and solve

2 // (seek implicit coefficients in QuadraticParameter.cpp

3 // by substituting in the parametric form of the other.

4 // The downside of this approach is that early rejects are difficult to come by.

5 // http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormu la.html#step

6

7

8 #include "CubicUtilities.h"

9 #include "CurveIntersection.h"

10 #include "Intersections.h"

11 #include "QuadraticParameterization.h"

12 #include "QuarticRoot.h"

13 #include "QuadraticUtilities.h"

14 #include "TSearch.h"

15

16 #ifdef SK_DEBUG

17 #include "LineUtilities.h"

18 #endif

19

20 /* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F

21 * and given x = at^2 + bt + c (the parameterized form)

22 * y = dt^2 + et + f

23 * then

24 * 0 = A(at^2+bt+c)(at^2+bt+c)+B(at^2+bt+c)(dt^2+et+f)+C(dt^2+et+f)(dt^2+et+f)+D (at^2+bt+c)+E(dt^2+et+f)+F

25 */

26

27 static int findRoots(const QuadImplicitForm& i, const Quadratic& q2, double root s[4],

28 bool oneHint, int firstCubicRoot) {

29 double a, b, c;

30 set_abc(&q2[0].x, a, b, c);

31 double d, e, f;

32 set_abc(&q2[0].y, d, e, f);

33 const double t4 = i.x2() * a * a

34 + i.xy() * a * d

35 + i.y2() * d * d;

36 const double t3 = 2 * i.x2() * a * b

37 + i.xy() * (a * e + b * d)

38 + 2 * i.y2() * d * e;

39 const double t2 = i.x2() * (b * b + 2 * a * c)

40 + i.xy() * (c * d + b * e + a * f)

41 + i.y2() * (e * e + 2 * d * f)

42 + i.x() * a

43 + i.y() * d;

44 const double t1 = 2 * i.x2() * b * c

45 + i.xy() * (c * e + b * f)

46 + 2 * i.y2() * e * f

47 + i.x() * b

48 + i.y() * e;

49 const double t0 = i.x2() * c * c

50 + i.xy() * c * f

51 + i.y2() * f * f

52 + i.x() * c

53 + i.y() * f

54 + i.c();

55 int rootCount = reducedQuarticRoots(t4, t3, t2, t1, t0, oneHint, roots);

56 if (rootCount >= 0) {

57 return rootCount;

58 }

59 return quarticRootsReal(firstCubicRoot, t4, t3, t2, t1, t0, roots);

60 }

61

62 static int addValidRoots(const double roots[4], const int count, double valid[4] ) {

63 int result = 0;

64 int index;

65 for (index = 0; index < count; ++index) {

66 if (!approximately_zero_or_more(roots[index]) \|\| !approximately_one_or_l ess(roots[index])) {

67 continue;

68 }

69 double t = 1 - roots[index];

70 if (approximately_less_than_zero(t)) {

71 t = 0;

72 } else if (approximately_greater_than_one(t)) {

73 t = 1;

74 }

75 valid[result++] = t;

76 }

77 return result;

78 }

79

80 static bool onlyEndPtsInCommon(const Quadratic& q1, const Quadratic& q2, Interse ctions& i) {

81 // the idea here is to see at minimum do a quick reject by rotating all points

82 // to either side of the line formed by connecting the endpoints

83 // if the opposite curves points are on the line or on the other side, the

84 // curves at most intersect at the endpoints

85 for (int oddMan = 0; oddMan < 3; ++oddMan) {

86 const _Point* endPt[2];

87 for (int opp = 1; opp < 3; ++opp) {

88 int end = oddMan ^ opp;

89 if (end == 3) {

90 end = opp;

91 }

92 endPt[opp - 1] = &q1[end];

93 }

94 double origX = endPt[0]->x;

95 double origY = endPt[0]->y;

96 double adj = endPt[1]->x - origX;

97 double opp = endPt[1]->y - origY;

98 double sign = (q1[oddMan].y - origY) * adj - (q1[oddMan].x - origX) * op p;

99 if (approximately_zero(sign)) {

100 goto tryNextHalfPlane;

101 }

102 for (int n = 0; n < 3; ++n) {

103 double test = (q2[n].y - origY) * adj - (q2[n].x - origX) * opp;

104 if (test * sign > 0) {

105 goto tryNextHalfPlane;

106 }

107 }

108 for (int i1 = 0; i1 < 3; i1 += 2) {

109 for (int i2 = 0; i2 < 3; i2 += 2) {

110 if (q1[i1] == q2[i2]) {

111 i.insert(i1 >> 1, i2 >> 1, q1[i1]);

112 }

113 }

114 }

115 SkASSERT(i.fUsed < 3);

116 return true;

117 tryNextHalfPlane:

118 ;

119 }

120 return false;

121 }

122

123 // returns false if there's more than one intercept or the intercept doesn't mat ch the point

124 // returns true if the intercept was successfully added or if the

125 // original quads need to be subdivided

126 static bool addIntercept(const Quadratic& q1, const Quadratic& q2, double tMin, double tMax,

127 Intersections& i, bool* subDivide) {

128 double tMid = (tMin + tMax) / 2;

129 _Point mid;

130 xy_at_t(q2, tMid, mid.x, mid.y);

131 _Line line;

132 line[0] = line[1] = mid;

133 _Vector dxdy = dxdy_at_t(q2, tMid);

134 line[0] -= dxdy;

135 line[1] += dxdy;

136 Intersections rootTs;

137 int roots = intersect(q1, line, rootTs);

138 if (roots == 0) {

139 if (subDivide) {

140 *subDivide = true;

141 }

142 return true;

143 }

144 if (roots == 2) {

145 return false;

146 }

147 _Point pt2;

148 xy_at_t(q1, rootTs.fT[0][0], pt2.x, pt2.y);

149 if (!pt2.approximatelyEqualHalf(mid)) {

150 return false;

151 }

152 i.insertSwap(rootTs.fT[0][0], tMid, pt2);

153 return true;

154 }

155

156 static bool isLinearInner(const Quadratic& q1, double t1s, double t1e, const Qua dratic& q2,

157 double t2s, double t2e, Intersections& i, bool* subDivide) {

158 Quadratic hull;

159 sub_divide(q1, t1s, t1e, hull);

160 _Line line = {hull[2], hull[0]};

161 const _Line* testLines[] = { &line, (const _Line) &hull[0], (const _Line) &hull[1] };

162 size_t testCount = sizeof(testLines) / sizeof(testLines[0]);

163 SkTDArray<double> tsFound;

164 for (size_t index = 0; index < testCount; ++index) {

165 Intersections rootTs;

166 int roots = intersect(q2, *testLines[index], rootTs);

167 for (int idx2 = 0; idx2 < roots; ++idx2) {

168 double t = rootTs.fT[0][idx2];

169 #ifdef SK_DEBUG

170 _Point qPt, lPt;

171 xy_at_t(q2, t, qPt.x, qPt.y);

172 xy_at_t(*testLines[index], rootTs.fT[1][idx2], lPt.x, lPt.y);

173 SkASSERT(qPt.approximatelyEqual(lPt));

174 #endif

175 if (approximately_negative(t - t2s) \|\| approximately_positive(t - t2 e)) {

176 continue;

177 }

178 *tsFound.append() = rootTs.fT[0][idx2];

179 }

180 }

181 int tCount = tsFound.count();

182 if (!tCount) {

183 return true;

184 }

185 double tMin, tMax;

186 if (tCount == 1) {

187 tMin = tMax = tsFound[0];

188 } else if (tCount > 1) {

189 QSort<double>(tsFound.begin(), tsFound.end() - 1);

190 tMin = tsFound[0];

191 tMax = tsFound[tsFound.count() - 1];

192 }

193 _Point end;

194 xy_at_t(q2, t2s, end.x, end.y);

195 bool startInTriangle = point_in_hull(hull, end);

196 if (startInTriangle) {

197 tMin = t2s;

198 }

199 xy_at_t(q2, t2e, end.x, end.y);

200 bool endInTriangle = point_in_hull(hull, end);

201 if (endInTriangle) {

202 tMax = t2e;

203 }

204 int split = 0;

205 _Vector dxy1, dxy2;

206 if (tMin != tMax \|\| tCount > 2) {

207 dxy2 = dxdy_at_t(q2, tMin);

208 for (int index = 1; index < tCount; ++index) {

209 dxy1 = dxy2;

210 dxy2 = dxdy_at_t(q2, tsFound[index]);

211 double dot = dxy1.dot(dxy2);

212 if (dot < 0) {

213 split = index - 1;

214 break;

215 }

216 }

217

218 }

219 if (split == 0) { // there's one point

220 if (addIntercept(q1, q2, tMin, tMax, i, subDivide)) {

221 return true;

222 }

223 i.swap();

224 return isLinearInner(q2, tMin, tMax, q1, t1s, t1e, i, subDivide);

225 }

226 // At this point, we have two ranges of t values -- treat each separately at the split

227 bool result;

228 if (addIntercept(q1, q2, tMin, tsFound[split - 1], i, subDivide)) {

229 result = true;

230 } else {

231 i.swap();

232 result = isLinearInner(q2, tMin, tsFound[split - 1], q1, t1s, t1e, i, su bDivide);

233 }

234 if (addIntercept(q1, q2, tsFound[split], tMax, i, subDivide)) {

235 result = true;

236 } else {

237 i.swap();

238 result \|= isLinearInner(q2, tsFound[split], tMax, q1, t1s, t1e, i, subDi vide);

239 }

240 return result;

241 }

242

243 static double flatMeasure(const Quadratic& q) {

244 _Vector mid = q[1] - q[0];

245 _Vector dxy = q[2] - q[0];

246 double length = dxy.length(); // OPTIMIZE: get rid of sqrt

247 return fabs(mid.cross(dxy) / length);

248 }

249

250 // FIXME ? should this measure both and then use the quad that is the flattest a s the line?

251 static bool isLinear(const Quadratic& q1, const Quadratic& q2, Intersections& i) {

252 double measure = flatMeasure(q1);

253 // OPTIMIZE: (get rid of sqrt) use approximately_zero

254 if (!approximately_zero_sqrt(measure)) {

255 return false;

256 }

257 return isLinearInner(q1, 0, 1, q2, 0, 1, i, NULL);

258 }

259

260 // FIXME: if flat measure is sufficiently large, then probably the quartic solut ion failed

261 static void relaxedIsLinear(const Quadratic& q1, const Quadratic& q2, Intersecti ons& i) {

262 double m1 = flatMeasure(q1);

263 double m2 = flatMeasure(q2);

264 #ifdef SK_DEBUG

265 double min = SkTMin(m1, m2);

266 if (min > 5) {

267 SkDebugf("%s maybe not flat enough.. %1.9g\n", __FUNCTION__, min);

268 }

269 #endif

270 i.reset();

271 const Quadratic& rounder = m2 < m1 ? q1 : q2;

272 const Quadratic& flatter = m2 < m1 ? q2 : q1;

273 bool subDivide = false;

274 isLinearInner(flatter, 0, 1, rounder, 0, 1, i, &subDivide);

275 if (subDivide) {

276 QuadraticPair pair;

277 chop_at(flatter, pair, 0.5);

278 Intersections firstI, secondI;

279 relaxedIsLinear(pair.first(), rounder, firstI);

280 for (int index = 0; index < firstI.used(); ++index) {

281 i.insert(firstI.fT[0][index] * 0.5, firstI.fT[1][index], firstI.fPt[ index]);

282 }

283 relaxedIsLinear(pair.second(), rounder, secondI);

284 for (int index = 0; index < secondI.used(); ++index) {

285 i.insert(0.5 + secondI.fT[0][index] * 0.5, secondI.fT[1][index], sec ondI.fPt[index]);

286 }

287 }

288 if (m2 < m1) {

289 i.swapPts();

290 }

291 }

292

293 #if 0

294 static void unsortableExpanse(const Quadratic& q1, const Quadratic& q2, Intersec tions& i) {

295 const Quadratic* qs[2] = { &q1, &q2 };

296 // need t values for start and end of unsortable expanse on both curves

297 // try projecting lines parallel to the end points

298 i.fT[0][0] = 0;

299 i.fT[0][1] = 1;

300 int flip = -1; // undecided

301 for (int qIdx = 0; qIdx < 2; qIdx++) {

302 for (int t = 0; t < 2; t++) {

303 _Point dxdy;

304 dxdy_at_t(*qs[qIdx], t, dxdy);

305 _Line perp;

306 perp[0] = perp[1] = (*qs[qIdx])[t == 0 ? 0 : 2];

307 perp[0].x += dxdy.y;

308 perp[0].y -= dxdy.x;

309 perp[1].x -= dxdy.y;

310 perp[1].y += dxdy.x;

311 Intersections hitData;

312 int hits = intersectRay(*qs[qIdx ^ 1], perp, hitData);

313 SkASSERT(hits <= 1);

314 if (hits) {

315 if (flip < 0) {

316 _Point dxdy2;

317 dxdy_at_t(*qs[qIdx ^ 1], hitData.fT[0][0], dxdy2);

318 double dot = dxdy.dot(dxdy2);

319 flip = dot < 0;

320 i.fT[1][0] = flip;

321 i.fT[1][1] = !flip;

322 }

323 i.fT[qIdx ^ 1][t ^ flip] = hitData.fT[0][0];

324 }

325 }

326 }

327 i.fUnsortable = true; // failed, probably coincident or near-coincident

328 i.fUsed = 2;

329 }

330 #endif

331

332 // each time through the loop, this computes values it had from the last loop

333 // if i == j == 1, the center values are still good

334 // otherwise, for i != 1 or j != 1, four of the values are still good

335 // and if i == 1 ^ j == 1, an additional value is good

336 static bool binarySearch(const Quadratic& quad1, const Quadratic& quad2, double& t1Seed,

337 double& t2Seed, _Point& pt) {

338 double tStep = ROUGH_EPSILON;

339 _Point t1[3], t2[3];

340 int calcMask = ~0;

341 do {

342 if (calcMask & (1 << 1)) t1[1] = xy_at_t(quad1, t1Seed);

343 if (calcMask & (1 << 4)) t2[1] = xy_at_t(quad2, t2Seed);

344 if (t1[1].approximatelyEqual(t2[1])) {

345 pt = t1[1];

346 #if ONE_OFF_DEBUG

347 SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) == (%1.9g,%1.9g)\n", __ FUNCTION__,

348 t1Seed, t2Seed, t1[1].x, t1[1].y, t1[2].x, t1[2].y);

349 #endif

350 return true;

351 }

352 if (calcMask & (1 << 0)) t1[0] = xy_at_t(quad1, t1Seed - tStep);

353 if (calcMask & (1 << 2)) t1[2] = xy_at_t(quad1, t1Seed + tStep);

354 if (calcMask & (1 << 3)) t2[0] = xy_at_t(quad2, t2Seed - tStep);

355 if (calcMask & (1 << 5)) t2[2] = xy_at_t(quad2, t2Seed + tStep);

356 double dist[3][3];

357 // OPTIMIZE: using calcMask value permits skipping some distance calcuat ions

358 // if prior loop's results are moved to correct slot for reuse

359 dist[1][1] = t1[1].distanceSquared(t2[1]);

360 int best_i = 1, best_j = 1;

361 for (int i = 0; i < 3; ++i) {

362 for (int j = 0; j < 3; ++j) {

363 if (i == 1 && j == 1) {

364 continue;

365 }

366 dist[i][j] = t1[i].distanceSquared(t2[j]);

367 if (dist[best_i][best_j] > dist[i][j]) {

368 best_i = i;

369 best_j = j;

370 }

371 }

372 }

373 if (best_i == 1 && best_j == 1) {

374 tStep /= 2;

375 if (tStep < FLT_EPSILON_HALF) {

376 break;

377 }

378 calcMask = (1 << 0) \| (1 << 2) \| (1 << 3) \| (1 << 5);

379 continue;

380 }

381 if (best_i == 0) {

382 t1Seed -= tStep;

383 t1[2] = t1[1];

384 t1[1] = t1[0];

385 calcMask = 1 << 0;

386 } else if (best_i == 2) {

387 t1Seed += tStep;

388 t1[0] = t1[1];

389 t1[1] = t1[2];

390 calcMask = 1 << 2;

391 } else {

392 calcMask = 0;

393 }

394 if (best_j == 0) {

395 t2Seed -= tStep;

396 t2[2] = t2[1];

397 t2[1] = t2[0];

398 calcMask \|= 1 << 3;

399 } else if (best_j == 2) {

400 t2Seed += tStep;

401 t2[0] = t2[1];

402 t2[1] = t2[2];

403 calcMask \|= 1 << 5;

404 }

405 } while (true);

406 #if ONE_OFF_DEBUG

407 SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) != (%1.9g,%1.9g) %s\n", __FUNCT ION__,

408 t1Seed, t2Seed, t1[1].x, t1[1].y, t1[2].x, t1[2].y);

409 #endif

410 return false;

411 }

412

413 bool intersect2(const Quadratic& q1, const Quadratic& q2, Intersections& i) {

414 // if the quads share an end point, check to see if they overlap

415

416 if (onlyEndPtsInCommon(q1, q2, i)) {

417 return i.intersected();

418 }

419 if (onlyEndPtsInCommon(q2, q1, i)) {

420 i.swapPts();

421 return i.intersected();

422 }

423 // see if either quad is really a line

424 if (isLinear(q1, q2, i)) {

425 return i.intersected();

426 }

427 if (isLinear(q2, q1, i)) {

428 i.swapPts();

429 return i.intersected();

430 }

431 QuadImplicitForm i1(q1);

432 QuadImplicitForm i2(q2);

433 if (i1.implicit_match(i2)) {

434 // FIXME: compute T values

435 // compute the intersections of the ends to find the coincident span

436 bool useVertical = fabs(q1[0].x - q1[2].x) < fabs(q1[0].y - q1[2].y);

437 double t;

438 if ((t = axialIntersect(q1, q2[0], useVertical)) >= 0) {

439 i.insertCoincident(t, 0, q2[0]);

440 }

441 if ((t = axialIntersect(q1, q2[2], useVertical)) >= 0) {

442 i.insertCoincident(t, 1, q2[2]);

443 }

444 useVertical = fabs(q2[0].x - q2[2].x) < fabs(q2[0].y - q2[2].y);

445 if ((t = axialIntersect(q2, q1[0], useVertical)) >= 0) {

446 i.insertCoincident(0, t, q1[0]);

447 }

448 if ((t = axialIntersect(q2, q1[2], useVertical)) >= 0) {

449 i.insertCoincident(1, t, q1[2]);

450 }

451 SkASSERT(i.coincidentUsed() <= 2);

452 return i.coincidentUsed() > 0;

453 }

454 int index;

455 bool useCubic = q1[0] == q2[0] \|\| q1[0] == q2[2] \|\| q1[2] == q2[0];

456 double roots1[4];

457 int rootCount = findRoots(i2, q1, roots1, useCubic, 0);

458 // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1

459 double roots1Copy[4];

460 int r1Count = addValidRoots(roots1, rootCount, roots1Copy);

461 _Point pts1[4];

462 for (index = 0; index < r1Count; ++index) {

463 xy_at_t(q1, roots1Copy[index], pts1[index].x, pts1[index].y);

464 }

465 double roots2[4];

466 int rootCount2 = findRoots(i1, q2, roots2, useCubic, 0);

467 double roots2Copy[4];

468 int r2Count = addValidRoots(roots2, rootCount2, roots2Copy);

469 _Point pts2[4];

470 for (index = 0; index < r2Count; ++index) {

471 xy_at_t(q2, roots2Copy[index], pts2[index].x, pts2[index].y);

472 }

473 if (r1Count == r2Count && r1Count <= 1) {

474 if (r1Count == 1) {

475 if (pts1[0].approximatelyEqualHalf(pts2[0])) {

476 i.insert(roots1Copy[0], roots2Copy[0], pts1[0]);

477 } else if (pts1[0].moreRoughlyEqual(pts2[0])) {

478 // experiment: see if a different cubic solution provides the co rrect quartic answer

479 #if 0

480 for (int cu1 = 0; cu1 < 3; ++cu1) {

481 rootCount = findRoots(i2, q1, roots1, useCubic, cu1);

482 r1Count = addValidRoots(roots1, rootCount, roots1Copy);

483 if (r1Count == 0) {

484 continue;

485 }

486 for (int cu2 = 0; cu2 < 3; ++cu2) {

487 if (cu1 == 0 && cu2 == 0) {

488 continue;

489 }

490 rootCount2 = findRoots(i1, q2, roots2, useCubic, cu2);

491 r2Count = addValidRoots(roots2, rootCount2, roots2Copy);

492 if (r2Count == 0) {

493 continue;

494 }

495 SkASSERT(r1Count == 1 && r2Count == 1);

496 SkDebugf("*** [%d,%d] (%1.9g,%1.9g) %s (%1.9g,%1.9g)\n", cu1, cu2,

497 pts1[0].x, pts1[0].y, pts1[0].approximatelyEqual Half(pts2[0])

498 ? "==" : "!=", pts2[0].x, pts2[0].y);

499 }

500 }

501 #endif

502 // experiment: try to find intersection by chasing t

503 rootCount = findRoots(i2, q1, roots1, useCubic, 0);

504 r1Count = addValidRoots(roots1, rootCount, roots1Copy);

505 rootCount2 = findRoots(i1, q2, roots2, useCubic, 0);

506 r2Count = addValidRoots(roots2, rootCount2, roots2Copy);

507 if (binarySearch(q1, q2, roots1Copy[0], roots2Copy[0], pts1[0])) {

508 i.insert(roots1Copy[0], roots2Copy[0], pts1[0]);

509 }

510 }

511 }

512 return i.intersected();

513 }

514 int closest[4];

515 double dist[4];

516 bool foundSomething = false;

517 for (index = 0; index < r1Count; ++index) {

518 dist[index] = DBL_MAX;

519 closest[index] = -1;

520 for (int ndex2 = 0; ndex2 < r2Count; ++ndex2) {

521 if (!pts2[ndex2].approximatelyEqualHalf(pts1[index])) {

522 continue;

523 }

524 double dx = pts2[ndex2].x - pts1[index].x;

525 double dy = pts2[ndex2].y - pts1[index].y;

526 double distance = dx * dx + dy * dy;

527 if (dist[index] <= distance) {

528 continue;

529 }

530 for (int outer = 0; outer < index; ++outer) {

531 if (closest[outer] != ndex2) {

532 continue;

533 }

534 if (dist[outer] < distance) {

535 goto next;

536 }

537 closest[outer] = -1;

538 }

539 dist[index] = distance;

540 closest[index] = ndex2;

541 foundSomething = true;

542 next:

543 ;

544 }

545 }

546 if (r1Count && r2Count && !foundSomething) {

547 relaxedIsLinear(q1, q2, i);

548 return i.intersected();

549 }

550 int used = 0;

551 do {

552 double lowest = DBL_MAX;

553 int lowestIndex = -1;

554 for (index = 0; index < r1Count; ++index) {

555 if (closest[index] < 0) {

556 continue;

557 }

558 if (roots1Copy[index] < lowest) {

559 lowestIndex = index;

560 lowest = roots1Copy[index];

561 }

562 }

563 if (lowestIndex < 0) {

564 break;

565 }

566 i.insert(roots1Copy[lowestIndex], roots2Copy[closest[lowestIndex]],

567 pts1[lowestIndex]);

568 closest[lowestIndex] = -1;

569 } while (++used < r1Count);

570 i.fFlip = false;

571 return i.intersected();

572 }

OLD	NEW

« no previous file with comments | « experimental/Intersection/QuadraticBounds.cpp ('k') | experimental/Intersection/QuadraticIntersection.cpp » ('j') | no next file with comments »