1 package org.opentrafficsim.core.geometry;
2
3 import org.djunits.value.vdouble.scalar.Angle;
4 import org.djunits.value.vdouble.scalar.Direction;
5 import org.djutils.draw.line.PolyLine2d;
6 import org.djutils.draw.point.OrientedPoint2d;
7 import org.djutils.draw.point.Point2d;
8 import org.djutils.exceptions.Throw;
9 import org.djutils.exceptions.Try;
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35 public class ContinuousClothoid implements ContinuousLine
36 {
37
38
39 private static final double ANGLE_TOLERANCE = 2.0 * Math.PI / 3600.0;
40
41
42 private static final double SECANT_TOLERANCE = 1e-8;
43
44
45 private final OrientedPoint2d startPoint;
46
47
48 private final OrientedPoint2d endPoint;
49
50
51 private final double startCurvature;
52
53
54 private final double endCurvature;
55
56
57 private final double length;
58
59
60
61
62
63 private final double a;
64
65
66 private final double alphaMin;
67
68
69 private final double alphaMax;
70
71
72 private final double[] t0;
73
74
75 private final double[] n0;
76
77
78 private final boolean opposite;
79
80
81 private final boolean reflected;
82
83
84 private final ContinuousStraight straight;
85
86
87 private final ContinuousArc arc;
88
89
90 private boolean shiftDetermined;
91
92
93 private double shiftX;
94
95
96 private double shiftY;
97
98
99 private double dShiftX;
100
101
102 private double dShiftY;
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124 public ContinuousClothoid(final OrientedPoint2d startPoint, final OrientedPoint2d endPoint)
125 {
126 Throw.whenNull(startPoint, "Start point may not be null.");
127 Throw.whenNull(endPoint, "End point may not be null.");
128 this.startPoint = startPoint;
129 this.endPoint = endPoint;
130
131 double dx = endPoint.x - startPoint.x;
132 double dy = endPoint.y - startPoint.y;
133 double d2 = Math.hypot(dx, dy);
134 double d = Math.atan2(dy, dx);
135
136 double phi1 = normalizeAngle(d - startPoint.dirZ);
137 double phi2 = normalizeAngle(endPoint.dirZ - d);
138 double phi1Abs = Math.abs(phi1);
139 double phi2Abs = Math.abs(phi2);
140
141 if (phi1Abs < ANGLE_TOLERANCE && phi2Abs < ANGLE_TOLERANCE)
142 {
143
144 this.length = Math.hypot(endPoint.x - startPoint.x, endPoint.y - startPoint.y);
145 this.a = Double.POSITIVE_INFINITY;
146 this.startCurvature = 0.0;
147 this.endCurvature = 0.0;
148 this.straight = new ContinuousStraight(startPoint, this.length);
149 this.arc = null;
150 this.alphaMin = 0.0;
151 this.alphaMax = 0.0;
152 this.t0 = null;
153 this.n0 = null;
154 this.opposite = false;
155 this.reflected = false;
156 return;
157 }
158 else if (Math.abs(phi2 - phi1) < ANGLE_TOLERANCE)
159 {
160
161 double r = .5 * d2 / Math.sin(phi1);
162 double cosStartDirection = Math.cos(startPoint.dirZ);
163 double sinStartDirection = Math.sin(startPoint.dirZ);
164 double ang = Math.PI / 2.0;
165 double cosAng = Math.cos(ang);
166 double sinAng = Math.sin(ang);
167 double x0 = startPoint.x - r * (cosStartDirection * cosAng + sinStartDirection * sinAng);
168 double y0 = startPoint.y - r * (cosStartDirection * -sinAng + sinStartDirection * cosAng);
169 double from = Math.atan2(startPoint.y - y0, startPoint.x - x0);
170 double to = Math.atan2(endPoint.y - y0, endPoint.x - x0);
171 if (r < 0 && to > from)
172 {
173 to = to - 2.0 * Math.PI;
174 }
175 else if (r > 0 && to < from)
176 {
177 to = to + 2.0 * Math.PI;
178 }
179 Angle angle = Angle.instantiateSI(Math.abs(to - from));
180 this.length = angle.si * Math.abs(r);
181 this.a = 0.0;
182 this.startCurvature = 1.0 / r;
183 this.endCurvature = 1.0 / r;
184 this.straight = null;
185 this.arc = new ContinuousArc(startPoint, Math.abs(r), r > 0.0, angle);
186 this.alphaMin = 0.0;
187 this.alphaMax = 0.0;
188 this.t0 = null;
189 this.n0 = null;
190 this.opposite = false;
191 this.reflected = false;
192 return;
193 }
194 this.straight = null;
195 this.arc = null;
196
197
198
199 if (phi2Abs < phi1Abs)
200 {
201 this.opposite = true;
202 double phi3 = phi1;
203 phi1 = -phi2;
204 phi2 = -phi3;
205 dx = -dx;
206 dy = -dy;
207 }
208 else
209 {
210 this.opposite = false;
211 }
212
213
214 this.reflected = phi2 < 0 || phi2 > Math.PI;
215 if (this.reflected)
216 {
217 phi1 = -phi1;
218 phi2 = -phi2;
219 }
220
221
222 double[] cs = Fresnel.fresnel(alphaToT(phi1 + phi2));
223 double h = cs[1] * Math.cos(phi1) - cs[0] * Math.sin(phi1);
224 boolean cShape = 0 < phi1 && phi1 < phi2 && phi2 < Math.PI && h < 0;
225 double theta = getTheta(phi1, phi2, cShape);
226 double aSign = cShape ? -1.0 : 1.0;
227 double thetaSign = -aSign;
228
229 double v1 = theta + phi1 + phi2;
230 double v2 = theta + phi1;
231 double[] cs0 = Fresnel.fresnel(alphaToT(theta));
232 double[] cs1 = Fresnel.fresnel(alphaToT(v1));
233 this.a = d2 / ((cs1[1] + aSign * cs0[1]) * Math.sin(v2) + (cs1[0] + aSign * cs0[0]) * Math.cos(v2));
234
235 dx /= d2;
236 dy /= d2;
237 if (this.reflected)
238 {
239
240 this.t0 = new double[] {Math.cos(-v2) * dx + Math.sin(-v2) * dy, -Math.sin(-v2) * dx + Math.cos(-v2) * dy};
241 this.n0 = new double[] {-this.t0[1], this.t0[0]};
242 }
243 else
244 {
245 this.t0 = new double[] {Math.cos(v2) * dx + Math.sin(v2) * dy, -Math.sin(v2) * dx + Math.cos(v2) * dy};
246 this.n0 = new double[] {this.t0[1], -this.t0[0]};
247 }
248
249 this.alphaMin = thetaSign * theta;
250 this.alphaMax = v1;
251 double sign = (this.reflected ? -1.0 : 1.0);
252 double curveMin = Math.PI * alphaToT(this.alphaMin) / this.a;
253 double curveMax = Math.PI * alphaToT(v1) / this.a;
254 this.startCurvature = sign * (this.opposite ? -curveMax : curveMin);
255 this.endCurvature = sign * (this.opposite ? -curveMin : curveMax);
256 this.length = this.a * (alphaToT(v1) - alphaToT(this.alphaMin));
257 }
258
259
260
261
262
263
264
265
266 public ContinuousClothoid(final OrientedPoint2d startPoint, final double a, final double startCurvature,
267 final double endCurvature)
268 {
269 Throw.whenNull(startPoint, "Start point may not be null.");
270 Throw.when(a <= 0.0, IllegalArgumentException.class, "A value must be above 0.");
271 this.startPoint = startPoint;
272
273 this.a = a * Math.sqrt(Math.PI);
274 this.length = a * a * Math.abs(endCurvature - startCurvature);
275 this.startCurvature = startCurvature;
276 this.endCurvature = endCurvature;
277
278 double l1 = a * a * startCurvature;
279 double l2 = a * a * endCurvature;
280 this.alphaMin = Math.abs(l1) * startCurvature / 2.0;
281 this.alphaMax = Math.abs(l2) * endCurvature / 2.0;
282
283 double ang = normalizeAngle(startPoint.dirZ) - Math.abs(this.alphaMin);
284 this.t0 = new double[] {Math.cos(ang), Math.sin(ang)};
285 this.n0 = new double[] {this.t0[1], -this.t0[0]};
286 Direction endDirection = Direction.instantiateSI(ang + Math.abs(this.alphaMax));
287 if (startCurvature > endCurvature)
288 {
289
290
291 double m = Math.tan(startPoint.dirZ + Math.PI / 2.0);
292
293
294 double onePlusMm = 1.0 + m * m;
295 double oneMinusMm = 1.0 - m * m;
296 double mmMinusOne = m * m - 1.0;
297 double twoM = 2.0 * m;
298 double t00 = this.t0[0];
299 double t01 = this.t0[1];
300 double n00 = this.n0[0];
301 double n01 = this.n0[1];
302 this.t0[0] = (oneMinusMm * t00 + 2 * m * t01) / onePlusMm;
303 this.t0[1] = (twoM * t00 + mmMinusOne * t01) / onePlusMm;
304 this.n0[0] = (oneMinusMm * n00 + 2 * m * n01) / onePlusMm;
305 this.n0[1] = (twoM * n00 + mmMinusOne * n01) / onePlusMm;
306
307 double ang2 = Math.atan2(this.t0[1], this.t0[0]);
308 endDirection = Direction.instantiateSI(ang2 - Math.abs(this.alphaMax) + Math.PI);
309 }
310 PolyLine2d line = flatten(new Flattener.NumSegments(1));
311 Point2d end = Try.assign(() -> line.get(line.size() - 1), "Line does not have an end point.");
312 this.endPoint = new OrientedPoint2d(end.x, end.y, endDirection.si);
313
314
315 this.straight = null;
316 this.arc = null;
317 this.opposite = false;
318 this.reflected = false;
319 }
320
321
322
323
324
325
326
327
328
329
330
331 public static ContinuousClothoid withLength(final OrientedPoint2d startPoint, final double length,
332 final double startCurvature, final double endCurvature)
333 {
334 Throw.when(length <= 0.0, IllegalArgumentException.class, "Length must be above 0.");
335 double a = Math.sqrt(length / Math.abs(endCurvature - startCurvature));
336 return new ContinuousClothoid(startPoint, a, startCurvature, endCurvature);
337 }
338
339
340
341
342
343
344 private static double normalizeAngle(final double angle)
345 {
346 double out = angle;
347 while (out > Math.PI)
348 {
349 out -= 2 * Math.PI;
350 }
351 while (out < -Math.PI)
352 {
353 out += 2 * Math.PI;
354 }
355 return out;
356 }
357
358
359
360
361
362
363 private static double alphaToT(final double alpha)
364 {
365 return alpha >= 0 ? Math.sqrt(alpha * 2.0 / Math.PI) : -Math.sqrt(-alpha * 2.0 / Math.PI);
366 }
367
368
369
370
371
372
373
374
375 private static double getTheta(final double phi1, final double phi2, final boolean cShape)
376 {
377 double sign, phiMin, phiMax;
378 if (cShape)
379 {
380 double lambda = (1 - Math.cos(phi1)) / (1 - Math.cos(phi2));
381 phiMin = 0.0;
382 phiMax = (lambda * lambda * (phi1 + phi2)) / (1 - (lambda * lambda));
383 sign = -1.0;
384 }
385 else
386 {
387 phiMin = Math.max(0, -phi1);
388 phiMax = Math.PI / 2 - phi1;
389 sign = 1;
390 }
391
392 double fMin = fTheta(phiMin, phi1, phi2, sign);
393 double fMax = fTheta(phiMax, phi1, phi2, sign);
394 if (fMin * fMax > 0)
395 {
396 throw new RuntimeException("f(phiMin) and f(phiMax) have the same sign, we cant find f(theta) = 0 between them.");
397 }
398
399
400 double x0 = phiMin;
401 double x1 = phiMax;
402 double x2 = 0;
403 for (int i = 0; i < 100; i++)
404 {
405 double f1 = fTheta(x1, phi1, phi2, sign);
406 x2 = x1 - f1 * (x1 - x0) / (f1 - fTheta(x0, phi1, phi2, sign));
407 x2 = Math.max(Math.min(x2, phiMax), phiMin);
408 x0 = x1;
409 x1 = x2;
410 if (Math.abs(x0 - x1) < SECANT_TOLERANCE || Math.abs(x0 / x1 - 1) < SECANT_TOLERANCE
411 || Math.abs(f1) < SECANT_TOLERANCE)
412 {
413 return x2;
414 }
415 }
416
417 return x2;
418 }
419
420
421
422
423
424
425
426
427
428
429
430 private static double fTheta(final double theta, final double phi1, final double phi2, final double sign)
431 {
432 double thetaPhi1 = theta + phi1;
433 double[] cs0 = Fresnel.fresnel(alphaToT(theta));
434 double[] cs1 = Fresnel.fresnel(alphaToT(thetaPhi1 + phi2));
435 return (cs1[1] + sign * cs0[1]) * Math.cos(thetaPhi1) - (cs1[0] + sign * cs0[0]) * Math.sin(thetaPhi1);
436 }
437
438
439 @Override
440 public OrientedPoint2d getStartPoint()
441 {
442 return this.startPoint;
443 }
444
445
446 @Override
447 public OrientedPoint2d getEndPoint()
448 {
449 return this.endPoint;
450 }
451
452
453 @Override
454 public double getStartCurvature()
455 {
456 return this.startCurvature;
457 }
458
459
460 @Override
461 public double getEndCurvature()
462 {
463 return this.endCurvature;
464 }
465
466
467 @Override
468 public double getStartRadius()
469 {
470 return 1.0 / this.startCurvature;
471 }
472
473
474 @Override
475 public double getEndRadius()
476 {
477 return 1.0 / this.endCurvature;
478 }
479
480
481
482
483
484 public double getA()
485 {
486
487
488 return this.a / Math.sqrt(Math.PI);
489 }
490
491
492
493
494 private void assureShift()
495 {
496 if (this.shiftDetermined)
497 {
498 return;
499 }
500
501 OrientedPoint2d p1 = this.opposite ? this.endPoint : this.startPoint;
502 OrientedPoint2d p2 = this.opposite ? this.startPoint : this.endPoint;
503
504
505 double[] csMin = Fresnel.fresnel(alphaToT(this.alphaMin));
506 double xMin = this.a * (csMin[0] * this.t0[0] - csMin[1] * this.n0[0]);
507 double yMin = this.a * (csMin[0] * this.t0[1] - csMin[1] * this.n0[1]);
508 this.shiftX = p1.x - xMin;
509 this.shiftY = p1.y - yMin;
510
511
512 if (p2 != null)
513 {
514 double[] csMax = Fresnel.fresnel(alphaToT(this.alphaMax));
515 double xMax = this.a * (csMax[0] * this.t0[0] - csMax[1] * this.n0[0]);
516 double yMax = this.a * (csMax[0] * this.t0[1] - csMax[1] * this.n0[1]);
517 this.dShiftX = p2.x - (xMax + this.shiftX);
518 this.dShiftY = p2.y - (yMax + this.shiftY);
519 }
520 else
521 {
522 this.dShiftX = 0.0;
523 this.dShiftY = 0.0;
524 }
525
526 this.shiftDetermined = true;
527 }
528
529
530
531
532
533
534
535 private Point2d getPoint(final double fraction, final double offset)
536 {
537 double f = this.opposite ? 1.0 - fraction : fraction;
538 double alpha = this.alphaMin + f * (this.alphaMax - this.alphaMin);
539 double[] cs = Fresnel.fresnel(alphaToT(alpha));
540 double x = this.shiftX + this.a * (cs[0] * this.t0[0] - cs[1] * this.n0[0]) + f * this.dShiftX;
541 double y = this.shiftY + this.a * (cs[0] * this.t0[1] - cs[1] * this.n0[1]) + f * this.dShiftY;
542 double d = getDirection(alpha) + Math.PI / 2;
543 return new Point2d(x + Math.cos(d) * offset, y + Math.sin(d) * offset);
544 }
545
546
547
548
549
550
551 private double getDirection(final double alpha)
552 {
553 double rot = Math.atan2(this.t0[1], this.t0[0]);
554
555 rot += this.reflected ? -Math.abs(alpha) : Math.abs(alpha);
556 if (this.opposite)
557 {
558 rot += Math.PI;
559 }
560 return normalizeAngle(rot);
561 }
562
563
564 @Override
565 public PolyLine2d flatten(final Flattener flattener)
566 {
567 Throw.whenNull(flattener, "Flattener may not be null.");
568 if (this.straight != null)
569 {
570 return this.straight.flatten(flattener);
571 }
572 if (this.arc != null)
573 {
574 return this.arc.flatten(flattener);
575 }
576 assureShift();
577 return flattener.flatten(new FlattableLine()
578 {
579
580 @Override
581 public Point2d get(final double fraction)
582 {
583 return getPoint(fraction, 0.0);
584 }
585
586
587 @Override
588 public double getDirection(final double fraction)
589 {
590 return ContinuousClothoid.this.getDirection(ContinuousClothoid.this.alphaMin
591 + fraction * (ContinuousClothoid.this.alphaMax - ContinuousClothoid.this.alphaMin));
592 }
593 });
594 }
595
596
597 @Override
598 public PolyLine2d flattenOffset(final FractionalLengthData offsets, final Flattener flattener)
599 {
600 Throw.whenNull(offsets, "Offsets may not be null.");
601 Throw.whenNull(flattener, "Flattener may not be null.");
602 if (this.straight != null)
603 {
604 return this.straight.flattenOffset(offsets, flattener);
605 }
606 if (this.arc != null)
607 {
608 return this.arc.flattenOffset(offsets, flattener);
609 }
610 assureShift();
611 return flattener.flatten(new FlattableLine()
612 {
613
614 @Override
615 public Point2d get(final double fraction)
616 {
617 return getPoint(fraction, offsets.get(fraction));
618 }
619
620
621 @Override
622 public double getDirection(final double fraction)
623 {
624 return ContinuousClothoid.this.getDirection(ContinuousClothoid.this.alphaMin
625 + fraction * (ContinuousClothoid.this.alphaMax - ContinuousClothoid.this.alphaMin));
626 }
627 });
628 }
629
630
631 @Override
632 public double getLength()
633 {
634 return this.length;
635 }
636
637
638
639
640
641
642 public String getAppliedShape()
643 {
644 return this.straight == null ? (this.arc == null ? "Clothoid" : "Arc") : "Straight";
645 }
646
647
648 @Override
649 public String toString()
650 {
651 return "ContinuousClothoid [startPoint=" + this.startPoint + ", endPoint=" + this.endPoint + ", startCurvature="
652 + this.startCurvature + ", endCurvature=" + this.endCurvature + ", length=" + this.length + "]";
653 }
654
655 }