GCC Code Coverage Report


Directory: ./
File: src/2geom/line.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 53 269 19.7%
Functions: 6 21 28.6%
Branches: 33 514 6.4%

Line Branch Exec Source
1 /*
2 * Infinite Straight Line
3 *
4 * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>
5 * Nathan Hurst
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it either under the terms of the GNU Lesser General Public
9 * License version 2.1 as published by the Free Software Foundation
10 * (the "LGPL") or, at your option, under the terms of the Mozilla
11 * Public License Version 1.1 (the "MPL"). If you do not alter this
12 * notice, a recipient may use your version of this file under either
13 * the MPL or the LGPL.
14 *
15 * You should have received a copy of the LGPL along with this library
16 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 * You should have received a copy of the MPL along with this library
19 * in the file COPYING-MPL-1.1
20 *
21 * The contents of this file are subject to the Mozilla Public License
22 * Version 1.1 (the "License"); you may not use this file except in
23 * compliance with the License. You may obtain a copy of the License at
24 * http://www.mozilla.org/MPL/
25 *
26 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
27 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
28 * the specific language governing rights and limitations.
29 */
30
31 #include <algorithm>
32 #include <optional>
33 #include <2geom/line.h>
34 #include <2geom/math-utils.h>
35
36 namespace Geom
37 {
38
39 /**
40 * @class Line
41 * @brief Infinite line on a plane.
42 *
43 * A line is specified as two points through which it passes. Lines can be interpreted as functions
44 * \f$ f: (-\infty, \infty) \to \mathbb{R}^2\f$. Zero corresponds to the first (origin) point,
45 * one corresponds to the second (final) point. All other points are computed as a linear
46 * interpolation between those two: \f$p = (1-t) a + t b\f$. Many such functions have the same
47 * image and therefore represent the same lines; for example, adding \f$b-a\f$ to both points
48 * yields the same line.
49 *
50 * 2Geom can represent the same line in many ways by design: using a different representation
51 * would lead to precision loss. For example, a line from (1e30, 1e30) to (10,0) would actually
52 * evaluate to (0,0) at time 1 if it was stored as origin and normalized versor,
53 * or origin and angle.
54 *
55 * @ingroup Primitives
56 */
57
58 /** @brief Set the line by solving the line equation.
59 * A line is a set of points that satisfies the line equation
60 * \f$Ax + By + C = 0\f$. This function changes the line so that its points
61 * satisfy the line equation with the given coefficients. */
62 25 void Line::setCoefficients (Coord a, Coord b, Coord c)
63 {
64 // degenerate case
65
3/4
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
25 if (a == 0 && b == 0) {
66 if (c != 0) {
67 THROW_LOGICALERROR("the passed coefficients give the empty set");
68 }
69 _initial = Point(0,0);
70 _final = Point(0,0);
71 return;
72 }
73
74 // The way final / initial points are set based on coefficients is somewhat unusual.
75 // This is done to make sure that calling coefficients() will give back
76 // (almost) the same values.
77
78 // vertical case
79
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 22 times.
25 if (a == 0) {
80 // b must be nonzero
81 3 _initial = Point(-b/2, -c / b);
82 3 _final = _initial;
83 3 _final[X] = b/2;
84 3 return;
85 }
86
87 // horizontal case
88
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 10 times.
22 if (b == 0) {
89 12 _initial = Point(-c / a, a/2);
90 12 _final = _initial;
91 12 _final[Y] = -a/2;
92 12 return;
93 }
94
95 // This gives reasonable results regardless of the magnitudes of a, b and c.
96 10 _initial = Point(-b/2,a/2);
97 10 _final = Point(b/2,-a/2);
98
99 10 Point offset(-c/(2*a), -c/(2*b));
100
101 10 _initial += offset;
102 10 _final += offset;
103 }
104
105 26874 void Line::coefficients(Coord &a, Coord &b, Coord &c) const
106 {
107
1/2
✓ Branch 3 taken 26874 times.
✗ Branch 4 not taken.
26874 Point v = vector().cw();
108 26874 a = v[X];
109 26874 b = v[Y];
110 26874 c = cross(_initial, _final);
111 26874 }
112
113 /** @brief Get the implicit line equation coefficients.
114 * Note that conversion to implicit form always causes loss of
115 * precision when dealing with lines that start far from the origin
116 * and end very close to it. It is recommended to normalize the line
117 * before converting it to implicit form.
118 * @return Vector with three values corresponding to the A, B and C
119 * coefficients of the line equation for this line. */
120 std::vector<Coord> Line::coefficients() const
121 {
122 std::vector<Coord> c(3);
123 coefficients(c[0], c[1], c[2]);
124 return c;
125 }
126
127 /** @brief Find intersection with an axis-aligned line.
128 * @param v Coordinate of the axis-aligned line
129 * @param d Which axis the coordinate is on. X means a vertical line, Y means a horizontal line.
130 * @return Time values at which this line intersects the query line. */
131 std::vector<Coord> Line::roots(Coord v, Dim2 d) const {
132 std::vector<Coord> result;
133 Coord r = root(v, d);
134 if (std::isfinite(r)) {
135 result.push_back(r);
136 }
137 return result;
138 }
139
140 Coord Line::root(Coord v, Dim2 d) const
141 {
142 assert(d == X || d == Y);
143 Point vs = vector();
144 if (vs[d] != 0) {
145 return (v - _initial[d]) / vs[d];
146 } else {
147 return nan("");
148 }
149 }
150
151 std::optional<LineSegment> Line::clip(Rect const &r) const
152 {
153 Point v = vector();
154 // handle horizontal and vertical lines first,
155 // since the root-based code below will break for them
156 for (unsigned i = 0; i < 2; ++i) {
157 Dim2 d = (Dim2) i;
158 Dim2 o = other_dimension(d);
159 if (v[d] != 0) continue;
160 if (r[d].contains(_initial[d])) {
161 Point a, b;
162 a[o] = r[o].min();
163 b[o] = r[o].max();
164 a[d] = b[d] = _initial[d];
165 if (v[o] > 0) {
166 return LineSegment(a, b);
167 } else {
168 return LineSegment(b, a);
169 }
170 } else {
171 return std::nullopt;
172 }
173 }
174
175 Interval xpart(root(r[X].min(), X), root(r[X].max(), X));
176 Interval ypart(root(r[Y].min(), Y), root(r[Y].max(), Y));
177 if (!xpart.isFinite() || !ypart.isFinite()) {
178 return std::nullopt;
179 }
180
181 OptInterval common = xpart & ypart;
182 if (common) {
183 Point p1 = pointAt(common->min()), p2 = pointAt(common->max());
184 LineSegment result(r.clamp(p1), r.clamp(p2));
185 return result;
186 } else {
187 return std::nullopt;
188 }
189
190 /* old implementation using coefficients:
191
192 if (fabs(b) > fabs(a)) {
193 p0 = Point(r[X].min(), (-c - a*r[X].min())/b);
194 if (p0[Y] < r[Y].min())
195 p0 = Point((-c - b*r[Y].min())/a, r[Y].min());
196 if (p0[Y] > r[Y].max())
197 p0 = Point((-c - b*r[Y].max())/a, r[Y].max());
198 p1 = Point(r[X].max(), (-c - a*r[X].max())/b);
199 if (p1[Y] < r[Y].min())
200 p1 = Point((-c - b*r[Y].min())/a, r[Y].min());
201 if (p1[Y] > r[Y].max())
202 p1 = Point((-c - b*r[Y].max())/a, r[Y].max());
203 } else {
204 p0 = Point((-c - b*r[Y].min())/a, r[Y].min());
205 if (p0[X] < r[X].min())
206 p0 = Point(r[X].min(), (-c - a*r[X].min())/b);
207 if (p0[X] > r[X].max())
208 p0 = Point(r[X].max(), (-c - a*r[X].max())/b);
209 p1 = Point((-c - b*r[Y].max())/a, r[Y].max());
210 if (p1[X] < r[X].min())
211 p1 = Point(r[X].min(), (-c - a*r[X].min())/b);
212 if (p1[X] > r[X].max())
213 p1 = Point(r[X].max(), (-c - a*r[X].max())/b);
214 }
215 return LineSegment(p0, p1); */
216 }
217
218 /** @brief Get a time value corresponding to a point.
219 * @param p Point on the line. If the point is not on the line,
220 * the returned value will be meaningless.
221 * @return Time value t such that \f$f(t) = p\f$.
222 * @see timeAtProjection */
223 40070 Coord Line::timeAt(Point const &p) const
224 {
225
1/2
✓ Branch 2 taken 40070 times.
✗ Branch 3 not taken.
40070 Point v = vector();
226 // degenerate case
227
4/6
✓ Branch 1 taken 23 times.
✓ Branch 2 taken 40047 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 23 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 40070 times.
40070 if (v[X] == 0 && v[Y] == 0) {
228 return 0;
229 }
230
231 // use the coordinate that will give better precision
232
2/2
✓ Branch 2 taken 38257 times.
✓ Branch 3 taken 1813 times.
40070 if (fabs(v[X]) > fabs(v[Y])) {
233 38257 return (p[X] - _initial[X]) / v[X];
234 } else {
235 1813 return (p[Y] - _initial[Y]) / v[Y];
236 }
237 }
238
239 /** @brief Create a transformation that maps one line to another.
240 * This will return a transformation \f$A\f$ such that
241 * \f$L_1(t) \cdot A = L_2(t)\f$, where \f$L_1\f$ is this line
242 * and \f$L_2\f$ is the line passed as the parameter. The returned
243 * transformation will preserve angles. */
244 Affine Line::transformTo(Line const &other) const
245 {
246 Affine result = Translate(-_initial);
247 result *= Rotate(angle_between(vector(), other.vector()));
248 result *= Scale(other.vector().length() / vector().length());
249 result *= Translate(other._initial);
250 return result;
251 }
252
253 20016 std::vector<ShapeIntersection> Line::intersect(Line const &other) const
254 {
255 20016 std::vector<ShapeIntersection> result;
256
257
1/2
✓ Branch 2 taken 20016 times.
✗ Branch 3 not taken.
20016 Point v1 = vector();
258
1/2
✓ Branch 2 taken 20016 times.
✗ Branch 3 not taken.
20016 Point v2 = other.vector();
259 20016 Coord cp = cross(v1, v2);
260
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 20005 times.
20016 if (cp == 0) return result;
261
262 20005 Point odiff = other.initialPoint() - initialPoint();
263 20005 Coord t1 = cross(odiff, v2) / cp;
264 20005 Coord t2 = cross(odiff, v1) / cp;
265
1/2
✓ Branch 1 taken 20005 times.
✗ Branch 2 not taken.
20005 result.emplace_back(*this, other, t1, t2);
266 20005 return result;
267 }
268
269 std::vector<ShapeIntersection> Line::intersect(Ray const &r) const
270 {
271 Line other(r);
272 std::vector<ShapeIntersection> result = intersect(other);
273 filter_ray_intersections(result, false, true);
274 return result;
275 }
276
277 20004 std::vector<ShapeIntersection> Line::intersect(LineSegment const &ls) const
278 {
279
1/2
✓ Branch 2 taken 20004 times.
✗ Branch 3 not taken.
20004 Line other(ls);
280
1/2
✓ Branch 1 taken 20004 times.
✗ Branch 2 not taken.
20004 std::vector<ShapeIntersection> result = intersect(other);
281
1/2
✓ Branch 1 taken 20004 times.
✗ Branch 2 not taken.
20004 filter_line_segment_intersections(result, false, true);
282 40008 return result;
283 }
284
285
286
287 20004 void filter_line_segment_intersections(std::vector<ShapeIntersection> &xs, bool a, bool b)
288 {
289 20004 Interval unit(0, 1);
290 20004 std::vector<ShapeIntersection>::reverse_iterator i = xs.rbegin(), last = xs.rend();
291
2/2
✓ Branch 1 taken 20004 times.
✓ Branch 2 taken 20004 times.
40008 while (i != last) {
292
7/14
✗ Branch 0 not taken.
✓ Branch 1 taken 20004 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 20004 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 20004 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 2 times.
✓ Branch 15 taken 20002 times.
✓ Branch 16 taken 2 times.
✓ Branch 17 taken 20002 times.
20004 if ((a && !unit.contains(i->first)) || (b && !unit.contains(i->second))) {
293
1/2
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
2 xs.erase((++i).base());
294 } else {
295 20002 ++i;
296 }
297 }
298 20004 }
299
300 void filter_ray_intersections(std::vector<ShapeIntersection> &xs, bool a, bool b)
301 {
302 Interval unit(0, 1);
303 std::vector<ShapeIntersection>::reverse_iterator i = xs.rbegin(), last = xs.rend();
304 while (i != last) {
305 if ((a && i->first < 0) || (b && i->second < 0)) {
306 xs.erase((++i).base());
307 } else {
308 ++i;
309 }
310 }
311 }
312
313 namespace detail
314 {
315
316 inline
317 OptCrossing intersection_impl(Point const &v1, Point const &o1,
318 Point const &v2, Point const &o2)
319 {
320 Coord cp = cross(v1, v2);
321 if (cp == 0) return OptCrossing();
322
323 Point odiff = o2 - o1;
324
325 Crossing c;
326 c.ta = cross(odiff, v2) / cp;
327 c.tb = cross(odiff, v1) / cp;
328 return c;
329 }
330
331
332 OptCrossing intersection_impl(Ray const& r1, Line const& l2, unsigned int i)
333 {
334 using std::swap;
335
336 OptCrossing crossing =
337 intersection_impl(r1.vector(), r1.origin(),
338 l2.vector(), l2.origin() );
339
340 if (crossing) {
341 if (crossing->ta < 0) {
342 return OptCrossing();
343 } else {
344 if (i != 0) {
345 swap(crossing->ta, crossing->tb);
346 }
347 return crossing;
348 }
349 }
350 if (are_near(r1.origin(), l2)) {
351 THROW_INFINITESOLUTIONS();
352 } else {
353 return OptCrossing();
354 }
355 }
356
357
358 OptCrossing intersection_impl( LineSegment const& ls1,
359 Line const& l2,
360 unsigned int i )
361 {
362 using std::swap;
363
364 OptCrossing crossing =
365 intersection_impl(ls1.finalPoint() - ls1.initialPoint(),
366 ls1.initialPoint(),
367 l2.vector(),
368 l2.origin() );
369
370 if (crossing) {
371 if ( crossing->getTime(0) < 0
372 || crossing->getTime(0) > 1 )
373 {
374 return OptCrossing();
375 } else {
376 if (i != 0) {
377 swap((*crossing).ta, (*crossing).tb);
378 }
379 return crossing;
380 }
381 }
382 if (are_near(ls1.initialPoint(), l2)) {
383 THROW_INFINITESOLUTIONS();
384 } else {
385 return OptCrossing();
386 }
387 }
388
389
390 OptCrossing intersection_impl( LineSegment const& ls1,
391 Ray const& r2,
392 unsigned int i )
393 {
394 using std::swap;
395
396 Point direction = ls1.finalPoint() - ls1.initialPoint();
397 OptCrossing crossing =
398 intersection_impl( direction,
399 ls1.initialPoint(),
400 r2.vector(),
401 r2.origin() );
402
403 if (crossing) {
404 if ( (crossing->getTime(0) < 0)
405 || (crossing->getTime(0) > 1)
406 || (crossing->getTime(1) < 0) )
407 {
408 return OptCrossing();
409 } else {
410 if (i != 0) {
411 swap(crossing->ta, crossing->tb);
412 }
413 return crossing;
414 }
415 }
416
417 if ( are_near(r2.origin(), ls1) ) {
418 bool eqvs = (dot(direction, r2.vector()) > 0);
419 if ( are_near(ls1.initialPoint(), r2.origin()) && !eqvs) {
420 crossing->ta = crossing->tb = 0;
421 return crossing;
422 } else if ( are_near(ls1.finalPoint(), r2.origin()) && eqvs) {
423 if (i == 0) {
424 crossing->ta = 1;
425 crossing->tb = 0;
426 } else {
427 crossing->ta = 0;
428 crossing->tb = 1;
429 }
430 return crossing;
431 } else {
432 THROW_INFINITESOLUTIONS();
433 }
434 } else if ( are_near(ls1.initialPoint(), r2) ) {
435 THROW_INFINITESOLUTIONS();
436 } else {
437 OptCrossing no_crossing;
438 return no_crossing;
439 }
440 }
441
442 } // end namespace detail
443
444
445
446 OptCrossing intersection(Line const& l1, Line const& l2)
447 {
448 OptCrossing c = detail::intersection_impl(
449 l1.vector(), l1.origin(),
450 l2.vector(), l2.origin());
451
452 if (!c && distance(l1.origin(), l2) == 0) {
453 THROW_INFINITESOLUTIONS();
454 }
455 return c;
456 }
457
458 OptCrossing intersection(Ray const& r1, Ray const& r2)
459 {
460 OptCrossing crossing =
461 detail::intersection_impl( r1.vector(), r1.origin(),
462 r2.vector(), r2.origin() );
463
464 if (crossing)
465 {
466 if ( crossing->ta < 0
467 || crossing->tb < 0 )
468 {
469 OptCrossing no_crossing;
470 return no_crossing;
471 }
472 else
473 {
474 return crossing;
475 }
476 }
477
478 if ( are_near(r1.origin(), r2) || are_near(r2.origin(), r1) )
479 {
480 if ( are_near(r1.origin(), r2.origin())
481 && !are_near(r1.vector(), r2.vector()) )
482 {
483 crossing->ta = crossing->tb = 0;
484 return crossing;
485 }
486 else
487 {
488 THROW_INFINITESOLUTIONS();
489 }
490 }
491 else
492 {
493 OptCrossing no_crossing;
494 return no_crossing;
495 }
496 }
497
498
499 OptCrossing intersection( LineSegment const& ls1, LineSegment const& ls2 )
500 {
501 Point direction1 = ls1.finalPoint() - ls1.initialPoint();
502 Point direction2 = ls2.finalPoint() - ls2.initialPoint();
503 OptCrossing crossing =
504 detail::intersection_impl( direction1,
505 ls1.initialPoint(),
506 direction2,
507 ls2.initialPoint() );
508
509 if (crossing)
510 {
511 if ( crossing->getTime(0) < 0
512 || crossing->getTime(0) > 1
513 || crossing->getTime(1) < 0
514 || crossing->getTime(1) > 1 )
515 {
516 OptCrossing no_crossing;
517 return no_crossing;
518 }
519 else
520 {
521 return crossing;
522 }
523 }
524
525 bool eqvs = (dot(direction1, direction2) > 0);
526 if ( are_near(ls2.initialPoint(), ls1) )
527 {
528 if ( are_near(ls1.initialPoint(), ls2.initialPoint()) && !eqvs )
529 {
530 crossing->ta = crossing->tb = 0;
531 return crossing;
532 }
533 else if ( are_near(ls1.finalPoint(), ls2.initialPoint()) && eqvs )
534 {
535 crossing->ta = 1;
536 crossing->tb = 0;
537 return crossing;
538 }
539 else
540 {
541 THROW_INFINITESOLUTIONS();
542 }
543 }
544 else if ( are_near(ls2.finalPoint(), ls1) )
545 {
546 if ( are_near(ls1.finalPoint(), ls2.finalPoint()) && !eqvs )
547 {
548 crossing->ta = crossing->tb = 1;
549 return crossing;
550 }
551 else if ( are_near(ls1.initialPoint(), ls2.finalPoint()) && eqvs )
552 {
553 crossing->ta = 0;
554 crossing->tb = 1;
555 return crossing;
556 }
557 else
558 {
559 THROW_INFINITESOLUTIONS();
560 }
561 }
562 else
563 {
564 OptCrossing no_crossing;
565 return no_crossing;
566 }
567 }
568
569 Line make_angle_bisector_line(Line const& l1, Line const& l2)
570 {
571 OptCrossing crossing;
572 try
573 {
574 crossing = intersection(l1, l2);
575 }
576 catch(InfiniteSolutions const &e)
577 {
578 return l1;
579 }
580 if (!crossing)
581 {
582 THROW_RANGEERROR("passed lines are parallel");
583 }
584 Point O = l1.pointAt(crossing->ta);
585 Point A = l1.pointAt(crossing->ta + 1);
586 double angle = angle_between(l1.vector(), l2.vector());
587 Point B = (angle > 0) ? l2.pointAt(crossing->tb + 1)
588 : l2.pointAt(crossing->tb - 1);
589
590 return make_angle_bisector_line(A, O, B);
591 }
592
593
594
595
596 } // end namespace Geom
597
598
599
600 /*
601 Local Variables:
602 mode:c++
603 c-file-style:"stroustrup"
604 c-file-offsets:((innamespace . 0)(substatement-open . 0))
605 indent-tabs-mode:nil
606 c-brace-offset:0
607 fill-column:99
608 End:
609 vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 :
610 */
611