Inkscape
Vector Graphics Editor
Loading...
Searching...
No Matches
line.cpp
Go to the documentation of this file.
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
36namespace Geom
37{
38
63{
64 // degenerate case
65 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 if (a == 0) {
80 // b must be nonzero
81 _initial = Point(-b/2, -c / b);
83 _final[X] = b/2;
84 return;
85 }
86
87 // horizontal case
88 if (b == 0) {
89 _initial = Point(-c / a, a/2);
91 _final[Y] = -a/2;
92 return;
93 }
94
95 // This gives reasonable results regardless of the magnitudes of a, b and c.
96 _initial = Point(-b/2,a/2);
97 _final = Point(b/2,-a/2);
98
99 Point offset(-c/(2*a), -c/(2*b));
100
101 _initial += offset;
102 _final += offset;
103}
104
106{
107 Point v = vector().cw();
108 a = v[X];
109 b = v[Y];
111}
112
120std::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
131std::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
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
151std::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
223Coord Line::timeAt(Point const &p) const
224{
225 Point v = vector();
226 // degenerate case
227 if (v[X] == 0 && v[Y] == 0) {
228 return 0;
229 }
230
231 // use the coordinate that will give better precision
232 if (fabs(v[X]) > fabs(v[Y])) {
233 return (p[X] - _initial[X]) / v[X];
234 } else {
235 return (p[Y] - _initial[Y]) / v[Y];
236 }
237}
238
244Affine Line::transformTo(Line const &other) const
245{
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
253std::vector<ShapeIntersection> Line::intersect(Line const &other) const
254{
255 std::vector<ShapeIntersection> result;
256
257 Point v1 = vector();
258 Point v2 = other.vector();
259 Coord cp = cross(v1, v2);
260 if (cp == 0) return result;
261
262 Point odiff = other.initialPoint() - initialPoint();
263 Coord t1 = cross(odiff, v2) / cp;
264 Coord t2 = cross(odiff, v1) / cp;
265 result.emplace_back(*this, other, t1, t2);
266 return result;
267}
268
269std::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
277std::vector<ShapeIntersection> Line::intersect(LineSegment const &ls) const
278{
279 Line other(ls);
280 std::vector<ShapeIntersection> result = intersect(other);
282 return result;
283}
284
285
286
287void filter_line_segment_intersections(std::vector<ShapeIntersection> &xs, bool a, bool b)
288{
289 Interval unit(0, 1);
290 std::vector<ShapeIntersection>::reverse_iterator i = xs.rbegin(), last = xs.rend();
291 while (i != last) {
292 if ((a && !unit.contains(i->first)) || (b && !unit.contains(i->second))) {
293 xs.erase((++i).base());
294 } else {
295 ++i;
296 }
297 }
298}
299
300void 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
313namespace detail
314{
315
316inline
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
332OptCrossing intersection_impl(Ray const& r1, Line const& l2, unsigned int i)
333{
334 using std::swap;
335
336 OptCrossing crossing =
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
359 Line const& l2,
360 unsigned int i )
361{
362 using std::swap;
363
364 OptCrossing crossing =
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
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
446OptCrossing intersection(Line const& l1, Line const& l2)
447{
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
458OptCrossing intersection(Ray const& r1, Ray const& r2)
459{
460 OptCrossing crossing =
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
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
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*/
3x3 matrix representing an affine transformation.
Definition affine.h:70
Point finalPoint() const override
Retrieve the end of the curve.
Point initialPoint() const override
Retrieve the start of the curve.
constexpr bool contains(C val) const
Check whether the interval includes this number.
CPoint clamp(CPoint const &p) const
Clamp point to the rectangle.
CPoint min() const
Get the corner of the rectangle with smallest coordinate values.
CPoint max() const
Get the corner of the rectangle with largest coordinate values.
Range of real numbers that is never empty.
Definition interval.h:59
bool isFinite() const
Definition interval.h:94
Infinite line on a plane.
Definition line.h:53
void setCoefficients(double a, double b, double c)
Set the coefficients of the line equation.
Definition line.cpp:62
std::vector< double > coefficients() const
Get the coefficients of the line equation as a vector.
Definition line.cpp:120
Coord timeAt(Point const &p) const
Get a time value corresponding to a point.
Definition line.cpp:223
Coord root(Coord v, Dim2 d) const
Definition line.cpp:140
Point origin() const
Get the line's origin point.
Definition line.h:128
Point vector() const
Get the line's raw direction vector.
Definition line.h:132
std::vector< Coord > roots(Coord v, Dim2 d) const
Find intersection with an axis-aligned line.
Definition line.cpp:131
std::optional< LineSegment > clip(Rect const &r) const
Return the portion of the line that is inside the given rectangle.
Definition line.cpp:151
Affine transformTo(Line const &other) const
Create a transformation that maps one line to another.
Definition line.cpp:244
Point _final
Definition line.h:56
Point initialPoint() const
Definition line.h:225
std::vector< ShapeIntersection > intersect(Line const &other) const
Definition line.cpp:253
Point pointAt(Coord t) const
Definition line.h:231
Point _initial
Definition line.h:55
Range of real numbers that can be empty.
Definition interval.h:199
Two-dimensional point that doubles as a vector.
Definition point.h:66
Coord length() const
Compute the distance from origin.
Definition point.h:118
constexpr Point cw() const
Return a point like this point but rotated +90 degrees.
Definition point.h:137
Straight ray from a specific point to infinity.
Definition ray.h:53
Point origin() const
Definition ray.h:68
Point vector() const
Definition ray.h:69
Axis aligned, non-empty rectangle.
Definition rect.h:92
Rotation around the origin.
Definition transforms.h:187
Scaling from the origin.
Definition transforms.h:150
Translation by a vector.
Definition transforms.h:115
RootCluster root
Css & result
double c[8][4]
BezierCurveN< 1 > LineSegment
Line segment.
Dim2
2D axis enumeration (X or Y).
Definition coord.h:48
constexpr Dim2 other_dimension(Dim2 d)
Get the other (perpendicular) dimension.
Definition coord.h:52
double Coord
Floating point type used to store coordinates.
Definition coord.h:76
@ Y
Definition coord.h:48
@ X
Definition coord.h:48
Infinite straight line.
Low level math functions and compatibility wrappers.
double offset
OptCrossing intersection_impl(Ray const &r1, Line const &l2, unsigned int i)
Definition line.cpp:332
Various utility functions.
Definition affine.h:22
Coord length(LineSegment const &seg)
void filter_ray_intersections(std::vector< ShapeIntersection > &xs, bool a=false, bool b=true)
Definition line.cpp:300
MultiDegree< n > max(MultiDegree< n > const &p, MultiDegree< n > const &q)
Returns the maximal degree appearing in the two arguments for each variables.
Definition sbasisN.h:158
Angle distance(Angle const &a, Angle const &b)
Definition angle.h:163
std::optional< Crossing > OptCrossing
Definition crossing.h:64
OptCrossing intersection(Ray const &r1, Line const &l2)
Definition line.h:545
double angle_between(Line const &l1, Line const &l2)
Definition line.h:456
bool contains(Path const &p, Point const &i, bool evenodd=true)
void filter_line_segment_intersections(std::vector< ShapeIntersection > &xs, bool a=false, bool b=true)
Removes intersections outside of the unit interval.
Definition line.cpp:287
Piecewise< SBasis > cross(Piecewise< D2< SBasis > > const &a, Piecewise< D2< SBasis > > const &b)
Piecewise< SBasis > min(SBasis const &f, SBasis const &g)
Return the more negative of the two functions pointwise.
Line make_angle_bisector_line(Point const &A, Point const &O, Point const &B)
Definition line.h:504
T dot(D2< T > const &a, D2< T > const &b)
Definition d2.h:355
bool are_near(Affine const &a1, Affine const &a2, Coord eps=EPSILON)