View Javadoc
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      *            &lt;strong&gt;must&lt;/strong&gt; 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      *            &lt;strong&gt;must&lt;/strong&gt; 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 &lt; 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&lt;Point2d&gt;; 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 &lt; 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&lt;OtsPoint3d&gt;.
312      * @param pointList List&lt;OtsPoint3d&gt;; 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 &lt; 0 or i &gt; 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 &lt; 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 }