-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy path0002.html
More file actions
242 lines (233 loc) · 10.6 KB
/
0002.html
File metadata and controls
242 lines (233 loc) · 10.6 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<link rel="alternate" type="application/atom+xml" title="Atom Feed" href="feed.xml">
<link rel="stylesheet" type="text/css" href="style.css">
<title>intersect_lines() [Part 1] - Lars' Website</title>
</head>
<body>
<header>
<nav>
<h1><a href="index.html">Lars' website</a></h1>
<a href="archive.html" title="archive">Archive</a>
<a href="feed.xml" title="feed">Feed</a>
</nav>
</header>
<article>
<header>
<h2>intersect_lines() [Part 1]</h2>
<div>First published: <time>2021-10-10</time></div>
</header>
<section>
<p>
Recently I looking to get some exercise working with Projective Geometric Algebra.
So I scrolled through my personal codebase looking for something to reimplement, when I saw at my 2D line intersection code (listing 0).
I couldn't remember where or how I got it, and I had probably ported it between languages at least twice since writing it.
</p>
<code id="listing0">
<pre>
<span>inline int line_line_intersection(</span>
<span> Vec2 line0_a, Vec2 line0_b,</span>
<span> Vec2 line1_a, Vec2 line1_b,</span>
<span> Vec2 *result)</span>
<span>{</span>
<span> auto d0 = line0_b - line0_a;</span>
<span> auto d1 = line1_b - line1_a;</span>
<span> auto d2 = line0_a - line1_a;</span>
<span> auto det = cross(d0, d1);</span>
<span> auto s = cross(d0, d2) / det;</span>
<span> auto t = cross(d1, d2) / det;</span>
<span> if (0.0f < s && s < 1.0f && 0.0f < t && t < 1.0f) {</span>
<span> if (result) *result = line0_a + d0 * t;</span>
<span> return 1;</span>
<span> }</span>
<span> return 0;</span>
<span>}</span>
</pre>
<figcaption>Listing 0: It works, but why? And yes that's a 2D cross product.</figcaption>
</code>
<p>
So I started rewriting it, using the obvious operations to join the points into lines and meet the lines to get the intersection.
One caveat being, that the resulting point needs to be normalized, meaning that the e<sub>01</sub> (x) and e<sub>20</sub> (y) component have to be divided by e<sub>12</sub> (z).
But because I'm actually interesting in intersecting line segments, the code also has to check that the intersection happens between the original points, not outside.
The elegant way seemed to be to join the lines with each of the points of the other line, resulting in two oriented areas per line.
If the sign of this areas differed, the points were on opposite sides, if this was true for both line segments, the must be intersecting.
The more efficient code though turned out to be taking the difference from the intersection point to each of the four initial points and then taking the dot product of the the differences per line.
If the intersection was between the points, the difference would point in opposite directions, making the product of dot products negative.
If both products are negative, the intersection happened inside the segments.
</p>
<code id="equation0">
<pre>
Given points a, b, c, d
line l = a join b
line m = c join d
point p = normalize(l meet m)
bool intersects = (a - p) dot (b - p) < 0 && (c - p) dot (d - p) < 0
</pre>
<figcaption>Equation 0: pseudomath for finding a line segment intersection with PGA</figcaption>
</code>
<p>
But after implementing all that, I counted the number of calculations that each algorithm had to do and was surprised to find, that the PGA version was pretty bad.
While my existing code has 11 subs/adds, 8 muls and 2 divs throughout those vector operations, the PGA version has 17 adds/subs, 14 muls and 2 divs.
If you want to count yourself, the "unrolled" code of both versions is below.
</p>
<p>
I started to wonder why the original version was so much better, and how It worked in the first place.
So after building some visualization I found a neat way to explain the algorithm, which I will show in the rest of this article.
</p>
<code id="listing1">
<pre>
<span>inline int line_line_intersection(</span>
<span> Vec2 line0_a, Vec2 line0_b,</span>
<span> Vec2 line1_a, Vec2 line1_b,</span>
<span> Vec2 *result)</span>
<span>{</span>
<span> auto ab_x = line1_b.x - line1_a.x;</span>
<span> auto ab_y = line1_b.y - line1_a.y;</span>
<span> auto cd_x = line2_b.x - line2_a.x;</span>
<span> auto cd_y = line2_b.y - line2_a.y;</span>
<span> auto ca_x = line1_a.x - line2_a.x;</span>
<span> auto ca_y = line1_a.y - line2_a.y;</span>
<span> auto det = ab_x * cd_y - ab_y * cd_x;</span>
<span> auto s = (ab_x * ca_y - ab_y * ca_x) / det;</span>
<span> auto t = (cd_x * ca_y - cd_y * ca_x) / det;</span>
<span> if (0.0f <= s && s <= 1.0f && 0.0f <= t && t <= 1.0f) {</span>
<span> if (result) *result = {</span>
<span> line1_a.x + ab_x * t,</span>
<span> line1_a.y + ab_y * t</span>
<span> };</span>
<span> return 1;</span>
<span> }</span>
<span> return 0;</span>
<span>}</span>
</pre>
<figcaption>Listing 1: Same math as above, but without the vector operations.</figcaption>
</code>
<code id="listing2">
<pre>
<span>inline int line_line_intersection(</span>
<span> Vec2 line0_a, Vec2 line0_b,</span>
<span> Vec2 line1_a, Vec2 line1_b,</span>
<span> Vec2 *result)</span>
<span>{</span>
<span> // Join points a, b in line l</span>
<span> auto l_0 = line1_a.y * line1_b.x - line1_a.x * line1_b.y; // y_a*x_b-x_a*y_b</span>
<span> auto l_1 = line1_b.y - line1_a.y; // y_b-y_a</span>
<span> auto l_2 = line1_a.x - line1_b.x; // x_a-x_b</span>
<span></span>
<span> // Join points c, d in line m</span>
<span> auto m_0 = line2_a.y * line2_b.x - line2_a.x * line2_b.y; // y_c*x_d-x_c*y_d</span>
<span> auto m_1 = line2_b.y - line2_a.y; // y_d-y_c</span>
<span> auto m_2 = line2_a.x - line2_b.x; // x_c-x_d</span>
<span></span>
<span> // Meet lines l, m in point p</span>
<span> auto p_12 = m_2 * l_1 - m_1 * l_2;</span>
<span> auto p_01 = (m_1 * l_0 - m_0 * l_1) / p_12;</span>
<span> auto p_20 = (m_0 * l_2 - m_2 * l_0) / p_12;</span>
<span></span>
<span> auto s = (line1_a.x - p_20) * (line1_b.x - p_20) + (line1_a.y - p_01) * (line1_b.y - p_01); </span>
<span> auto t = (line2_a.x - p_20) * (line2_b.x - p_20) + (line2_a.y - p_01) * (line2_b.y - p_01); </span>
<span></span>
<span> if (s <= 0.0f && t <= 0.0f) {</span>
<span> if (result) {</span>
<span> *result = { p_20, p_01 };</span>
<span> }</span>
<span> return 1;</span>
<span> }</span>
<span> return 0;</span>
<span>}</span>
</pre>
<figcaption>Listing 2: The PGA version in all its glory.</figcaption>
</code>
</section>
<section>
<h3>The easiest way to find an intersection between line segments</h3>
<p>
I was about to say that this algorithm is so easy, you could do it by hand.
But then I was informed, that there exists an even easier way to tell if two lines intersect, called "by eye".
So let's just stick to computers and algorithms.
</p>
<code id="figure0">
<pre>
A │ C
╲ │ │
╲ │ B───┼─A
C───D ╲ │ │
B │ D
</pre>
<figcaption>Figure 0: Two non-intersecting and two intersecting lines, each made of two points A, B and C, D.</figcaption>
</code>
<p>
We have here two line segments, defined by their endpoints: A, B and C, D.
The first step to the algorithm is to turn the line segments into differences: A-B and D-C.
Another way to look at is that we move both line segment so that one of their ends lines up with the origin.
</p>
<code id="figure1">
<pre>
(A-B) │ O─────(B-A)
╲ │ │
╲ │ │
╲ │ │
O───(D-C) │ (D-C)
</pre>
<figcaption>Figure 1</figcaption>
</code>
<p>
Now we can duplicate both lines and line them up, so that we get a parallelogram.
Another way to look at it is, that we add the two differences, ending up with A-B+D-C.
</p>
<code id="figure2">
<pre>
(A-B)───(A-B+D-C) │ O─────(A-B)
╲ ╲ │ │ │
╲ ╲ │ │ │
╲ ╲ │ │ │
O───(D-C) │ (D-C)─────(A-B+D-C)
</pre>
<figcaption>Figure 2</figcaption>
</code>
<p>
And this is where the magic happens: if we now take a third difference A-C, we can check if it lands inside the parallelogram or not.
If it lands outside, as is the case on the left side, the line segments did not intersect.
But if they do intersect, as is the case on the right side, A-C will land inside the parallelogram.
</p>
<code id="figure3">
<pre>
(A-B)───(A-B+D-C) │ O─────(A-B)
╲ ╲ (A-C) │ │╲ │
╲ ╳ │ │ (A-C)
╲ ╱ ╲ │ │ │
O───(D-C) │ (D-C)─────(A-B+D-C)
</pre>
<figcaption>Figure 3</figcaption>
</code>
<p>
Now to be honest, I had to fiddle around a bit before finding which differences to take to get to this result, and they are not quite the same as in the original code.
But the math can be made to work either way and this was in my opinion the best way to visualize it.
</p>
</section>
<section>
<h3>Minkowski at it again</h3>
<p>
While writing this post, I realized, that what I'm doing here is very similar to taking the Minkowski difference of the two lines.
So there will be a part 2 in the future, exploring that side of the algorithm as well.
I also have yet to explained how to get from the parallelogram to the point of intersection, and I might also have some more to say about the PGA version.
Cheers!
</p>
</section>
<section>
<h2>Change-log</h2>
<ul>
<li>2022-06-07 - Switched from pure ASCII to Unicode box-drawing characters.</li>
<li>2021-10-12 - Improved formulations.</li>
<li>2021-10-10 - Initial post.</li>
</ul>
</section>
</article>
<footer>
<a href="changelog.html">Change-log</a>
</footer>
</body>
</html>