1// Copyright 2012 Google Inc.
2// Copyright 2020 Yevhenii Reizner
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/*
8Find the intersection of a line and cubic by solving for valid t values.
9
10Analogous to line-quadratic intersection, solve line-cubic intersection by
11representing the cubic as:
12 x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3
13 y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3
14and the line as:
15 y = i*x + j (if the line is more horizontal)
16or:
17 x = i*y + j (if the line is more vertical)
18
19Then using Mathematica, solve for the values of t where the cubic intersects the
20line:
21
22 (in) Resultant[
23 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x,
24 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x]
25 (out) -e + j +
26 3 e t - 3 f t -
27 3 e t^2 + 6 f t^2 - 3 g t^2 +
28 e t^3 - 3 f t^3 + 3 g t^3 - h t^3 +
29 i ( a -
30 3 a t + 3 b t +
31 3 a t^2 - 6 b t^2 + 3 c t^2 -
32 a t^3 + 3 b t^3 - 3 c t^3 + d t^3 )
33
34if i goes to infinity, we can rewrite the line in terms of x. Mathematica:
35
36 (in) Resultant[
37 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j,
38 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
39 (out) a - j -
40 3 a t + 3 b t +
41 3 a t^2 - 6 b t^2 + 3 c t^2 -
42 a t^3 + 3 b t^3 - 3 c t^3 + d t^3 -
43 i ( e -
44 3 e t + 3 f t +
45 3 e t^2 - 6 f t^2 + 3 g t^2 -
46 e t^3 + 3 f t^3 - 3 g t^3 + h t^3 )
47
48Solving this with Mathematica produces an expression with hundreds of terms;
49instead, use Numeric Solutions recipe to solve the cubic.
50
51The near-horizontal case, in terms of: Ax^3 + Bx^2 + Cx + D == 0
52 A = (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d) )
53 B = 3*(-( e - 2*f + g ) + i*( a - 2*b + c ) )
54 C = 3*(-(-e + f ) + i*(-a + b ) )
55 D = (-( e ) + i*( a ) + j )
56
57The near-vertical case, in terms of: Ax^3 + Bx^2 + Cx + D == 0
58 A = ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h) )
59 B = 3*( ( a - 2*b + c ) - i*( e - 2*f + g ) )
60 C = 3*( (-a + b ) - i*(-e + f ) )
61 D = ( ( a ) - i*( e ) - j )
62
63For horizontal lines:
64(in) Resultant[
65 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j,
66 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
67(out) e - j -
68 3 e t + 3 f t +
69 3 e t^2 - 6 f t^2 + 3 g t^2 -
70 e t^3 + 3 f t^3 - 3 g t^3 + h t^3
71*/
72
73use super::cubic64::{self, Cubic64};
74use super::point64::SearchAxis;
75use super::Scalar64;
76
77pub fn horizontal_intersect(cubic: &Cubic64, axis_intercept: f64, roots: &mut [f64; 3]) -> usize {
78 let (a, b, c, mut d) = cubic64::coefficients(&cubic.as_f64_slice()[1..]);
79 d -= axis_intercept;
80 let mut count = cubic64::roots_valid_t(a, b, c, d, roots);
81 let mut index = 0;
82 while index < count {
83 let calc_pt = cubic.point_at_t(roots[index]);
84 if !calc_pt.y.approximately_equal(axis_intercept) {
85 let mut extreme_ts = [0.0; 6];
86 let extrema = cubic64::find_extrema(&cubic.as_f64_slice()[1..], &mut extreme_ts);
87 count = cubic.search_roots(
88 extrema,
89 axis_intercept,
90 SearchAxis::Y,
91 &mut extreme_ts,
92 roots,
93 );
94 break;
95 }
96
97 index += 1;
98 }
99
100 count
101}
102
103pub fn vertical_intersect(cubic: &Cubic64, axis_intercept: f64, roots: &mut [f64; 3]) -> usize {
104 let (a, b, c, mut d) = cubic64::coefficients(&cubic.as_f64_slice());
105 d -= axis_intercept;
106 let mut count = cubic64::roots_valid_t(a, b, c, d, roots);
107 let mut index = 0;
108 while index < count {
109 let calc_pt = cubic.point_at_t(roots[index]);
110 if !calc_pt.x.approximately_equal(axis_intercept) {
111 let mut extreme_ts = [0.0; 6];
112 let extrema = cubic64::find_extrema(&cubic.as_f64_slice(), &mut extreme_ts);
113 count = cubic.search_roots(
114 extrema,
115 axis_intercept,
116 SearchAxis::X,
117 &mut extreme_ts,
118 roots,
119 );
120 break;
121 }
122
123 index += 1;
124 }
125
126 count
127}
128