1 package org.opentrafficsim.core.geometry;
2
3 import java.awt.geom.Line2D;
4 import java.util.ArrayList;
5 import java.util.List;
6
7 import org.locationtech.jts.geom.Coordinate;
8 import org.locationtech.jts.geom.Geometry;
9 import org.locationtech.jts.linearref.LengthIndexedLine;
10 import org.locationtech.jts.operation.buffer.BufferParameters;
11 import org.opentrafficsim.base.logger.Cat;
12 import org.opentrafficsim.core.network.NetworkException;
13
14 import nl.tudelft.simulation.dsol.logger.SimLogger;
15
16
17
18
19
20
21
22
23
24
25
26 public final class OTSBufferingJTS
27 {
28
29 private static final int QUADRANTSEGMENTS = 16;
30
31
32
33
34 private OTSBufferingJTS()
35 {
36
37 }
38
39
40
41
42
43
44 private static double norm(final double angle)
45 {
46 double normalized = angle % (2 * Math.PI);
47 if (normalized < 0.0)
48 {
49 normalized += 2 * Math.PI;
50 }
51 return normalized;
52 }
53
54
55
56
57
58
59 private static double angle(final Coordinate c1, final Coordinate c2)
60 {
61 return norm(Math.atan2(c2.y - c1.y, c2.x - c1.x));
62 }
63
64
65
66
67
68
69
70
71
72
73
74
75 public static double distanceLineSegmentToPoint(final OTSPoint3D lineP1, final OTSPoint3D lineP2, final OTSPoint3D point)
76 {
77 return closestPointOnSegmentToPoint(lineP1, lineP2, point).distanceSI(point);
78 }
79
80
81
82
83
84
85
86
87
88
89
90
91 public static OTSPoint3D closestPointOnSegmentToPoint(final OTSPoint3D lineP1, final OTSPoint3D lineP2,
92 final OTSPoint3D point)
93 {
94 double dX = lineP2.x - lineP1.x;
95 double dY = lineP2.y - lineP1.y;
96 if ((0 == dX) && (0 == dY))
97 {
98 return lineP1;
99 }
100 final double u = ((point.x - lineP1.x) * dX + (point.y - lineP1.y) * dY) / (dX * dX + dY * dY);
101 if (u < 0)
102 {
103 return lineP1;
104 }
105 else if (u > 1)
106 {
107 return lineP2;
108 }
109 else
110 {
111 return new OTSPoint3D(lineP1.x + u * dX, lineP1.y + u * dY);
112 }
113 }
114
115
116
117
118
119
120
121 public static OTSLine3D offsetLine(final OTSLine3D referenceLine, final double offset)
122 {
123 try
124 {
125 double bufferOffset = Math.abs(offset);
126 final double precision = 0.00001;
127 if (bufferOffset < precision)
128 {
129 return referenceLine;
130 }
131 final double circlePrecision = 0.001;
132 List<OTSPoint3D> points = new ArrayList<>();
133
134 OTSPoint3D prevPoint = referenceLine.get(0);
135 Double prevAngle = null;
136 for (int index = 0; index < referenceLine.size() - 1; index++)
137 {
138 OTSPoint3D nextPoint = referenceLine.get(index + 1);
139 double angle = Math.atan2(nextPoint.y - prevPoint.y, nextPoint.x - prevPoint.x);
140 OTSPoint3D segmentFrom =
141 new OTSPoint3D(prevPoint.x - Math.sin(angle) * offset, prevPoint.y + Math.cos(angle) * offset);
142 OTSPoint3D segmentTo =
143 new OTSPoint3D(nextPoint.x - Math.sin(angle) * offset, nextPoint.y + Math.cos(angle) * offset);
144 if (index > 0)
145 {
146 double deltaAngle = angle - prevAngle;
147 if (Math.abs(deltaAngle) > Math.PI)
148 {
149 deltaAngle -= Math.signum(deltaAngle) * 2 * Math.PI;
150 }
151 if (deltaAngle * offset > 0)
152 {
153
154
155 OTSPoint3D pPoint = null;
156 for (int i = 0; i < points.size(); i++)
157 {
158 OTSPoint3D p = points.get(i);
159 if (Double.isNaN(p.z))
160 {
161 continue;
162 }
163 if (null != pPoint)
164 {
165 double pAngle = Math.atan2(p.y - pPoint.y, p.x - pPoint.x);
166 double totalAngle = angle - pAngle;
167 if (Math.abs(totalAngle) > Math.PI)
168 {
169 totalAngle += Math.signum(totalAngle) * 2 * Math.PI;
170 }
171 if (Math.abs(totalAngle) > 0.01)
172 {
173
174
175 OTSPoint3D intermediatePoint =
176 intersectionOfLineSegments(pPoint, p, segmentFrom, segmentTo);
177 if (null != intermediatePoint)
178 {
179
180 intermediatePoint =
181 new OTSPoint3D(intermediatePoint.x, intermediatePoint.y, Double.NaN);
182
183
184 points.add(intermediatePoint);
185 }
186 }
187 }
188 pPoint = p;
189 }
190 }
191 else
192 {
193
194
195
196 int numSegments = 1;
197 if (Math.abs(deltaAngle) > Math.PI / 2)
198 {
199 numSegments = 2;
200 }
201 for (; numSegments < 1000; numSegments *= 2)
202 {
203 double maxError = bufferOffset * (1 - Math.abs(Math.cos(deltaAngle / numSegments / 2)));
204 if (maxError < circlePrecision)
205 {
206 break;
207 }
208 }
209
210 for (int additionalPoint = 1; additionalPoint < numSegments; additionalPoint++)
211 {
212 double intermediateAngle =
213 (additionalPoint * angle + (numSegments - additionalPoint) * prevAngle) / numSegments;
214 if (prevAngle * angle < 0 && Math.abs(prevAngle) > Math.PI / 2 && Math.abs(angle) > Math.PI / 2)
215 {
216 intermediateAngle += Math.PI;
217 }
218 OTSPoint3D intermediatePoint = new OTSPoint3D(prevPoint.x - Math.sin(intermediateAngle) * offset,
219 prevPoint.y + Math.cos(intermediateAngle) * offset);
220
221
222 points.add(intermediatePoint);
223 }
224 }
225 }
226 points.add(segmentFrom);
227 points.add(segmentTo);
228 prevPoint = nextPoint;
229 prevAngle = angle;
230 }
231
232
233
234 for (int index = 1; index < points.size() - 1; index++)
235 {
236 OTSPoint3D checkPoint = points.get(index);
237 prevPoint = null;
238 boolean tooClose = false;
239 boolean somewhereAtCorrectDistance = false;
240 for (int i = 0; i < referenceLine.size(); i++)
241 {
242 OTSPoint3D p = referenceLine.get(i);
243 if (null != prevPoint)
244 {
245 OTSPoint3D closestPoint = closestPointOnSegmentToPoint(prevPoint, p, checkPoint);
246 if (closestPoint != referenceLine.get(0) && closestPoint != referenceLine.get(referenceLine.size() - 1))
247 {
248 double distance = closestPoint.horizontalDistanceSI(checkPoint);
249 if (distance < bufferOffset - circlePrecision)
250 {
251
252
253 tooClose = true;
254 break;
255 }
256 else if (distance < bufferOffset + precision)
257 {
258 somewhereAtCorrectDistance = true;
259 }
260 }
261 }
262 prevPoint = p;
263 }
264 if (tooClose || !somewhereAtCorrectDistance)
265 {
266
267 points.remove(index);
268 index--;
269 }
270 }
271
272 for (int index = 0; index < points.size(); index++)
273 {
274 OTSPoint3D p = points.get(index);
275 if (Double.isNaN(p.z))
276 {
277 points.set(index, new OTSPoint3D(p.x, p.y, 0));
278 }
279 }
280 return OTSLine3D.createAndCleanOTSLine3D(points);
281 }
282 catch (OTSGeometryException exception)
283 {
284 SimLogger.always().error(exception, "Exception in offsetLine - should never happen");
285 return null;
286 }
287 }
288
289
290
291
292
293
294
295
296
297 private static OTSPoint3D intersectionOfLineSegments(final OTSPoint3D line1P1, final OTSPoint3D line1P2,
298 final OTSPoint3D line2P1, final OTSPoint3D line2P2)
299 {
300 double denominator =
301 (line2P2.y - line2P1.y) * (line1P2.x - line1P1.x) - (line2P2.x - line2P1.x) * (line1P2.y - line1P1.y);
302 if (denominator == 0f)
303 {
304 return null;
305 }
306 double uA = ((line2P2.x - line2P1.x) * (line1P1.y - line2P1.y) - (line2P2.y - line2P1.y) * (line1P1.x - line2P1.x))
307 / denominator;
308 if ((uA < 0f) || (uA > 1f))
309 {
310 return null;
311 }
312 double uB = ((line1P2.x - line1P1.x) * (line1P1.y - line2P1.y) - (line1P2.y - line1P1.y) * (line1P1.x - line2P1.x))
313 / denominator;
314 if (uB < 0 || uB > 1)
315 {
316 return null;
317 }
318 return new OTSPoint3D(line1P1.x + uA * (line1P2.x - line1P1.x), line1P1.y + uA * (line1P2.y - line1P1.y), 0);
319 }
320
321
322
323
324
325
326
327
328 @SuppressWarnings("checkstyle:methodlength")
329 public static OTSLine3D offsetGeometryOLD(final OTSLine3D referenceLine, final double offset) throws OTSGeometryException
330 {
331 Coordinate[] referenceCoordinates = referenceLine.getCoordinates();
332
333 double bufferOffset = Math.abs(offset);
334 final double precision = 0.000001;
335 if (bufferOffset < precision)
336 {
337
338 return new OTSLine3D(referenceCoordinates);
339 }
340 Geometry geometryLine = referenceLine.getLineString();
341 Coordinate[] bufferCoordinates =
342 geometryLine.buffer(bufferOffset, QUADRANTSEGMENTS, BufferParameters.CAP_FLAT).getCoordinates();
343
344
345
346
347
348 Coordinate sC0 = referenceCoordinates[0];
349 Coordinate sC1 = referenceCoordinates[1];
350 Coordinate eCm1 = referenceCoordinates[referenceCoordinates.length - 1];
351 Coordinate eCm2 = referenceCoordinates[referenceCoordinates.length - 2];
352
353 double expectedStartAngle = norm(angle(sC0, sC1) + Math.signum(offset) * Math.PI / 2.0);
354 double expectedEndAngle = norm(angle(eCm2, eCm1) + Math.signum(offset) * Math.PI / 2.0);
355 Coordinate sExpected = new Coordinate(sC0.x + bufferOffset * Math.cos(expectedStartAngle),
356 sC0.y + bufferOffset * Math.sin(expectedStartAngle));
357 Coordinate eExpected = new Coordinate(eCm1.x + bufferOffset * Math.cos(expectedEndAngle),
358 eCm1.y + bufferOffset * Math.sin(expectedEndAngle));
359
360
361 double dS = Double.MAX_VALUE;
362 double dE = Double.MAX_VALUE;
363 int sIndex = -1;
364 int eIndex = -1;
365 for (int i = 0; i < bufferCoordinates.length; i++)
366 {
367 Coordinate c = bufferCoordinates[i];
368 double dsc = c.distance(sExpected);
369 double dec = c.distance(eExpected);
370 if (dsc < dS)
371 {
372 dS = dsc;
373 sIndex = i;
374 }
375 if (dec < dE)
376 {
377 dE = dec;
378 eIndex = i;
379 }
380 }
381
382 if (sIndex == -1)
383 {
384 throw new OTSGeometryException("offsetGeometry: startIndex not found for line " + referenceLine);
385 }
386 if (eIndex == -1)
387 {
388 throw new OTSGeometryException("offsetGeometry: endIndex not found for line " + referenceLine);
389 }
390 if (dS > 0.01)
391 {
392 SimLogger.filter(Cat.CORE).trace(referenceLine.toExcel() + "\n\n\n\n" + new OTSLine3D(bufferCoordinates).toExcel()
393 + "\n\n\n\n" + sExpected + "\n" + eExpected);
394 throw new OTSGeometryException("offsetGeometry: startDistance too big (" + dS + ") for line " + referenceLine);
395 }
396 if (dE > 0.01)
397 {
398 throw new OTSGeometryException("offsetGeometry: endDistance too big (" + dE + ") for line " + referenceLine);
399 }
400
401
402 boolean ok = true;
403 int i = sIndex;
404 Coordinate lastC = null;
405 List<OTSPoint3D> result = new ArrayList<>();
406 while (ok)
407 {
408 Coordinate c = bufferCoordinates[i];
409 if (lastC != null && close(c, lastC, sC0, eCm1))
410 {
411 ok = false;
412 break;
413 }
414 result.add(new OTSPoint3D(c));
415 if (i == eIndex)
416 {
417 return OTSLine3D.createAndCleanOTSLine3D(result);
418 }
419 i = (i == bufferCoordinates.length - 1) ? 0 : i + 1;
420 lastC = c;
421 }
422
423
424 ok = true;
425 i = sIndex;
426 lastC = null;
427 result = new ArrayList<>();
428 while (ok)
429 {
430 Coordinate c = bufferCoordinates[i];
431 if (lastC != null && close(c, lastC, sC0, eCm1))
432 {
433 ok = false;
434 break;
435 }
436 result.add(new OTSPoint3D(c));
437 if (i == eIndex)
438 {
439 return OTSLine3D.createAndCleanOTSLine3D(result);
440 }
441 i = (i == 0) ? bufferCoordinates.length - 1 : i - 1;
442 lastC = c;
443 }
444
445
446
447 throw new OTSGeometryException("offsetGeometry: could not find offset in either direction for line " + referenceLine);
448 }
449
450
451
452
453
454
455
456
457 private static boolean close(final Coordinate lineC1, final Coordinate lineC2, final Coordinate... check)
458 {
459 Line2D.Double line = new Line2D.Double(lineC1.x, lineC1.y, lineC2.x, lineC2.y);
460 for (Coordinate c : check)
461 {
462 if (line.ptSegDist(c.x, c.y) < 0.01)
463 {
464 return true;
465 }
466 }
467 return false;
468 }
469
470
471
472
473
474
475
476
477
478
479 public static OTSLine3D offsetLine(final OTSLine3D referenceLine, final double offsetAtStart, final double offsetAtEnd)
480 throws OTSGeometryException
481 {
482
483
484
485 OTSLine3D offsetLineAtStart = offsetLine(referenceLine, offsetAtStart);
486 if (offsetAtStart == offsetAtEnd)
487 {
488 return offsetLineAtStart;
489 }
490
491
492 OTSLine3D offsetLineAtEnd = offsetLine(referenceLine, offsetAtEnd);
493
494
495 Geometry startGeometry = offsetLineAtStart.getLineString();
496 Geometry endGeometry = offsetLineAtEnd.getLineString();
497 LengthIndexedLine first = new LengthIndexedLine(startGeometry);
498 double firstLength = startGeometry.getLength();
499 LengthIndexedLine second = new LengthIndexedLine(endGeometry);
500 double secondLength = endGeometry.getLength();
501 ArrayList<Coordinate> out = new ArrayList<Coordinate>();
502 Coordinate[] firstCoordinates = startGeometry.getCoordinates();
503 Coordinate[] secondCoordinates = endGeometry.getCoordinates();
504 int firstIndex = 0;
505 int secondIndex = 0;
506 Coordinate prevCoordinate = null;
507 final double tooClose = 0.05;
508 while (firstIndex < firstCoordinates.length && secondIndex < secondCoordinates.length)
509 {
510 double firstRatio = firstIndex < firstCoordinates.length ? first.indexOf(firstCoordinates[firstIndex]) / firstLength
511 : Double.MAX_VALUE;
512 double secondRatio = secondIndex < secondCoordinates.length
513 ? second.indexOf(secondCoordinates[secondIndex]) / secondLength : Double.MAX_VALUE;
514 double ratio;
515 if (firstRatio < secondRatio)
516 {
517 ratio = firstRatio;
518 firstIndex++;
519 }
520 else
521 {
522 ratio = secondRatio;
523 secondIndex++;
524 }
525 Coordinate firstCoordinate = first.extractPoint(ratio * firstLength);
526 Coordinate secondCoordinate = second.extractPoint(ratio * secondLength);
527 Coordinate resultCoordinate = new Coordinate((1 - ratio) * firstCoordinate.x + ratio * secondCoordinate.x,
528 (1 - ratio) * firstCoordinate.y + ratio * secondCoordinate.y);
529 if (null == prevCoordinate || resultCoordinate.distance(prevCoordinate) > tooClose)
530 {
531 out.add(resultCoordinate);
532 prevCoordinate = resultCoordinate;
533 }
534 }
535 Coordinate[] resultCoordinates = new Coordinate[out.size()];
536 for (int index = 0; index < out.size(); index++)
537 {
538 resultCoordinates[index] = out.get(index);
539 }
540 return new OTSLine3D(resultCoordinates);
541 }
542
543
544
545
546
547
548 public static void main(final String[] args) throws NetworkException, OTSGeometryException
549 {
550
551
552
553
554
555 OTSLine3D line = new OTSLine3D(new OTSPoint3D[] { new OTSPoint3D(-579.253, 60.157, 4.710),
556 new OTSPoint3D(-579.253, 60.144, 4.712), new OTSPoint3D(-579.253, 60.144, 0.000),
557 new OTSPoint3D(-579.251, 60.044, 0.000), new OTSPoint3D(-579.246, 59.944, 0.000),
558 new OTSPoint3D(-579.236, 59.845, 0.000), new OTSPoint3D(-579.223, 59.746, 0.000),
559 new OTSPoint3D(-579.206, 59.647, 0.000), new OTSPoint3D(-579.185, 59.549, 0.000),
560 new OTSPoint3D(-579.161, 59.452, 0.000), new OTSPoint3D(-579.133, 59.356, 0.000),
561 new OTSPoint3D(-579.101, 59.261, 0.000), new OTSPoint3D(-579.066, 59.168, 0.000),
562 new OTSPoint3D(-579.028, 59.075, 0.000), new OTSPoint3D(-578.986, 58.985, 0.000),
563 new OTSPoint3D(-578.940, 58.896, 0.000), new OTSPoint3D(-578.891, 58.809, 0.000),
564 new OTSPoint3D(-578.839, 58.723, 0.000), new OTSPoint3D(-578.784, 58.640, 0.000),
565 new OTSPoint3D(-578.725, 58.559, 0.000), new OTSPoint3D(-578.664, 58.480, 0.000),
566 new OTSPoint3D(-578.599, 58.403, 0.000), new OTSPoint3D(-578.532, 58.329, 0.000),
567 new OTSPoint3D(-578.462, 58.258, 0.000), new OTSPoint3D(-578.390, 58.189, 0.000),
568 new OTSPoint3D(-578.314, 58.123, 0.000), new OTSPoint3D(-578.237, 58.060, 0.000),
569 new OTSPoint3D(-578.157, 58.000, 0.000), new OTSPoint3D(-578.075, 57.943, 0.000),
570 new OTSPoint3D(-577.990, 57.889, 0.000), new OTSPoint3D(-577.904, 57.839, 0.000),
571 new OTSPoint3D(-577.816, 57.791, 0.000), new OTSPoint3D(-577.726, 57.747, 0.000),
572 new OTSPoint3D(-577.635, 57.707, 0.000), new OTSPoint3D(-577.542, 57.670, 0.000),
573 new OTSPoint3D(-577.448, 57.636, 0.000), new OTSPoint3D(-577.352, 57.606, 0.000),
574 new OTSPoint3D(-577.256, 57.580, 0.000), new OTSPoint3D(-577.159, 57.557, 0.000),
575 new OTSPoint3D(-577.060, 57.538, 0.000), new OTSPoint3D(-576.962, 57.523, 0.000),
576 new OTSPoint3D(-576.862, 57.512, 0.000), new OTSPoint3D(-576.763, 57.504, 0.000),
577 new OTSPoint3D(-576.663, 57.500, 0.000), new OTSPoint3D(-576.623, 57.500, 6.278),
578 new OTSPoint3D(-576.610, 57.500, 6.280), new OTSPoint3D(-567.499, 57.473, 6.280) });
579 System.out.println(line.toExcel());
580 System.out.println(OTSBufferingJTS.offsetGeometryOLD(line, -1.831));
581 }
582 }