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