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 OTSPoint3DOTSPoint3D.html#OTSPoint3D">OTSPoint3DOTSPoint3D.html#OTSPoint3D">OTSPoint3D lineP1, final OTSPoint3DOTSPoint3D.html#OTSPoint3D">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 OTSPoint3DOTSPoint3D.html#OTSPoint3D">OTSPoint3Dint3D">OTSPoint3D closestPointOnSegmentToPoint(final OTSPoint3DOTSPoint3D.html#OTSPoint3D">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 OTSLine3DSLine3D.html#OTSLine3D">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.html#OTSPoint3D">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 OTSPoint3DTSPoint3D.html#OTSPoint3D">OTSPoint3DPoint3D">OTSPoint3D intersectionOfLineSegments(final OTSPoint3DTSPoint3D.html#OTSPoint3D">OTSPoint3D line1P1, final OTSPoint3D line1P2,
298             final OTSPoint3DTSPoint3D.html#OTSPoint3D">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.html#OTSLine3D">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 OTSLine3DSLine3D.html#OTSLine3D">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         OTSLine3Dtry/OTSLine3D.html#OTSLine3D">OTSLine3D line = new OTSLine3D(new OTSPoint3Dometry/OTSPoint3D.html#OTSPoint3D">OTSPoint3D[] {new OTSPoint3D(-579.253, 60.157, 4.710),
556                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-579.253, 60.144, 4.712), new OTSPoint3D(-579.253, 60.144, 0.000),
557                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-579.251, 60.044, 0.000), new OTSPoint3D(-579.246, 59.944, 0.000),
558                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-579.236, 59.845, 0.000), new OTSPoint3D(-579.223, 59.746, 0.000),
559                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-579.206, 59.647, 0.000), new OTSPoint3D(-579.185, 59.549, 0.000),
560                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-579.161, 59.452, 0.000), new OTSPoint3D(-579.133, 59.356, 0.000),
561                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-579.101, 59.261, 0.000), new OTSPoint3D(-579.066, 59.168, 0.000),
562                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-579.028, 59.075, 0.000), new OTSPoint3D(-578.986, 58.985, 0.000),
563                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-578.940, 58.896, 0.000), new OTSPoint3D(-578.891, 58.809, 0.000),
564                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-578.839, 58.723, 0.000), new OTSPoint3D(-578.784, 58.640, 0.000),
565                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-578.725, 58.559, 0.000), new OTSPoint3D(-578.664, 58.480, 0.000),
566                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-578.599, 58.403, 0.000), new OTSPoint3D(-578.532, 58.329, 0.000),
567                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-578.462, 58.258, 0.000), new OTSPoint3D(-578.390, 58.189, 0.000),
568                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-578.314, 58.123, 0.000), new OTSPoint3D(-578.237, 58.060, 0.000),
569                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-578.157, 58.000, 0.000), new OTSPoint3D(-578.075, 57.943, 0.000),
570                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-577.990, 57.889, 0.000), new OTSPoint3D(-577.904, 57.839, 0.000),
571                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-577.816, 57.791, 0.000), new OTSPoint3D(-577.726, 57.747, 0.000),
572                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-577.635, 57.707, 0.000), new OTSPoint3D(-577.542, 57.670, 0.000),
573                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-577.448, 57.636, 0.000), new OTSPoint3D(-577.352, 57.606, 0.000),
574                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-577.256, 57.580, 0.000), new OTSPoint3D(-577.159, 57.557, 0.000),
575                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-577.060, 57.538, 0.000), new OTSPoint3D(-576.962, 57.523, 0.000),
576                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-576.862, 57.512, 0.000), new OTSPoint3D(-576.763, 57.504, 0.000),
577                 new OTSPoint3DOTSPoint3D">OTSPoint3D(-576.663, 57.500, 0.000), new OTSPoint3D(-576.623, 57.500, 6.278),
578                 new OTSPoint3DOTSPoint3D">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 }