1 package org.opentrafficsim.core.geometry;
2
3 import java.awt.geom.Line2D;
4 import java.awt.geom.Path2D;
5 import java.awt.geom.PathIterator;
6 import java.awt.geom.Point2D;
7 import java.io.Serializable;
8 import java.util.ArrayList;
9 import java.util.Arrays;
10 import java.util.List;
11
12 import org.djunits.unit.DirectionUnit;
13 import org.djunits.value.vdouble.scalar.Direction;
14 import org.djunits.value.vdouble.scalar.Length;
15 import org.djutils.draw.bounds.Bounds2d;
16 import org.djutils.draw.line.PolyLine2d;
17 import org.djutils.draw.line.Ray2d;
18 import org.djutils.draw.point.OrientedPoint2d;
19 import org.djutils.draw.point.Point2d;
20 import org.djutils.exceptions.Throw;
21 import org.djutils.exceptions.Try;
22
23 import nl.tudelft.simulation.dsol.animation.Locatable;
24
25 /**
26 * Line with underlying PolyLine2d, a cached length indexed line, a cached length, and a cached centroid (all calculated on
27 * first use). This class supports fractional projection.
28 * <p>
29 * Copyright (c) 2013-2024 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
30 * BSD-style license. See <a href="https://opentrafficsim.org/docs/license.html">OpenTrafficSim License</a>.
31 * </p>
32 * @author <a href="https://github.com/averbraeck">Alexander Verbraeck</a>
33 * @author <a href="https://tudelft.nl/staff/p.knoppers-1">Peter Knoppers</a>
34 * @author <a href="https://www.citg.tudelft.nl">Guus Tamminga</a>
35 * @author <a href="https://github.com/wjschakel">Wouter Schakel</a>
36 */
37 public class OtsLine2d implements Locatable, Serializable
38 {
39 /** */
40 private static final long serialVersionUID = 20150722L;
41
42 /** The 2d line. */
43 private PolyLine2d line2d;
44
45 /** The cumulative length of the line at point 'i'. */
46 private double[] lengthIndexedLine = null;
47
48 /** The cached length; will be calculated at time of construction. */
49 private Length length;
50
51 /** The cached centroid; will be calculated when needed for the first time. */
52 private Point2d centroid = null;
53
54 /** The cached bounds; will be calculated when needed for the first time. */
55 private Bounds2d bounds = null;
56
57 /** The cached helper points for fractional projection; will be calculated when needed for the first time. */
58 private Point2d[] fractionalHelperCenters = null;
59
60 /** The cached helper directions for fractional projection; will be calculated when needed for the first time. */
61 private Point2D.Double[] fractionalHelperDirections = null;
62
63 /** Intersection of unit offset lines of first two segments. */
64 private Point2d firstOffsetIntersection;
65
66 /** Intersection of unit offset lines of last two segments. */
67 private Point2d lastOffsetIntersection;
68
69 /** Precision for fractional projection algorithm. */
70 private static final double FRAC_PROJ_PRECISION = 2e-5 /* PK too fine 1e-6 */;
71
72 /** Radius at each vertex. */
73 private Length[] vertexRadii;
74
75 /**
76 * Construct a new OtsLine2d.
77 * @param points Point2d...; the array of points to construct this OtsLine2d from.
78 */
79 public OtsLine2d(final Point2d... points)
80 {
81 this(new PolyLine2d(points));
82 }
83
84 /**
85 * Creates a new OtsLine2d based on 2d information. Elevation will be 0.
86 * @param line2d PolyLine2d; 2d information.
87 */
88 public OtsLine2d(final PolyLine2d line2d)
89 {
90 init(line2d);
91 }
92
93 /**
94 * Construct a new OtsLine2d, and immediately make the length-indexed line.
95 * @param line2d PolyLine2d; the 2d line.
96 */
97 private void init(final PolyLine2d line2d)
98 {
99 this.lengthIndexedLine = new double[line2d.size()];
100 this.lengthIndexedLine[0] = 0.0;
101 for (int i = 1; i < line2d.size(); i++)
102 {
103 this.lengthIndexedLine[i] = this.lengthIndexedLine[i - 1] + line2d.get(i - 1).distance(line2d.get(i));
104 }
105 this.line2d = line2d;
106 this.length = Length.instantiateSI(this.lengthIndexedLine[this.lengthIndexedLine.length - 1]);
107 }
108
109 /**
110 * Construct parallel line.<br>
111 * @param offset double; offset distance from the reference line; positive is LEFT, negative is RIGHT
112 * @return OtsLine2d; the line that has the specified offset from the reference line
113 */
114 public final OtsLine2d offsetLine(final double offset)
115 {
116 return new OtsLine2d(this.line2d.offsetLine(offset));
117 }
118
119 /**
120 * Clean up a list of points that describe a polyLine by removing points that lie within epsilon distance of a more
121 * straightened version of the line. <br>
122 * @param epsilon double; maximal deviation
123 * @param useHorizontalDistance boolean; if true; the horizontal distance is used; if false; the 3D distance is used
124 * @return OtsLine2d; a new OtsLine2d containing all the remaining points
125 */
126 @Deprecated
127 public final OtsLine2d noiseFilterRamerDouglasPeucker(final double epsilon, final boolean useHorizontalDistance)
128 {
129 // Apply the Ramer-Douglas-Peucker algorithm to the buffered points.
130 // Adapted from https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm
131 double maxDeviation = 0;
132 int splitIndex = -1;
133 int pointCount = size();
134 // Find the point with largest deviation from the straight line from start point to end point
135 for (int i = 1; i < pointCount - 1; i++)
136 {
137 Point2d point = this.line2d.get(i);
138 Point2d closest = point.closestPointOnLine(this.line2d.get(0), this.line2d.get(pointCount - 1));
139 double deviation = useHorizontalDistance ? closest.distance(point) : closest.distance(point);
140 if (deviation > maxDeviation)
141 {
142 splitIndex = i;
143 maxDeviation = deviation;
144 }
145 }
146 if (maxDeviation <= epsilon)
147 {
148 // All intermediate points can be dropped. Return a new list containing only the first and last point.
149 return new OtsLine2d(this.line2d.get(0), this.line2d.get(pointCount - 1));
150 }
151 // The largest deviation is larger than epsilon.
152 // Split the polyLine at the point with the maximum deviation. Process each sub list recursively and concatenate the
153 // results
154 List<Point2d> points = this.line2d.getPointList();
155 OtsLine2d first = new OtsLine2d(points.subList(0, splitIndex + 1).toArray(new Point2d[splitIndex + 1]))
156 .noiseFilterRamerDouglasPeucker(epsilon, useHorizontalDistance);
157 OtsLine2d second = new OtsLine2d(
158 points.subList(splitIndex, this.line2d.size()).toArray(new Point2d[this.line2d.size() - splitIndex]))
159 .noiseFilterRamerDouglasPeucker(epsilon, useHorizontalDistance);
160 return concatenate(epsilon, first, second);
161 }
162
163 /**
164 * Returns a 2d representation of this line.
165 * @return PolyLine2d; Returns a 2d representation of this line.
166 */
167 public PolyLine2d getLine2d()
168 {
169 return this.line2d;
170 }
171
172 /**
173 * Create a line at linearly varying offset from this line. The offset may change linearly from its initial value at the
174 * start of the reference line to its final offset value at the end of the reference line.
175 * @param offsetAtStart double; offset at the start of the reference line (positive value is Left, negative value is Right)
176 * @param offsetAtEnd double; offset at the end of the reference line (positive value is Left, negative value is Right)
177 * @return OtsLine2d; the OtsLine2d of the line at linearly changing offset of the reference line
178 */
179 public final OtsLine2d offsetLine(final double offsetAtStart, final double offsetAtEnd)
180 {
181 return new OtsLine2d(this.line2d.offsetLine(offsetAtStart, offsetAtEnd));
182 }
183
184 /**
185 * Create a line at linearly varying offset from this line. The offset may change linearly from its initial value at the
186 * start of the reference line via a number of intermediate offsets at intermediate positions to its final offset value at
187 * the end of the reference line.
188 * @param relativeFractions double[]; positional fractions for which the offsets have to be generated
189 * @param offsets double[]; offsets at the relative positions (positive value is Left, negative value is Right)
190 * @return Geometry; the Geometry of the line at linearly changing offset of the reference line
191 * @throws OtsGeometryException when this method fails to create the offset line
192 */
193 public final OtsLine2d offsetLine(final double[] relativeFractions, final double[] offsets) throws OtsGeometryException
194 {
195 return new OtsLine2d(OtsGeometryUtil.offsetLine(this.line2d, relativeFractions, offsets));
196 }
197
198 /**
199 * Concatenate several OtsLine2d instances.
200 * @param lines OtsLine2d...; OtsLine2d... one or more OtsLine2d. The last point of the first
201 * <strong>must</strong> match the first of the second, etc.
202 * @return OtsLine2d
203 */
204 public static OtsLine2d concatenate(final OtsLine2d... lines)
205 {
206 return concatenate(0.0, lines);
207 }
208
209 /**
210 * Concatenate two OtsLine2d instances. This method is separate for efficiency reasons.
211 * @param toleranceSI double; the tolerance between the end point of a line and the first point of the next line
212 * @param line1 OtsLine2d; first line
213 * @param line2 OtsLine2d; second line
214 * @return OtsLine2d
215 */
216 public static OtsLine2d concatenate(final double toleranceSI, final OtsLine2d line1, final OtsLine2d line2)
217 {
218 return new OtsLine2d(PolyLine2d.concatenate(toleranceSI, line1.line2d, line2.line2d));
219 }
220
221 /**
222 * Concatenate several OtsLine2d instances.
223 * @param toleranceSI double; the tolerance between the end point of a line and the first point of the next line
224 * @param lines OtsLine2d...; OtsLine2d... one or more OtsLine2d. The last point of the first
225 * <strong>must</strong> match the first of the second, etc.
226 * @return OtsLine2d
227 */
228 public static OtsLine2d concatenate(final double toleranceSI, final OtsLine2d... lines)
229 {
230 List<PolyLine2d> lines2d = new ArrayList<>();
231 for (OtsLine2d line : lines)
232 {
233 lines2d.add(line.line2d);
234 }
235 return new OtsLine2d(PolyLine2d.concatenate(toleranceSI, lines2d.toArray(new PolyLine2d[lines.length])));
236 }
237
238 /**
239 * Construct a new OtsLine2d with all points of this OtsLine2d in reverse order.
240 * @return OtsLine2d; the new OtsLine2d
241 */
242 public final OtsLine2d reverse()
243 {
244 return new OtsLine2d(this.line2d.reverse());
245 }
246
247 /**
248 * Construct a new OtsLine2d covering the indicated fraction of this OtsLine2d.
249 * @param start double; starting point, valid range [0..<cite>end</cite>)
250 * @param end double; ending point, valid range (<cite>start</cite>..1]
251 * @return OtsLine2d; the new OtsLine2d
252 */
253 public final OtsLine2d extractFractional(final double start, final double end)
254 {
255 return extract(start * this.length.si, end * this.length.si);
256 }
257
258 /**
259 * Create a new OtsLine2d that covers a sub-section of this OtsLine2d.
260 * @param start Length; the length along this OtsLine2d where the sub-section starts, valid range [0..<cite>end</cite>)
261 * @param end Length; length along this OtsLine2d where the sub-section ends, valid range
262 * (<cite>start</cite>..<cite>length</cite> (length is the length of this OtsLine2d)
263 * @return OtsLine2d; the selected sub-section
264 */
265 public final OtsLine2d extract(final Length start, final Length end)
266 {
267 return extract(start.si, end.si);
268 }
269
270 /**
271 * Create a new OtsLine2d that covers a sub-section of this OtsLine2d.
272 * @param start double; length along this OtsLine2d where the sub-section starts, valid range [0..<cite>end</cite>)
273 * @param end double; length along this OtsLine2d where the sub-section ends, valid range
274 * (<cite>start</cite>..<cite>length</cite> (length is the length of this OtsLine2d)
275 * @return OtsLine2d; the selected sub-section
276 */
277 public final OtsLine2d extract(final double start, final double end)
278 {
279 return new OtsLine2d(this.line2d.extract(start, end));
280 }
281
282 /**
283 * Create an OtsLine2d, while cleaning repeating successive points.
284 * @param points Point2d...; the coordinates of the line as OtsPoint3d
285 * @return the line
286 * @throws OtsGeometryException when number of points < 2
287 */
288 public static OtsLine2d createAndCleanOtsLine2d(final Point2d... points) throws OtsGeometryException
289 {
290 if (points.length < 2)
291 {
292 throw new OtsGeometryException(
293 "Degenerate OtsLine2d; has " + points.length + " point" + (points.length != 1 ? "s" : ""));
294 }
295 return createAndCleanOtsLine2d(new ArrayList<>(Arrays.asList(points)));
296 }
297
298 /**
299 * Create an OtsLine2d, while cleaning repeating successive points.
300 * @param pointList List<Point2d>; list of the coordinates of the line as OtsPoint3d; any duplicate points in this
301 * list are removed (this method may modify the provided list)
302 * @return OtsLine2d; the line
303 * @throws OtsGeometryException when number of non-equal points < 2
304 */
305 public static OtsLine2d createAndCleanOtsLine2d(final List<Point2d> pointList) throws OtsGeometryException
306 {
307 return new OtsLine2d(new PolyLine2d(true, pointList));
308 }
309
310 /**
311 * Construct a new OtsLine2d from a List<OtsPoint3d>.
312 * @param pointList List<OtsPoint3d>; the list of points to construct this OtsLine2d from.
313 * @throws OtsGeometryException when the provided points do not constitute a valid line (too few points or identical
314 * adjacent points)
315 */
316 public OtsLine2d(final List<Point2d> pointList) throws OtsGeometryException
317 {
318 this(pointList.toArray(new Point2d[pointList.size()]));
319 }
320
321 /**
322 * Construct a new OtsShape (closed shape) from a Path2D. Elevation will be 0.
323 * @param path Path2D; the Path2D to construct this OtsLine2d from.
324 * @throws OtsGeometryException when the provided points do not constitute a valid line (too few points or identical
325 * adjacent points)
326 */
327 public OtsLine2d(final Path2D path) throws OtsGeometryException
328 {
329 List<Point2d> pl = new ArrayList<>();
330 for (PathIterator pi = path.getPathIterator(null); !pi.isDone(); pi.next())
331 {
332 double[] p = new double[6];
333 int segType = pi.currentSegment(p);
334 if (segType == PathIterator.SEG_MOVETO || segType == PathIterator.SEG_LINETO)
335 {
336 pl.add(new Point2d(p[0], p[1]));
337 }
338 else if (segType == PathIterator.SEG_CLOSE)
339 {
340 if (!pl.get(0).equals(pl.get(pl.size() - 1)))
341 {
342 pl.add(new Point2d(pl.get(0).x, pl.get(0).y));
343 }
344 break;
345 }
346 }
347 init(new PolyLine2d(pl.toArray(new Point2d[pl.size() - 1])));
348 }
349
350 /**
351 * Return the number of points in this OtsLine2d. This is the number of points in horizontal plane.
352 * @return the number of points on the line
353 */
354 public final int size()
355 {
356 return this.line2d.size();
357 }
358
359 /**
360 * Return the first point of this OtsLine2d.
361 * @return the first point on the line
362 */
363 public final Point2d getFirst()
364 {
365 return this.line2d.getFirst();
366 }
367
368 /**
369 * Return the last point of this OtsLine2d.
370 * @return the last point on the line
371 */
372 public final Point2d getLast()
373 {
374 return this.line2d.getLast();
375 }
376
377 /**
378 * Return one point of this OtsLine2d.
379 * @param i int; the index of the point to retrieve
380 * @return Point2d; the i-th point of the line
381 * @throws OtsGeometryException when i < 0 or i > the number of points
382 */
383 public final Point2d get(final int i) throws OtsGeometryException
384 {
385 if (i < 0 || i > size() - 1)
386 {
387 throw new OtsGeometryException("OtsLine2d.get(i=" + i + "); i<0 or i>=size(), which is " + size());
388 }
389 return this.line2d.get(i);
390 }
391
392 /**
393 * Return the length of this OtsLine2d in meters. (Assuming that the coordinates of the points constituting this line are
394 * expressed in meters.)
395 * @return the length of the line
396 */
397 public final Length getLength()
398 {
399 return this.length;
400 }
401
402 /**
403 * Return an array of OtsPoint3d that represents this OtsLine2d.
404 * @return the points of this line
405 */
406 public final Point2d[] getPoints()
407 {
408 return this.line2d.getPointList().toArray(new Point2d[this.line2d.size()]);
409 }
410
411 /**
412 * Get the location at a position on the line, with its direction. Position can be below 0 or more than the line length. In
413 * that case, the position will be extrapolated in the direction of the line at its start or end.
414 * @param position Length; the position on the line for which to calculate the point on, before, of after the line
415 * @return a directed point
416 */
417 public final OrientedPoint2d getLocationExtended(final Length position)
418 {
419 return getLocationExtendedSI(position.getSI());
420 }
421
422 /**
423 * Get the location at a position on the line, with its direction. Position can be below 0 or more than the line length. In
424 * that case, the position will be extrapolated in the direction of the line at its start or end.
425 * @param positionSI double; the position on the line for which to calculate the point on, before, of after the line, in SI
426 * units
427 * @return a directed point
428 */
429 public final synchronized OrientedPoint2d getLocationExtendedSI(final double positionSI)
430 {
431 Ray2d ray = this.line2d.getLocationExtended(positionSI);
432 return new OrientedPoint2d(ray.x, ray.y, ray.phi);
433 }
434
435 /**
436 * Get the location at a fraction of the line, with its direction. Fraction should be between 0.0 and 1.0.
437 * @param fraction double; the fraction for which to calculate the point on the line
438 * @return a directed point
439 * @throws OtsGeometryException when fraction less than 0.0 or more than 1.0.
440 */
441 public final OrientedPoint2d getLocationFraction(final double fraction) throws OtsGeometryException
442 {
443 if (fraction < 0.0 || fraction > 1.0)
444 {
445 throw new OtsGeometryException("getLocationFraction for line: fraction < 0.0 or > 1.0. fraction = " + fraction);
446 }
447 return getLocationSI(fraction * this.length.si);
448 }
449
450 /**
451 * Get the location at a fraction of the line, with its direction. Fraction should be between 0.0 and 1.0.
452 * @param fraction double; the fraction for which to calculate the point on the line
453 * @param tolerance double; the delta from 0.0 and 1.0 that will be forgiven
454 * @return a directed point
455 * @throws OtsGeometryException when fraction less than 0.0 or more than 1.0.
456 */
457 public final OrientedPoint2d getLocationFraction(final double fraction, final double tolerance) throws OtsGeometryException
458 {
459 if (fraction < -tolerance || fraction > 1.0 + tolerance)
460 {
461 throw new OtsGeometryException(
462 "getLocationFraction for line: fraction < 0.0 - tolerance or > 1.0 + tolerance; fraction = " + fraction);
463 }
464 double f = fraction < 0 ? 0.0 : fraction > 1.0 ? 1.0 : fraction;
465 return getLocationSI(f * this.length.si);
466 }
467
468 /**
469 * Get the location at a fraction of the line (or outside the line), with its direction.
470 * @param fraction double; the fraction for which to calculate the point on the line
471 * @return a directed point
472 */
473 public final OrientedPoint2d getLocationFractionExtended(final double fraction)
474 {
475 return getLocationExtendedSI(fraction * this.length.si);
476 }
477
478 /**
479 * Get the location at a position on the line, with its direction. Position should be between 0.0 and line length.
480 * @param position Length; the position on the line for which to calculate the point on the line
481 * @return a directed point
482 * @throws OtsGeometryException when position less than 0.0 or more than line length.
483 */
484 public final OrientedPoint2d getLocation(final Length position) throws OtsGeometryException
485 {
486 return getLocationSI(position.getSI());
487 }
488
489 /**
490 * Binary search for a position on the line.
491 * @param pos double; the position to look for.
492 * @return the index below the position; the position is between points[index] and points[index+1]
493 * @throws OtsGeometryException when index could not be found
494 */
495 private int find(final double pos) throws OtsGeometryException
496 {
497 if (pos == 0)
498 {
499 return 0;
500 }
501
502 int lo = 0;
503 int hi = this.lengthIndexedLine.length - 1;
504 while (lo <= hi)
505 {
506 if (hi == lo)
507 {
508 return lo;
509 }
510 int mid = lo + (hi - lo) / 2;
511 if (pos < this.lengthIndexedLine[mid])
512 {
513 hi = mid - 1;
514 }
515 else if (pos > this.lengthIndexedLine[mid + 1])
516 {
517 lo = mid + 1;
518 }
519 else
520 {
521 return mid;
522 }
523 }
524 throw new OtsGeometryException(
525 "Could not find position " + pos + " on line with length indexes: " + Arrays.toString(this.lengthIndexedLine));
526 }
527
528 /**
529 * Get the location at a position on the line, with its direction. Position should be between 0.0 and line length.
530 * @param positionSI double; the position on the line for which to calculate the point on the line
531 * @return a directed point
532 * @throws OtsGeometryException when position less than 0.0 or more than line length.
533 */
534 public final OrientedPoint2d getLocationSI(final double positionSI) throws OtsGeometryException
535 {
536 Ray2d ray = Try.assign(() -> this.line2d.getLocation(positionSI), OtsGeometryException.class, "Position not on line.");
537 return new OrientedPoint2d(ray.x, ray.y, ray.phi);
538 }
539
540 /**
541 * Truncate a line at the given length (less than the length of the line, and larger than zero) and return a new line.
542 * @param lengthSI double; the location where to truncate the line
543 * @return a new OtsLine2d truncated at the exact position where line.getLength() == lengthSI
544 * @throws OtsGeometryException when position less than 0.0 or more than line length.
545 */
546 public final OtsLine2d truncate(final double lengthSI) throws OtsGeometryException
547 {
548 return new OtsLine2d(this.line2d.truncate(lengthSI));
549 }
550
551 /**
552 * Returns the fractional position along this line of the orthogonal projection of point (x, y) on this line. If the point
553 * is not orthogonal to the closest line segment, the nearest point is selected.
554 * @param x double; x-coordinate of point to project
555 * @param y double; y-coordinate of point to project
556 * @return fractional position along this line of the orthogonal projection on this line of a point
557 */
558 public final double projectOrthogonal(final double x, final double y)
559 {
560 Point2d closest = this.line2d.closestPointOnPolyLine(new Point2d(x, y));
561 return this.line2d.projectOrthogonalFractionalExtended(closest);
562 }
563
564 /**
565 * Returns the fractional projection of a point to a line. The projection works by taking slices in space per line segment
566 * as shown below. A point is always projected to the nearest segment, but not necessarily to the closest point on that
567 * segment. The slices in space are analogous to a Voronoi diagram, but for the line segments instead of points. If
568 * fractional projection fails, the orthogonal projection is returned.<br>
569 * <br>
570 * The point 'A' is projected to point 'B' on the 3rd segment of line 'C-D'. The line from 'A' to 'B' extends towards point
571 * 'E', which is the intersection of lines 'E-F' and 'E-G'. Line 'E-F' cuts the first bend of the 3rd segment (at point 'H')
572 * in half, while the line 'E-G' cuts the second bend of the 3rd segment (at point 'I') in half.
573 *
574 * <pre>
575 * ____________________________ G .
576 * . | | . .
577 * . | . . . . helper lines | . .
578 * . | _.._.._ projection line | I. .
579 * . |____________________________| _.'._ . L
580 * F. _.' . '-. .
581 * .. B _.' . '-.
582 * . . _.\ . . D
583 * . . _.' : . .
584 * J . . _.' \ . .
585 * .. . _.' : . M
586 * . . ..-' \ .
587 * . . /H. A .
588 * . . / . .
589 * C _________/ . .
590 * . . . .
591 * K . . . .
592 * . . . .
593 * . . . . N
594 * . . . .
595 * . . . .
596 * . . . .
597 * . . . .
598 * . .E
599 * . .
600 * . .
601 * . .
602 * </pre>
603 *
604 * Fractional projection may fail in three cases.
605 * <ol>
606 * <li>Numerical difficulties at slight bend, orthogonal projection returns the correct point.</li>
607 * <li>Fractional projection is possible only to segments that aren't the nearest segment(s).</li>
608 * <li>Fractional projection is possible for no segment.</li>
609 * </ol>
610 * In the latter two cases the projection is undefined and a orthogonal projection is returned if
611 * {@code orthoFallback = true}, or {@code NaN} if {@code orthoFallback = false}.
612 * @param start Direction; direction in first point
613 * @param end Direction; direction in last point
614 * @param x double; x-coordinate of point to project
615 * @param y double; y-coordinate of point to project
616 * @param fallback FractionalFallback; fallback method for when fractional projection fails
617 * @return fractional position along this line of the fractional projection on that line of a point
618 */
619 public final synchronized double projectFractional(final Direction start, final Direction end, final double x,
620 final double y, final FractionalFallback fallback)
621 {
622
623 // prepare
624 double minDistance = Double.POSITIVE_INFINITY;
625 double minSegmentFraction = 0;
626 int minSegment = -1;
627 Point2d point = new Point2d(x, y);
628
629 // determine helpers (centers and directions)
630 determineFractionalHelpers(start, end);
631
632 // get distance of point to each segment
633 double[] d = new double[size() - 1];
634 double minD = Double.POSITIVE_INFINITY;
635 for (int i = 0; i < size() - 1; i++)
636 {
637 d[i] = Line2D.ptSegDist(this.line2d.get(i).x, this.line2d.get(i).y, this.line2d.get(i + 1).x,
638 this.line2d.get(i + 1).y, x, y);
639 minD = d[i] < minD ? d[i] : minD;
640 }
641
642 // loop over segments for projection
643 double distance;
644 for (int i = 0; i < size() - 1; i++)
645 {
646 // skip if not the closest segment, note that often two segments are equally close in their shared end point
647 if (d[i] > minD + FRAC_PROJ_PRECISION)
648 {
649 continue;
650 }
651 Point2d center = this.fractionalHelperCenters[i];
652 Point2d p;
653 if (center != null)
654 {
655 // get intersection of line "center - (x, y)" and the segment
656 p = intersectionOfLines(center, point, this.line2d.get(i), this.line2d.get(i + 1));
657 if (p == null || (x < center.x + FRAC_PROJ_PRECISION && center.x + FRAC_PROJ_PRECISION < p.x)
658 || (x > center.x - FRAC_PROJ_PRECISION && center.x - FRAC_PROJ_PRECISION > p.x)
659 || (y < center.y + FRAC_PROJ_PRECISION && center.y + FRAC_PROJ_PRECISION < p.y)
660 || (y > center.y - FRAC_PROJ_PRECISION && center.y - FRAC_PROJ_PRECISION > p.y))
661 {
662 // projected point may not be 'beyond' segment center (i.e. center may not be between (x, y) and (p.x, p.y)
663 continue;
664 }
665 }
666 else
667 {
668 // parallel helper lines, project along direction
669 Point2d offsetPoint =
670 new Point2d(x + this.fractionalHelperDirections[i].x, y + this.fractionalHelperDirections[i].y);
671 p = intersectionOfLines(point, offsetPoint, this.line2d.get(i), this.line2d.get(i + 1));
672 }
673 double segLength = this.line2d.get(i).distance(this.line2d.get(i + 1)) + FRAC_PROJ_PRECISION;
674 if (p == null || this.line2d.get(i).distance(p) > segLength || this.line2d.get(i + 1).distance(p) > segLength)
675 {
676 // intersection must be on the segment
677 // in case of p == null, the length of the fractional helper direction falls away due to precision
678 continue;
679 }
680 // distance from (x, y) to intersection on segment
681 double dx = x - p.x;
682 double dy = y - p.y;
683 distance = Math.hypot(dx, dy);
684 // distance from start of segment to point on segment
685 if (distance < minDistance)
686 {
687 dx = p.x - this.line2d.get(i).x;
688 dy = p.y - this.line2d.get(i).y;
689 double dFrac = Math.hypot(dx, dy);
690 // fraction to point on segment
691 minDistance = distance;
692 minSegmentFraction = dFrac / (this.lengthIndexedLine[i + 1] - this.lengthIndexedLine[i]);
693 minSegment = i;
694 }
695 }
696
697 // return
698 if (minSegment == -1)
699
700 {
701 /*
702 * If fractional projection fails (x, y) is either outside of the applicable area for fractional projection, or is
703 * inside an area where numerical difficulties arise (i.e. far away outside of very slight bend which is considered
704 * parallel).
705 */
706 // CategoryLogger.info(Cat.CORE, "projectFractional failed to project " + point + " on " + this
707 // + "; using fallback approach");
708 return fallback.getFraction(this, x, y);
709 }
710
711 double segLen = this.lengthIndexedLine[minSegment + 1] - this.lengthIndexedLine[minSegment];
712 return (this.lengthIndexedLine[minSegment] + segLen * minSegmentFraction) / this.length.si;
713
714 }
715
716 /**
717 * Fallback method for when fractional projection fails as the point is beyond the line or from numerical limitations.
718 * <p>
719 * Copyright (c) 2013-2024 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved.
720 * <br>
721 * BSD-style license. See <a href="https://opentrafficsim.org/docs/license.html">OpenTrafficSim License</a>.
722 * </p>
723 * @author <a href="https://github.com/averbraeck">Alexander Verbraeck</a>
724 * @author <a href="https://tudelft.nl/staff/p.knoppers-1">Peter Knoppers</a>
725 * @author <a href="https://github.com/wjschakel">Wouter Schakel</a>
726 */
727 public enum FractionalFallback
728 {
729 /** Orthogonal projection. */
730 ORTHOGONAL
731 {
732 @Override
733 double getFraction(final OtsLine2d line, final double x, final double y)
734 {
735 return line.projectOrthogonal(x, y);
736 }
737 },
738
739 /** Distance to nearest end point. */
740 ENDPOINT
741 {
742 @Override
743 double getFraction(final OtsLine2d line, final double x, final double y)
744 {
745 Point2d point = new Point2d(x, y);
746 double dStart = point.distance(line.getFirst());
747 double dEnd = point.distance(line.getLast());
748 if (dStart < dEnd)
749 {
750 return -dStart / line.length.si;
751 }
752 else
753 {
754 return (dEnd + line.length.si) / line.length.si;
755 }
756 }
757 },
758
759 /** NaN value. */
760 NaN
761 {
762 @Override
763 double getFraction(final OtsLine2d line, final double x, final double y)
764 {
765 return Double.NaN;
766 }
767 };
768
769 /**
770 * Returns fraction for when fractional projection fails as the point is beyond the line or from numerical limitations.
771 * @param line OtsLine2d; line
772 * @param x double; x coordinate of point
773 * @param y double; y coordinate of point
774 * @return double; fraction for when fractional projection fails
775 */
776 abstract double getFraction(OtsLine2d line, double x, double y);
777
778 }
779
780 /**
781 * Determines all helpers (points and/or directions) for fractional projection and stores fixed information in properties
782 * while returning the first and last center points (.
783 * @param start Direction; direction in first point
784 * @param end Direction; direction in last point
785 */
786 private synchronized void determineFractionalHelpers(final Direction start, final Direction end)
787 {
788
789 final int n = size() - 1;
790
791 // calculate fixed helpers if not done yet
792 if (this.fractionalHelperCenters == null)
793 {
794 this.fractionalHelperCenters = new Point2d[n];
795 this.fractionalHelperDirections = new Point2D.Double[n];
796 if (size() > 2)
797 {
798 // intersection of parallel lines of first and second segment
799 PolyLine2d prevOfsSeg = unitOffsetSegment(0);
800 PolyLine2d nextOfsSeg = unitOffsetSegment(1);
801 Point2d parStartPoint;
802 parStartPoint = intersectionOfLines(prevOfsSeg.get(0), prevOfsSeg.get(1), nextOfsSeg.get(0), nextOfsSeg.get(1));
803 if (parStartPoint == null || prevOfsSeg.get(1).distance(nextOfsSeg.get(0)) < Math
804 .min(prevOfsSeg.get(1).distance(parStartPoint), nextOfsSeg.get(0).distance(parStartPoint)))
805 {
806 parStartPoint = new Point2d((prevOfsSeg.get(1).x + nextOfsSeg.get(0).x) / 2,
807 (prevOfsSeg.get(1).y + nextOfsSeg.get(0).y) / 2);
808 }
809 // remember the intersection of the first two unit offset segments
810 this.firstOffsetIntersection = parStartPoint;
811 // loop segments
812 for (int i = 1; i < size() - 2; i++)
813 {
814 prevOfsSeg = nextOfsSeg;
815 nextOfsSeg = unitOffsetSegment(i + 1);
816 Point2d parEndPoint;
817 parEndPoint =
818 intersectionOfLines(prevOfsSeg.get(0), prevOfsSeg.get(1), nextOfsSeg.get(0), nextOfsSeg.get(1));
819 if (parEndPoint == null || prevOfsSeg.get(1).distance(nextOfsSeg.get(0)) < Math
820 .min(prevOfsSeg.get(1).distance(parEndPoint), nextOfsSeg.get(0).distance(parEndPoint)))
821 {
822 parEndPoint = new Point2d((prevOfsSeg.get(1).x + nextOfsSeg.get(0).x) / 2,
823 (prevOfsSeg.get(1).y + nextOfsSeg.get(0).y) / 2);
824 }
825 // center = intersections of helper lines
826 this.fractionalHelperCenters[i] =
827 intersectionOfLines(this.line2d.get(i), parStartPoint, this.line2d.get(i + 1), parEndPoint);
828 if (this.fractionalHelperCenters[i] == null)
829 {
830 // parallel helper lines, parallel segments or /\/ cause parallel helper lines, use direction
831 this.fractionalHelperDirections[i] = new Point2D.Double(parStartPoint.x - this.line2d.get(i).x,
832 parStartPoint.y - this.line2d.get(i).y);
833 }
834 parStartPoint = parEndPoint;
835 }
836 // remember the intersection of the last two unit offset segments
837 this.lastOffsetIntersection = parStartPoint;
838 }
839 }
840
841 // use directions at start and end to get unit offset points to the left at a distance of 1
842 double ang = (start == null
843 ? Math.atan2(this.line2d.get(1).y - this.line2d.get(0).y, this.line2d.get(1).x - this.line2d.get(0).x)
844 : start.getInUnit(DirectionUnit.DEFAULT)) + Math.PI / 2; // start.si + Math.PI / 2;
845 Point2d p1 = new Point2d(this.line2d.get(0).x + Math.cos(ang), this.line2d.get(0).y + Math.sin(ang));
846 ang = (end == null
847 ? Math.atan2(this.line2d.get(n).y - this.line2d.get(n - 1).y, this.line2d.get(n).x - this.line2d.get(n - 1).x)
848 : end.getInUnit(DirectionUnit.DEFAULT)) + Math.PI / 2; // end.si + Math.PI / 2;
849 Point2d p2 = new Point2d(this.line2d.get(n).x + Math.cos(ang), this.line2d.get(n).y + Math.sin(ang));
850
851 // calculate first and last center (i.e. intersection of unit offset segments), which depend on inputs 'start' and 'end'
852 if (size() > 2)
853 {
854 this.fractionalHelperCenters[0] =
855 intersectionOfLines(this.line2d.get(0), p1, this.line2d.get(1), this.firstOffsetIntersection);
856 this.fractionalHelperCenters[n - 1] =
857 intersectionOfLines(this.line2d.get(n - 1), this.lastOffsetIntersection, this.line2d.get(n), p2);
858 if (this.fractionalHelperCenters[n - 1] == null)
859 {
860 // parallel helper lines, use direction for projection
861 this.fractionalHelperDirections[n - 1] =
862 new Point2D.Double(p2.x - this.line2d.get(n).x, p2.y - this.line2d.get(n).y);
863 }
864 }
865 else
866 {
867 // only a single segment
868 this.fractionalHelperCenters[0] = intersectionOfLines(this.line2d.get(0), p1, this.line2d.get(1), p2);
869 }
870 if (this.fractionalHelperCenters[0] == null)
871 {
872 // parallel helper lines, use direction for projection
873 this.fractionalHelperDirections[0] = new Point2D.Double(p1.x - this.line2d.get(0).x, p1.y - this.line2d.get(0).y);
874 }
875
876 }
877
878 /**
879 * This method is used, rather than {@code Point2d.intersectionOfLines()} because this method will return {@code null} if
880 * the determinant < 0.0000001, rather than determinant &eq; 0.0. The benefit of this is that intersections are not so
881 * far away, that any calculations with them cause underflow or overflow issues.
882 * @param line1P1 Point2d; point 1 of line 1.
883 * @param line1P2 Point2d; point 2 of line 1.
884 * @param line2P1 Point2d; point 1 of line 2.
885 * @param line2P2 Point2d; point 2 of line 2.
886 * @return Point2d; intersection of lines, or {@code null} for (nearly) parallel lines.
887 */
888 private Point2d intersectionOfLines(final Point2d line1P1, final Point2d line1P2, final Point2d line2P1,
889 final Point2d line2P2)
890 {
891 double l1p1x = line1P1.x;
892 double l1p1y = line1P1.y;
893 double l1p2x = line1P2.x - l1p1x;
894 double l1p2y = line1P2.y - l1p1y;
895 double l2p1x = line2P1.x - l1p1x;
896 double l2p1y = line2P1.y - l1p1y;
897 double l2p2x = line2P2.x - l1p1x;
898 double l2p2y = line2P2.y - l1p1y;
899 double determinant = (0 - l1p2x) * (l2p1y - l2p2y) - (0 - l1p2y) * (l2p1x - l2p2x);
900 if (Math.abs(determinant) < 0.0000001)
901 {
902 return null;
903 }
904 return new Point2d(l1p1x + (l1p2x * (l2p1x * l2p2y - l2p1y * l2p2x)) / determinant,
905 l1p1y + (l1p2y * (l2p1x * l2p2y - l2p1y * l2p2x)) / determinant);
906 }
907
908 /**
909 * Helper method for fractional projection which returns an offset line to the left of a segment at a distance of 1.
910 * @param segment int; segment number
911 * @return parallel line to the left of a segment at a distance of 1
912 */
913 private synchronized PolyLine2d unitOffsetSegment(final int segment)
914 {
915 return new PolyLine2d(this.line2d.get(segment), this.line2d.get(segment + 1)).offsetLine(1.0);
916 }
917
918 /**
919 * Returns the projected directional radius of the line at a given fraction. Negative values reflect right-hand curvature in
920 * the design-line direction. The radius is taken as the minimum of the radii at the vertices before and after the given
921 * fraction. The radius at a vertex is calculated as the radius of a circle that is equidistant from both edges connected to
922 * the vertex. The circle center is on a line perpendicular to the shortest edge, crossing through the middle of the
923 * shortest edge. This method ignores Z components.
924 * @param fraction double; fraction along the line, between 0.0 and 1.0 (both inclusive)
925 * @return Length; radius; the local radius; or si field set to Double.NaN if line is totally straight
926 * @throws OtsGeometryException fraction out of bounds
927 */
928 // TODO: move to djutils?
929 public synchronized Length getProjectedRadius(final double fraction) throws OtsGeometryException
930 {
931 Throw.when(fraction < 0.0 || fraction > 1.0, OtsGeometryException.class, "Fraction %f is out of bounds [0.0 ... 1.0]",
932 fraction);
933 if (this.vertexRadii == null)
934 {
935 this.vertexRadii = new Length[size() - 1];
936 }
937 int index = find(fraction * getLength().si);
938 if (index > 0 && this.vertexRadii[index] == null)
939 {
940 this.vertexRadii[index] = getProjectedVertexRadius(index);
941 }
942 if (index < size() - 2 && this.vertexRadii[index + 1] == null)
943 {
944 this.vertexRadii[index + 1] = getProjectedVertexRadius(index + 1);
945 }
946 if (index == 0)
947 {
948 if (this.vertexRadii.length < 2)
949 {
950 return Length.instantiateSI(Double.NaN);
951 }
952 return this.vertexRadii[1];
953 }
954 if (index == size() - 2)
955 {
956 return this.vertexRadii[size() - 2];
957 }
958 return Math.abs(this.vertexRadii[index].si) < Math.abs(this.vertexRadii[index + 1].si) ? this.vertexRadii[index]
959 : this.vertexRadii[index + 1];
960 }
961
962 /**
963 * Calculates the directional radius at a vertex. Negative values reflect right-hand curvature in the design-line direction.
964 * The radius at a vertex is calculated as the radius of a circle that is equidistant from both edges connected to the
965 * vertex. The circle center is on a line perpendicular to the shortest edge, crossing through the middle of the shortest
966 * edge. This function ignores Z components.
967 * @param index int; index of the vertex in range [1 ... size() - 2]
968 * @return Length; radius at the vertex
969 * @throws OtsGeometryException if the index is out of bounds
970 */
971 // TODO: move to djutils? Note, uses fractionalHelperCenters
972 public synchronized Length getProjectedVertexRadius(final int index) throws OtsGeometryException
973 {
974 Throw.when(index < 1 || index > size() - 2, OtsGeometryException.class, "Index %d is out of bounds [1 ... size() - 2].",
975 index);
976 determineFractionalHelpers(null, null);
977 double length1 = this.lengthIndexedLine[index] - this.lengthIndexedLine[index - 1];
978 double length2 = this.lengthIndexedLine[index + 1] - this.lengthIndexedLine[index];
979 int shortIndex = length1 < length2 ? index : index + 1;
980 // center of shortest edge
981 Point2d p1 = new Point2d(.5 * (this.line2d.get(shortIndex - 1).x + this.line2d.get(shortIndex).x),
982 .5 * (this.line2d.get(shortIndex - 1).y + this.line2d.get(shortIndex).y));
983 // perpendicular to shortest edge, line crossing p1
984 Point2d p2 = new Point2d(p1.x + (this.line2d.get(shortIndex).y - this.line2d.get(shortIndex - 1).y),
985 p1.y - (this.line2d.get(shortIndex).x - this.line2d.get(shortIndex - 1).x));
986 // vertex
987 Point2d p3 = this.line2d.get(index);
988 // point on line that splits angle between edges at vertex 50-50
989 Point2d p4 = this.fractionalHelperCenters[index];
990 if (p4 == null)
991 {
992 // parallel helper lines
993 p4 = new Point2d(p3.x + this.fractionalHelperDirections[index].x, p3.y + this.fractionalHelperDirections[index].y);
994 }
995 Point2d intersection = intersectionOfLines(p1, p2, p3, p4);
996 if (null == intersection)
997 {
998 return Length.instantiateSI(Double.NaN);
999 }
1000 // determine left or right
1001 double refLength = length1 < length2 ? length1 : length2;
1002 double radius = intersection.distance(p1);
1003 double i2p2 = intersection.distance(p2);
1004 if (radius < i2p2 && i2p2 > refLength)
1005 {
1006 // left as p1 is closer than p2 (which was placed to the right) and not on the perpendicular line
1007 return Length.instantiateSI(radius);
1008 }
1009 // right as not left
1010 return Length.instantiateSI(-radius);
1011 }
1012
1013 /**
1014 * Returns the length fraction at the vertex.
1015 * @param index int; index of vertex [0 ... size() - 1]
1016 * @return double; length fraction at the vertex
1017 * @throws OtsGeometryException if the index is out of bounds
1018 */
1019 public double getVertexFraction(final int index) throws OtsGeometryException
1020 {
1021 Throw.when(index < 0 || index > size() - 1, OtsGeometryException.class, "Index %d is out of bounds [0 %d].", index,
1022 size() - 1);
1023 return this.lengthIndexedLine[index] / this.length.si;
1024 }
1025
1026 /**
1027 * Retrieve the centroid of this OtsLine2d.
1028 * @return OtsPoint3d; the centroid of this OtsLine2d
1029 */
1030 public final Point2d getCentroid()
1031 {
1032 if (this.centroid == null)
1033 {
1034 this.centroid = this.line2d.getBounds().midPoint();
1035 }
1036 return this.centroid;
1037 }
1038
1039 /**
1040 * Get the bounding rectangle of this OtsLine2d.
1041 * @return Rectangle2D; the bounding rectangle of this OtsLine2d
1042 */
1043 public final Bounds2d getEnvelope()
1044 {
1045 return this.line2d.getBounds();
1046 }
1047
1048 /** {@inheritDoc} */
1049 @Override
1050 @SuppressWarnings("checkstyle:designforextension")
1051 public Point2d getLocation()
1052 {
1053 return getCentroid();
1054 }
1055
1056 /** {@inheritDoc} */
1057 @Override
1058 @SuppressWarnings("checkstyle:designforextension")
1059 public Bounds2d getBounds()
1060 {
1061 if (this.bounds == null)
1062 {
1063 Bounds2d envelope = getEnvelope();
1064 this.bounds = new Bounds2d(envelope.getDeltaX(), envelope.getDeltaY());
1065 }
1066 return this.bounds;
1067 }
1068
1069 /** {@inheritDoc} */
1070 @Override
1071 @SuppressWarnings("checkstyle:designforextension")
1072 public String toString()
1073 {
1074 return this.line2d.toString();
1075 }
1076
1077 /** {@inheritDoc} */
1078 @Override
1079 @SuppressWarnings("checkstyle:designforextension")
1080 public int hashCode()
1081 {
1082 return this.line2d.hashCode();
1083 }
1084
1085 /** {@inheritDoc} */
1086 @Override
1087 @SuppressWarnings({"checkstyle:designforextension", "checkstyle:needbraces"})
1088 public boolean equals(final Object obj)
1089 {
1090 if (!(obj instanceof OtsLine2d))
1091 {
1092 return false;
1093 }
1094 return this.line2d.equals(((OtsLine2d) obj).line2d);
1095 }
1096
1097 /**
1098 * Convert the 2D projection of this OtsLine2d to something that MS-Excel can plot.
1099 * @return excel XY plottable output
1100 */
1101 public final String toExcel()
1102 {
1103 return this.line2d.toExcel();
1104 }
1105
1106 /**
1107 * Convert the 2D projection of this OtsLine2d to Peter's plot format.
1108 * @return Peter's format plot output
1109 */
1110 public final String toPlot()
1111 {
1112 return this.line2d.toPlot();
1113 }
1114
1115 }