View Javadoc
1   package org.opentrafficsim.base.geometry;
2   
3   import java.awt.geom.Line2D;
4   import java.awt.geom.Point2D;
5   import java.util.ArrayList;
6   import java.util.Iterator;
7   import java.util.List;
8   
9   import org.djunits.value.vdouble.scalar.Direction;
10  import org.djunits.value.vdouble.scalar.Length;
11  import org.djutils.draw.bounds.Bounds;
12  import org.djutils.draw.line.PolyLine2d;
13  import org.djutils.draw.point.DirectedPoint2d;
14  import org.djutils.draw.point.Point2d;
15  import org.djutils.exceptions.Throw;
16  
17  import nl.tudelft.simulation.dsol.animation.Locatable;
18  
19  /**
20   * This class supports fractional projection, radius, and has locatable methods.
21   * <p>
22   * Copyright (c) 2013-2024 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
23   * BSD-style license. See <a href="https://opentrafficsim.org/docs/license.html">OpenTrafficSim License</a>.
24   * </p>
25   * @author <a href="https://github.com/averbraeck">Alexander Verbraeck</a>
26   * @author <a href="https://github.com/peter-knoppers">Peter Knoppers</a>
27   * @author <a href="https://www.citg.tudelft.nl">Guus Tamminga</a>
28   * @author <a href="https://github.com/wjschakel">Wouter Schakel</a>
29   */
30  public class OtsLine2d extends PolyLine2d implements Locatable
31  {
32      /** The cached typed length; will be calculated at time of construction. */
33      private final Length length;
34  
35      // Fractional projection
36  
37      /** The cached helper points for fractional projection; will be calculated when needed for the first time. */
38      private Point2d[] fractionalHelperCenters = null;
39  
40      /** The cached helper directions for fractional projection; will be calculated when needed for the first time. */
41      private Point2D.Double[] fractionalHelperDirections = null;
42  
43      /** Last start direction for which fractional helpers were calculated. */
44      private Direction lastStartDirection;
45  
46      /** Last start direction for which fractional helpers were calculated. */
47      private Direction lastEndDirection;
48  
49      /** Intersection of unit offset lines of first two segments. */
50      private Point2d firstOffsetIntersection;
51  
52      /** Intersection of unit offset lines of last two segments. */
53      private Point2d lastOffsetIntersection;
54  
55      /** Precision for fractional projection algorithm. */
56      private static final double FRAC_PROJ_PRECISION = 2e-5 /* PK too fine 1e-6 */;
57  
58      // Radii for curvature
59  
60      /** Radius at each vertex. */
61      private Length[] vertexRadii;
62  
63      /**
64       * Construct a new OtsLine2d.
65       * @param points the array of points to construct this OtsLine2d from.
66       */
67      public OtsLine2d(final Point2d... points)
68      {
69          super(points);
70          this.length = Length.ofSI(lengthAtIndex(size() - 1));
71      }
72  
73      /**
74       * Creates a new OtsLine2d based on 2d information.
75       * @param line2d 2d information.
76       */
77      public OtsLine2d(final PolyLine2d line2d)
78      {
79          super(line2d.iterator());
80          this.length = Length.ofSI(lengthAtIndex(size() - 1));
81      }
82  
83      /**
84       * Creates a new OtsLine2d based on point iterator.
85       * @param line2d point iterator.
86       */
87      public OtsLine2d(final Iterator<Point2d> line2d)
88      {
89          super(line2d);
90          this.length = Length.ofSI(lengthAtIndex(size() - 1));
91      }
92  
93      /**
94       * Construct a new OtsLine2d from a List&lt;OtsPoint3d&gt;.
95       * @param pointList the list of points to construct this OtsLine2d from.
96       */
97      public OtsLine2d(final List<Point2d> pointList)
98      {
99          super(pointList);
100         this.length = Length.ofSI(lengthAtIndex(size() - 1));
101     }
102 
103     /**
104      * Construct parallel line.
105      * @param offset offset distance from the reference line; positive is LEFT, negative is RIGHT
106      * @return the line that has the specified offset from the reference line
107      */
108     @Override
109     public final OtsLine2d offsetLine(final double offset)
110     {
111         return new OtsLine2d(super.offsetLine(offset));
112     }
113 
114     /**
115      * Create a line at linearly varying offset from this line. The offset may change linearly from its initial value at the
116      * start of the reference line to its final offset value at the end of the reference line.
117      * @param offsetAtStart offset at the start of the reference line (positive value is Left, negative value is Right)
118      * @param offsetAtEnd offset at the end of the reference line (positive value is Left, negative value is Right)
119      * @return the OtsLine2d of the line at linearly changing offset of the reference line
120      */
121     @Override
122     public final OtsLine2d offsetLine(final double offsetAtStart, final double offsetAtEnd)
123     {
124         return new OtsLine2d(super.offsetLine(offsetAtStart, offsetAtEnd).getPointList());
125     }
126 
127     /**
128      * Create a line at linearly varying offset from this line. The offset may change linearly from its initial value at the
129      * start of the reference line via a number of intermediate offsets at intermediate positions to its final offset value at
130      * the end of the reference line.
131      * @param relativeFractions positional fractions for which the offsets have to be generated
132      * @param offsets offsets at the relative positions (positive value is Left, negative value is Right)
133      * @return the Geometry of the line at linearly changing offset of the reference line
134      */
135     public final OtsLine2d offsetLine(final double[] relativeFractions, final double[] offsets)
136     {
137         return new OtsLine2d(OtsGeometryUtil.offsetLine(this, relativeFractions, offsets));
138     }
139 
140     /**
141      * Concatenate several OtsLine2d instances.
142      * @param lines OtsLine2d... one or more OtsLine2d. The last point of the first &lt;strong&gt;must&lt;/strong&gt; match the
143      *            first of the second, etc.
144      * @return OtsLine2d
145      */
146     public static OtsLine2d concatenate(final OtsLine2d... lines)
147     {
148         return concatenate(0.0, lines);
149     }
150 
151     /**
152      * Concatenate two OtsLine2d instances. This method is separate for efficiency reasons.
153      * @param toleranceSI the tolerance between the end point of a line and the first point of the next line
154      * @param line1 first line
155      * @param line2 second line
156      * @return OtsLine2d
157      */
158     public static OtsLine2d concatenate(final double toleranceSI, final OtsLine2d line1, final OtsLine2d line2)
159     {
160         return new OtsLine2d(PolyLine2d.concatenate(toleranceSI, line1, line2));
161     }
162 
163     /**
164      * Concatenate several OtsLine2d instances.
165      * @param toleranceSI the tolerance between the end point of a line and the first point of the next line
166      * @param lines OtsLine2d... one or more OtsLine2d. The last point of the first &lt;strong&gt;must&lt;/strong&gt; match the
167      *            first of the second, etc.
168      * @return OtsLine2d
169      */
170     public static OtsLine2d concatenate(final double toleranceSI, final OtsLine2d... lines)
171     {
172         List<PolyLine2d> lines2d = new ArrayList<>();
173         for (OtsLine2d line : lines)
174         {
175             lines2d.add(line);
176         }
177         return new OtsLine2d(PolyLine2d.concatenate(toleranceSI, lines2d.toArray(new PolyLine2d[lines.length])));
178     }
179 
180     /**
181      * Construct a new OtsLine2d with all points of this OtsLine2d in reverse order.
182      * @return the new OtsLine2d
183      */
184     @Override
185     public final OtsLine2d reverse()
186     {
187         return new OtsLine2d(super.reverse());
188     }
189 
190     /**
191      * Construct a new OtsLine2d covering the indicated fraction of this OtsLine2d.
192      * @param start starting point, valid range [0..<cite>end</cite>)
193      * @param end ending point, valid range (<cite>start</cite>..1]
194      * @return the new OtsLine2d
195      */
196     @Override
197     public OtsLine2d extractFractional(final double start, final double end)
198     {
199         return extract(start * this.length.si, end * this.length.si);
200     }
201 
202     /**
203      * Create a new OtsLine2d that covers a sub-section of this OtsLine2d.
204      * @param start the length along this OtsLine2d where the sub-section starts, valid range [0..<cite>end</cite>)
205      * @param end length along this OtsLine2d where the sub-section ends, valid range (<cite>start</cite>..<cite>length</cite>
206      *            (length is the length of this OtsLine2d)
207      * @return the selected sub-section
208      */
209     public final OtsLine2d extract(final Length start, final Length end)
210     {
211         return extract(start.si, end.si);
212     }
213 
214     /**
215      * Create a new OtsLine2d that covers a sub-section of this OtsLine2d.
216      * @param start length along this OtsLine2d where the sub-section starts, valid range [0..<cite>end</cite>)
217      * @param end length along this OtsLine2d where the sub-section ends, valid range (<cite>start</cite>..<cite>length</cite>
218      *            (length is the length of this OtsLine2d)
219      * @return the selected sub-section
220      */
221     @Override
222     public final OtsLine2d extract(final double start, final double end)
223     {
224         return new OtsLine2d(super.extract(start, end));
225     }
226 
227     /**
228      * Return the length of this OtsLine2d in meters. (Assuming that the coordinates of the points constituting this line are
229      * expressed in meters.)
230      * @return the length of the line
231      */
232     public final Length getTypedLength()
233     {
234         return this.length;
235     }
236 
237     /**
238      * Get the location at a position on the line, with its direction. Position can be below 0 or more than the line length. In
239      * that case, the position will be extrapolated in the direction of the line at its start or end.
240      * @param position the position on the line for which to calculate the point on, before, of after the line
241      * @return a directed point
242      */
243     public final DirectedPoint2d getLocationExtended(final Length position)
244     {
245         return getLocationExtended(position.si);
246     }
247 
248     /**
249      * Get the location at a position on the line, with its direction. Position can be below 0 or more than the line length. In
250      * that case, the position will be extrapolated in the direction of the line at its start or end.
251      * @param positionSI the position on the line for which to calculate the point on, before, of after the line, in SI units
252      * @return a directed point
253      */
254     public final DirectedPoint2d getLocationExtendedSI(final double positionSI)
255     {
256         return getLocationExtended(positionSI);
257     }
258 
259     /**
260      * Get the location at a fraction of the line, with its direction. Fraction should be between 0.0 and 1.0.
261      * @param fraction the fraction for which to calculate the point on the line
262      * @return a directed point
263      */
264     public final DirectedPoint2d getLocationPointFraction(final double fraction)
265     {
266         return getLocationFraction(fraction);
267     }
268 
269     /**
270      * Get the location at a fraction of the line, with its direction. Fraction should be between 0.0 and 1.0.
271      * @param fraction the fraction for which to calculate the point on the line
272      * @param tolerance the delta from 0.0 and 1.0 that will be forgiven
273      * @return a directed point
274      */
275     public final DirectedPoint2d getLocationPointFraction(final double fraction, final double tolerance)
276     {
277         return getLocationFraction(fraction, tolerance);
278     }
279 
280     /**
281      * Get the location at a fraction of the line (or outside the line), with its direction.
282      * @param fraction the fraction for which to calculate the point on the line
283      * @return a directed point
284      */
285     public final DirectedPoint2d getLocationPointFractionExtended(final double fraction)
286     {
287         return getLocationFractionExtended(fraction);
288     }
289 
290     /**
291      * Get the location at a position on the line, with its direction. Position should be between 0.0 and line length.
292      * @param position the position on the line for which to calculate the point on the line
293      * @return a directed point
294      */
295     public final DirectedPoint2d getLocation(final Length position)
296     {
297         return getLocation(position.si);
298     }
299 
300     /**
301      * Get the location at a position on the line, with its direction. Position should be between 0.0 and line length.
302      * @param positionSI the position on the line for which to calculate the point on the line
303      * @return a directed point
304      */
305     public final DirectedPoint2d getLocationSI(final double positionSI)
306     {
307         return getLocation(positionSI);
308     }
309 
310     /**
311      * Truncate a line at the given length (less than the length of the line, and larger than zero) and return a new line.
312      * @param lengthSI the location where to truncate the line
313      * @return a new OtsLine2d truncated at the exact position where line.getLength() == lengthSI
314      */
315     @Override
316     public final OtsLine2d truncate(final double lengthSI)
317     {
318         return new OtsLine2d(super.truncate(lengthSI));
319     }
320 
321     /**
322      * Returns the fractional position along this line of the orthogonal projection of point (x, y) on this line. If the point
323      * is not orthogonal to the closest line segment, the nearest point is selected.
324      * @param x x-coordinate of point to project
325      * @param y y-coordinate of point to project
326      * @return fractional position along this line of the orthogonal projection on this line of a point
327      */
328     public final double projectOrthogonalSnap(final double x, final double y)
329     {
330         Point2d closest = closestPointOnPolyLine(new Point2d(x, y));
331         return projectOrthogonalFractionalExtended(closest);
332     }
333 
334     /**
335      * Returns the fractional projection of a point to a line. The projection works by taking slices in space per line segment
336      * as shown below. A point is always projected to the nearest segment, but not necessarily to the closest point on that
337      * segment. The slices in space are analogous to a Voronoi diagram, but for the line segments instead of points. If
338      * fractional projection fails, a fallback projection is returned.<br>
339      * <br>
340      * 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
341      * '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')
342      * in half, while the line 'E-G' cuts the second bend of the 3rd segment (at point 'I') in half.
343      *
344      * <pre>
345      *            ____________________________     G                   .
346      * .         |                            |    .                 .
347      *   .       |  . . . .  helper lines     |    .               .
348      *     .     |  _.._.._  projection line  |   I.             .
349      *       .   |____________________________|  _.'._         .       L
350      *        F.                              _.'  .  '-.    .
351      *          ..                       B _.'     .     '-.
352      *           . .                    _.\        .     .  D
353      *            .  .               _.'   :       .   .
354      *     J       .   .          _.'      \       . .
355      *             ..    .     _.'          :      .                M
356      *            .  .     ..-'             \      .
357      *           .    .    /H.               A     .
358      *          .      .  /    .                   .
359      *        C _________/       .                 .
360      *        .          .         .               .
361      *   K   .            .          .             .
362      *      .              .           .           .
363      *     .                .            .         .           N
364      *    .                  .             .       .
365      *   .                    .              .     .
366      *  .                      .               .   .
367      * .                        .                . .
368      *                           .                 .E
369      *                            .                  .
370      *                             .                   .
371      *                              .                    .
372      * </pre>
373      *
374      * Fractional projection may fail in three cases.
375      * <ol>
376      * <li>Numerical difficulties at slight bend, orthogonal projection returns the correct point.</li>
377      * <li>Fractional projection is possible only to segments that aren't the nearest segment(s).</li>
378      * <li>Fractional projection is possible for no segment.</li>
379      * </ol>
380      * In the latter two cases the projection is undefined and the provided fallback is used to provide a point.
381      * @param start direction in first point
382      * @param end direction in last point
383      * @param x x-coordinate of point to project
384      * @param y y-coordinate of point to project
385      * @param fallback fallback method for when fractional projection fails
386      * @return fractional position along this line of the fractional projection on that line of a point
387      */
388     public final synchronized double projectFractional(final Direction start, final Direction end, final double x,
389             final double y, final FractionalFallback fallback)
390     {
391 
392         // prepare
393         double minDistance = Double.POSITIVE_INFINITY;
394         double minSegmentFraction = 0;
395         int minSegment = -1;
396         Point2d point = new Point2d(x, y);
397 
398         // determine helpers (centers and directions)
399         determineFractionalHelpers(start, end);
400 
401         // get distance of point to each segment
402         double[] d = new double[size() - 1];
403         double minD = Double.POSITIVE_INFINITY;
404         for (int i = 0; i < size() - 1; i++)
405         {
406             d[i] = Line2D.ptSegDist(get(i).x, get(i).y, get(i + 1).x, get(i + 1).y, x, y);
407             minD = d[i] < minD ? d[i] : minD;
408         }
409 
410         // loop over segments for projection
411         double distance;
412         for (int i = 0; i < size() - 1; i++)
413         {
414             // skip if not the closest segment, note that often two segments are equally close in their shared end point
415             if (d[i] > minD + FRAC_PROJ_PRECISION)
416             {
417                 continue;
418             }
419             Point2d center = this.fractionalHelperCenters[i];
420             Point2d p;
421             if (center != null)
422             {
423                 // get intersection of line "center - (x, y)" and the segment
424                 p = intersectionOfLines(center, point, get(i), get(i + 1));
425                 if (p == null || (x < center.x + FRAC_PROJ_PRECISION && center.x + FRAC_PROJ_PRECISION < p.x)
426                         || (x > center.x - FRAC_PROJ_PRECISION && center.x - FRAC_PROJ_PRECISION > p.x)
427                         || (y < center.y + FRAC_PROJ_PRECISION && center.y + FRAC_PROJ_PRECISION < p.y)
428                         || (y > center.y - FRAC_PROJ_PRECISION && center.y - FRAC_PROJ_PRECISION > p.y))
429                 {
430                     // projected point may not be 'beyond' segment center (i.e. center may not be between (x, y) and (p.x, p.y)
431                     continue;
432                 }
433             }
434             else
435             {
436                 // parallel helper lines, project along direction
437                 Point2d offsetPoint =
438                         new Point2d(x + this.fractionalHelperDirections[i].x, y + this.fractionalHelperDirections[i].y);
439                 p = intersectionOfLines(point, offsetPoint, get(i), get(i + 1));
440             }
441             double segLength = get(i).distance(get(i + 1)) + FRAC_PROJ_PRECISION;
442             if (p == null || get(i).distance(p) > segLength || get(i + 1).distance(p) > segLength)
443             {
444                 // intersection must be on the segment
445                 // in case of p == null, the length of the fractional helper direction falls away due to precision
446                 continue;
447             }
448             // distance from (x, y) to intersection on segment
449             double dx = x - p.x;
450             double dy = y - p.y;
451             distance = Math.hypot(dx, dy);
452             // distance from start of segment to point on segment
453             if (distance < minDistance)
454             {
455                 dx = p.x - get(i).x;
456                 dy = p.y - get(i).y;
457                 double dFrac = Math.hypot(dx, dy);
458                 // fraction to point on segment
459                 minDistance = distance;
460                 minSegmentFraction = dFrac / (lengthAtIndex(i + 1) - lengthAtIndex(i));
461                 minSegment = i;
462             }
463         }
464 
465         // return
466         if (minSegment == -1)
467 
468         {
469             /*
470              * If fractional projection fails (x, y) is either outside of the applicable area for fractional projection, or is
471              * inside an area where numerical difficulties arise (i.e. far away outside of very slight bend which is considered
472              * parallel).
473              */
474             return fallback.getFraction(this, x, y);
475         }
476 
477         double segLen = lengthAtIndex(minSegment + 1) - lengthAtIndex(minSegment);
478         return (lengthAtIndex(minSegment) + segLen * minSegmentFraction) / this.length.si;
479 
480     }
481 
482     /**
483      * Fallback method for when fractional projection fails as the point is beyond the line or from numerical limitations.
484      * <p>
485      * Copyright (c) 2013-2024 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved.
486      * <br>
487      * BSD-style license. See <a href="https://opentrafficsim.org/docs/license.html">OpenTrafficSim License</a>.
488      * </p>
489      * @author <a href="https://github.com/averbraeck">Alexander Verbraeck</a>
490      * @author <a href="https://github.com/peter-knoppers">Peter Knoppers</a>
491      * @author <a href="https://github.com/wjschakel">Wouter Schakel</a>
492      */
493     public enum FractionalFallback
494     {
495         /** Orthogonal projection. */
496         ORTHOGONAL
497         {
498             @Override
499             double getFraction(final OtsLine2d line, final double x, final double y)
500             {
501                 return line.projectOrthogonalSnap(x, y);
502             }
503         },
504 
505         /** Orthogonal projection with extension. */
506         ORTHOGONAL_EXTENDED
507         {
508             @Override
509             double getFraction(final OtsLine2d line, final double x, final double y)
510             {
511                 return line.projectOrthogonalFractionalExtended(new Point2d(x, y));
512             }
513         },
514 
515         /** Distance to nearest end point. */
516         ENDPOINT
517         {
518             @Override
519             double getFraction(final OtsLine2d line, final double x, final double y)
520             {
521                 Point2d point = new Point2d(x, y);
522                 double dStart = point.distance(line.getFirst());
523                 double dEnd = point.distance(line.getLast());
524                 if (dStart < dEnd)
525                 {
526                     return -dStart / line.length.si;
527                 }
528                 else
529                 {
530                     return (dEnd + line.length.si) / line.length.si;
531                 }
532             }
533         },
534 
535         /** NaN value. */
536         NaN
537         {
538             @Override
539             double getFraction(final OtsLine2d line, final double x, final double y)
540             {
541                 return Double.NaN;
542             }
543         };
544 
545         /**
546          * Returns fraction for when fractional projection fails as the point is beyond the line or from numerical limitations.
547          * @param line line
548          * @param x x coordinate of point
549          * @param y y coordinate of point
550          * @return fraction for when fractional projection fails
551          */
552         abstract double getFraction(OtsLine2d line, double x, double y);
553 
554     }
555 
556     /**
557      * Determines all helpers (points and/or directions) for fractional projection and stores fixed information in properties
558      * while returning the first and last center points (.
559      * @param start direction in first point
560      * @param end direction in last point
561      */
562     private synchronized void determineFractionalHelpers(final Direction start, final Direction end)
563     {
564 
565         final int n = size() - 1;
566 
567         // calculate fixed helpers if not done yet
568         if (this.fractionalHelperCenters == null)
569         {
570             this.fractionalHelperCenters = new Point2d[n];
571             this.fractionalHelperDirections = new Point2D.Double[n];
572             if (size() > 2)
573             {
574                 // intersection of parallel lines of first and second segment
575                 PolyLine2d prevOfsSeg = unitOffsetSegment(0);
576                 PolyLine2d nextOfsSeg = unitOffsetSegment(1);
577                 Point2d parStartPoint;
578                 parStartPoint = intersectionOfLines(prevOfsSeg.get(0), prevOfsSeg.get(1), nextOfsSeg.get(0), nextOfsSeg.get(1));
579                 if (parStartPoint == null || prevOfsSeg.get(1).distance(nextOfsSeg.get(0)) < Math
580                         .min(prevOfsSeg.get(1).distance(parStartPoint), nextOfsSeg.get(0).distance(parStartPoint)))
581                 {
582                     parStartPoint = new Point2d((prevOfsSeg.get(1).x + nextOfsSeg.get(0).x) / 2,
583                             (prevOfsSeg.get(1).y + nextOfsSeg.get(0).y) / 2);
584                 }
585                 // remember the intersection of the first two unit offset segments
586                 this.firstOffsetIntersection = parStartPoint;
587                 // loop segments
588                 for (int i = 1; i < size() - 2; i++)
589                 {
590                     prevOfsSeg = nextOfsSeg;
591                     nextOfsSeg = unitOffsetSegment(i + 1);
592                     Point2d parEndPoint;
593                     parEndPoint =
594                             intersectionOfLines(prevOfsSeg.get(0), prevOfsSeg.get(1), nextOfsSeg.get(0), nextOfsSeg.get(1));
595                     if (parEndPoint == null || prevOfsSeg.get(1).distance(nextOfsSeg.get(0)) < Math
596                             .min(prevOfsSeg.get(1).distance(parEndPoint), nextOfsSeg.get(0).distance(parEndPoint)))
597                     {
598                         parEndPoint = new Point2d((prevOfsSeg.get(1).x + nextOfsSeg.get(0).x) / 2,
599                                 (prevOfsSeg.get(1).y + nextOfsSeg.get(0).y) / 2);
600                     }
601                     // center = intersections of helper lines
602                     this.fractionalHelperCenters[i] = intersectionOfLines(get(i), parStartPoint, get(i + 1), parEndPoint);
603                     if (this.fractionalHelperCenters[i] == null)
604                     {
605                         // parallel helper lines, parallel segments or /\/ cause parallel helper lines, use direction
606                         this.fractionalHelperDirections[i] =
607                                 new Point2D.Double(parStartPoint.x - get(i).x, parStartPoint.y - get(i).y);
608                     }
609                     parStartPoint = parEndPoint;
610                 }
611                 // remember the intersection of the last two unit offset segments
612                 this.lastOffsetIntersection = parStartPoint;
613             }
614         }
615 
616         // skip when same directions previously used to calculate edge helpers
617         if (this.lastStartDirection != null && this.lastStartDirection.eq(start) && this.lastEndDirection != null
618                 && this.lastEndDirection.eq(end))
619         {
620             return;
621         }
622         this.lastStartDirection = start;
623         this.lastEndDirection = end;
624 
625         // use directions at start and end to get unit offset points to the left at a distance of 1
626         double ang = (start == null ? Math.atan2(get(1).y - get(0).y, get(1).x - get(0).x) : start.si) + Math.PI / 2;
627         Point2d p1 = new Point2d(get(0).x + Math.cos(ang), get(0).y + Math.sin(ang));
628         ang = (end == null ? Math.atan2(get(n).y - get(n - 1).y, get(n).x - get(n - 1).x) : end.si) + Math.PI / 2;
629         Point2d p2 = new Point2d(get(n).x + Math.cos(ang), get(n).y + Math.sin(ang));
630 
631         // calculate first and last center (i.e. intersection of unit offset segments), which depend on inputs 'start' and 'end'
632         if (size() > 2)
633         {
634             this.fractionalHelperCenters[0] = intersectionOfLines(get(0), p1, get(1), this.firstOffsetIntersection);
635             this.fractionalHelperCenters[n - 1] = intersectionOfLines(get(n - 1), this.lastOffsetIntersection, get(n), p2);
636             if (this.fractionalHelperCenters[n - 1] == null)
637             {
638                 // parallel helper lines, use direction for projection
639                 this.fractionalHelperDirections[n - 1] = new Point2D.Double(p2.x - get(n).x, p2.y - get(n).y);
640             }
641         }
642         else
643         {
644             // only a single segment
645             this.fractionalHelperCenters[0] = intersectionOfLines(get(0), p1, get(1), p2);
646         }
647         if (this.fractionalHelperCenters[0] == null)
648         {
649             // parallel helper lines, use direction for projection
650             this.fractionalHelperDirections[0] = new Point2D.Double(p1.x - get(0).x, p1.y - get(0).y);
651         }
652 
653     }
654 
655     /**
656      * This method is used, rather than {@code Point2d.intersectionOfLines()} because this method will return {@code null} if
657      * the determinant &lt; 0.0000001, rather than determinant &eq; 0.0. The benefit of this is that intersections are not so
658      * far away, that any calculations with them cause underflow or overflow issues.
659      * @param line1P1 point 1 of line 1.
660      * @param line1P2 point 2 of line 1.
661      * @param line2P1 point 1 of line 2.
662      * @param line2P2 point 2 of line 2.
663      * @return intersection of lines, or {@code null} for (nearly) parallel lines.
664      */
665     private Point2d intersectionOfLines(final Point2d line1P1, final Point2d line1P2, final Point2d line2P1,
666             final Point2d line2P2)
667     {
668         double l1p1x = line1P1.x;
669         double l1p1y = line1P1.y;
670         double l1p2x = line1P2.x - l1p1x;
671         double l1p2y = line1P2.y - l1p1y;
672         double l2p1x = line2P1.x - l1p1x;
673         double l2p1y = line2P1.y - l1p1y;
674         double l2p2x = line2P2.x - l1p1x;
675         double l2p2y = line2P2.y - l1p1y;
676         double determinant = (0 - l1p2x) * (l2p1y - l2p2y) - (0 - l1p2y) * (l2p1x - l2p2x);
677         if (Math.abs(determinant) < 0.0000001)
678         {
679             return null;
680         }
681         return new Point2d(l1p1x + (l1p2x * (l2p1x * l2p2y - l2p1y * l2p2x)) / determinant,
682                 l1p1y + (l1p2y * (l2p1x * l2p2y - l2p1y * l2p2x)) / determinant);
683     }
684 
685     /**
686      * Helper method for fractional projection which returns an offset line to the left of a segment at a distance of 1.
687      * @param segment segment number
688      * @return parallel line to the left of a segment at a distance of 1
689      */
690     private synchronized PolyLine2d unitOffsetSegment(final int segment)
691     {
692         return new PolyLine2d(get(segment), get(segment + 1)).offsetLine(1.0);
693     }
694 
695     /**
696      * Returns the projected directional radius of the line at a given fraction. Negative values reflect right-hand curvature in
697      * the design-line direction. The radius is taken as the minimum of the radii at the vertices before and after the given
698      * fraction. The radius at a vertex is calculated as the radius of a circle that is equidistant from both edges connected to
699      * the vertex. The circle center is on a line perpendicular to the shortest edge, crossing through the middle of the
700      * shortest edge. This method ignores Z components.
701      * @param fraction fraction along the line, between 0.0 and 1.0 (both inclusive)
702      * @return radius; the local radius; or si field set to Double.NaN if line is totally straight
703      * @throws IllegalArgumentException fraction out of bounds
704      */
705     public synchronized Length getProjectedRadius(final double fraction) throws IllegalArgumentException
706     {
707         Throw.when(fraction < 0.0 || fraction > 1.0, IllegalArgumentException.class,
708                 "Fraction %f is out of bounds [0.0 ... 1.0]", fraction);
709         if (this.vertexRadii == null)
710         {
711             this.vertexRadii = new Length[size() - 1];
712         }
713         int index = find(fraction * getLength());
714         if (index > 0 && this.vertexRadii[index] == null)
715         {
716             this.vertexRadii[index] = getProjectedVertexRadius(index);
717         }
718         if (index < size() - 2 && this.vertexRadii[index + 1] == null)
719         {
720             this.vertexRadii[index + 1] = getProjectedVertexRadius(index + 1);
721         }
722         if (index == 0)
723         {
724             if (this.vertexRadii.length < 2)
725             {
726                 return Length.ofSI(Double.NaN);
727             }
728             return this.vertexRadii[1];
729         }
730         if (index == size() - 2)
731         {
732             return this.vertexRadii[size() - 2];
733         }
734         return Math.abs(this.vertexRadii[index].si) < Math.abs(this.vertexRadii[index + 1].si) ? this.vertexRadii[index]
735                 : this.vertexRadii[index + 1];
736     }
737 
738     /**
739      * Calculates the directional radius at a vertex. Negative values reflect right-hand curvature in the design-line direction.
740      * The radius at a vertex is calculated as the radius of a circle that is equidistant from both edges connected to the
741      * vertex. The circle center is on a line perpendicular to the shortest edge, crossing through the middle of the shortest
742      * edge. This function ignores Z components.
743      * @param index index of the vertex in range [1 ... size() - 2]
744      * @return radius at the vertex
745      * @throws IndexOutOfBoundsException if the index is out of bounds
746      */
747     public synchronized Length getProjectedVertexRadius(final int index) throws IndexOutOfBoundsException
748     {
749         Throw.when(index < 1 || index > size() - 2, IndexOutOfBoundsException.class,
750                 "Index %d is out of bounds [1 ... size() - 2].", index);
751         determineFractionalHelpers(null, null);
752         double length1 = lengthAtIndex(index) - lengthAtIndex(index - 1);
753         double length2 = lengthAtIndex(index + 1) - lengthAtIndex(index);
754         int shortIndex = length1 < length2 ? index : index + 1;
755         // center of shortest edge
756         Point2d p1 =
757                 new Point2d(.5 * (get(shortIndex - 1).x + get(shortIndex).x), .5 * (get(shortIndex - 1).y + get(shortIndex).y));
758         // perpendicular to shortest edge, line crossing p1
759         Point2d p2 = new Point2d(p1.x + (get(shortIndex).y - get(shortIndex - 1).y),
760                 p1.y - (get(shortIndex).x - get(shortIndex - 1).x));
761         // vertex
762         Point2d p3 = get(index);
763         // point on line that splits angle between edges at vertex 50-50
764         Point2d p4 = this.fractionalHelperCenters[index];
765         if (p4 == null)
766         {
767             // parallel helper lines
768             p4 = new Point2d(p3.x + this.fractionalHelperDirections[index].x, p3.y + this.fractionalHelperDirections[index].y);
769         }
770         Point2d intersection = intersectionOfLines(p1, p2, p3, p4);
771         if (null == intersection)
772         {
773             return Length.ofSI(Double.NaN);
774         }
775         // determine left or right
776         double refLength = length1 < length2 ? length1 : length2;
777         double radius = intersection.distance(p1);
778         double i2p2 = intersection.distance(p2);
779         if (radius < i2p2 && i2p2 > refLength)
780         {
781             // left as p1 is closer than p2 (which was placed to the right) and not on the perpendicular line
782             return Length.ofSI(radius);
783         }
784         // right as not left
785         return Length.ofSI(-radius);
786     }
787 
788     /**
789      * Returns the length fraction at the vertex.
790      * @param index index of vertex [0 ... size() - 1]
791      * @return length fraction at the vertex
792      * @throws IndexOutOfBoundsException if the index is out of bounds
793      */
794     public double getVertexFraction(final int index) throws IndexOutOfBoundsException
795     {
796         Throw.when(index < 0 || index > size() - 1, IndexOutOfBoundsException.class, "Index %d is out of bounds [0 %d].", index,
797                 size() - 1);
798         return lengthAtIndex(index) / this.length.si;
799     }
800 
801     /**
802      * Retrieve the centroid of this OtsLine2d.
803      * @return the centroid of this OtsLine2d
804      */
805     public Point2d getCentroid()
806     {
807         return getAbsoluteBounds().midPoint();
808     }
809 
810     @Override
811     public Bounds<?, ?> getRelativeBounds()
812     {
813         return OtsShape.toRelativeTransform(getLocation()).transform(getAbsoluteBounds());
814     }
815 
816     @Override
817     @SuppressWarnings("checkstyle:designforextension")
818     public Point2d getLocation()
819     {
820         return getCentroid();
821     }
822 
823 }