https://www.particleincell.com/2013/cubic-line-intersection/
(zs Intersections, l0, l1, p0, p1, p2, p3 Point)
| 421 | |
| 422 | // https://www.particleincell.com/2013/cubic-line-intersection/ |
| 423 | func intersectionLineCube(zs Intersections, l0, l1, p0, p1, p2, p3 Point) Intersections { |
| 424 | if l0.Equals(l1) { |
| 425 | return zs // zero-length Close |
| 426 | } |
| 427 | |
| 428 | // write line as A.X = bias |
| 429 | A := Point{l1.Y - l0.Y, l0.X - l1.X} |
| 430 | bias := l0.Dot(A) |
| 431 | |
| 432 | a := A.Dot(p3.Sub(p0).Add(p1.Mul(3.0)).Sub(p2.Mul(3.0))) |
| 433 | b := A.Dot(p0.Mul(3.0).Sub(p1.Mul(6.0)).Add(p2.Mul(3.0))) |
| 434 | c := A.Dot(p1.Mul(3.0).Sub(p0.Mul(3.0))) |
| 435 | d := A.Dot(p0) - bias |
| 436 | |
| 437 | roots := []float64{} |
| 438 | r0, r1, r2 := solveCubicFormula(a, b, c, d) |
| 439 | if !math.IsNaN(r0) { |
| 440 | roots = append(roots, r0) |
| 441 | if !math.IsNaN(r1) { |
| 442 | roots = append(roots, r1) |
| 443 | if !math.IsNaN(r2) { |
| 444 | roots = append(roots, r2) |
| 445 | } |
| 446 | } |
| 447 | } |
| 448 | |
| 449 | dira := l1.Sub(l0).Angle() |
| 450 | horizontal := math.Abs(l1.Y-l0.Y) <= math.Abs(l1.X-l0.X) |
| 451 | for _, root := range roots { |
| 452 | if Interval(root, 0.0, 1.0) { |
| 453 | var s float64 |
| 454 | pos := cubicBezierPos(p0, p1, p2, p3, root) |
| 455 | if horizontal { |
| 456 | s = (pos.X - l0.X) / (l1.X - l0.X) |
| 457 | } else { |
| 458 | s = (pos.Y - l0.Y) / (l1.Y - l0.Y) |
| 459 | } |
| 460 | if Interval(s, 0.0, 1.0) { |
| 461 | deriv := cubicBezierDeriv(p0, p1, p2, p3, root) |
| 462 | dirb := deriv.Angle() |
| 463 | tangent := Equal(A.Dot(deriv), 0.0) |
| 464 | endpoint := Equal(root, 0.0) || Equal(root, 1.0) || Equal(s, 0.0) || Equal(s, 1.0) |
| 465 | if endpoint { |
| 466 | // deviate angle slightly at endpoint when aligned to properly set Into |
| 467 | deriv2 := cubicBezierDeriv2(p0, p1, p2, p3, root) |
| 468 | if (0.0 <= deriv.PerpDot(deriv2)) == (Equal(root, 0.0) || !Equal(root, 1.0) && Equal(s, 0.0)) { |
| 469 | dirb += Epsilon * 2.0 // t=0 and CCW, or t=1 and CW |
| 470 | } else { |
| 471 | dirb -= Epsilon * 2.0 // t=0 and CW, or t=1 and CCW |
| 472 | } |
| 473 | } else if angleEqual(dira, dirb) || angleEqual(dira, dirb+math.Pi) { |
| 474 | // directions are parallel but the paths do cross (inflection point) |
| 475 | // TODO: test better |
| 476 | deriv2 := cubicBezierDeriv2(p0, p1, p2, p3, root) |
| 477 | if Equal(deriv2.X, 0.0) && Equal(deriv2.Y, 0.0) { |
| 478 | deriv3 := cubicBezierDeriv3(p0, p1, p2, p3, root) |
| 479 | if 0.0 < deriv.PerpDot(deriv3) { |
| 480 | dirb += Epsilon * 2.0 |