LCOV - code coverage report
Current view: top level - src/netbuild - NBNodeShapeComputer.cpp (source / functions) Hit Total Coverage
Test: lcov.info Lines: 412 448 92.0 %
Date: 2024-05-07 15:28:01 Functions: 18 18 100.0 %

          Line data    Source code
       1             : /****************************************************************************/
       2             : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
       3             : // Copyright (C) 2001-2024 German Aerospace Center (DLR) and others.
       4             : // This program and the accompanying materials are made available under the
       5             : // terms of the Eclipse Public License 2.0 which is available at
       6             : // https://www.eclipse.org/legal/epl-2.0/
       7             : // This Source Code may also be made available under the following Secondary
       8             : // Licenses when the conditions for such availability set forth in the Eclipse
       9             : // Public License 2.0 are satisfied: GNU General Public License, version 2
      10             : // or later which is available at
      11             : // https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
      12             : // SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
      13             : /****************************************************************************/
      14             : /// @file    NBNodeShapeComputer.cpp
      15             : /// @author  Daniel Krajzewicz
      16             : /// @author  Jakob Erdmann
      17             : /// @author  Michael Behrisch
      18             : /// @date    Sept 2002
      19             : ///
      20             : // This class computes shapes of junctions
      21             : /****************************************************************************/
      22             : #include <config.h>
      23             : 
      24             : #include <algorithm>
      25             : #include <iterator>
      26             : #include <utils/geom/PositionVector.h>
      27             : #include <utils/options/OptionsCont.h>
      28             : #include <utils/geom/GeomHelper.h>
      29             : #include <utils/common/StdDefs.h>
      30             : #include <utils/common/MsgHandler.h>
      31             : #include <utils/common/UtilExceptions.h>
      32             : #include <utils/common/ToString.h>
      33             : #include <utils/iodevices/OutputDevice.h>
      34             : #include "NBNode.h"
      35             : #include "NBAlgorithms.h"
      36             : #include "NBNodeShapeComputer.h"
      37             : 
      38             : //#define DEBUG_NODE_SHAPE
      39             : //#define DEBUG_SMOOTH_CORNERS
      40             : //#define DEBUG_RADIUS
      41             : #define DEBUGCOND (myNode.getID() == "C")
      42             : 
      43             : 
      44             : #define EXT2 10.0
      45             : 
      46             : // foot and bicycle paths as well as pure service roads should not get large junctions
      47             : // railways also do have have junctions with sharp turns so can be excluded
      48             : const SVCPermissions NBNodeShapeComputer::SVC_LARGE_TURN(
      49             :     SVCAll & ~(SVC_BICYCLE | SVC_PEDESTRIAN | SVC_DELIVERY | SVC_RAIL_CLASSES));
      50             : 
      51             : // ===========================================================================
      52             : // method definitions
      53             : // ===========================================================================
      54      106598 : NBNodeShapeComputer::NBNodeShapeComputer(const NBNode& node) :
      55      106598 :     myNode(node),
      56      106598 :     myRadius(node.getRadius()) {
      57      106598 :     if (node.getEdges().size() > 4 && !NBNodeTypeComputer::isRailwayNode(&node)) {
      58       25706 :         EXT = 50;
      59             :     } else {
      60       80892 :         EXT = 100;
      61             :     }
      62      106598 : }
      63             : 
      64             : 
      65      106598 : NBNodeShapeComputer::~NBNodeShapeComputer() {}
      66             : 
      67             : 
      68             : const PositionVector
      69      106598 : NBNodeShapeComputer::compute(bool forceSmall) {
      70             : #ifdef DEBUG_NODE_SHAPE
      71             :     if (DEBUGCOND) {
      72             :         // annotate edges edges to make their ordering visible
      73             :         int i = 0;
      74             :         for (NBEdge* e : myNode.getEdges()) {
      75             :             e->setStreetName(toString(i));
      76             :             i++;
      77             :         }
      78             :     }
      79             : #endif
      80             :     // check whether the node is a dead end node or a node where only turning is possible
      81             :     //  in this case, we will use "computeNodeShapeSmall"
      82      106598 :     if (myNode.getEdges().size() == 1 || forceSmall) {
      83       13295 :         return computeNodeShapeSmall();
      84             :     }
      85       93303 :     if (myNode.getEdges().size() == 2 && myNode.getIncomingEdges().size() == 1) {
      86       29184 :         if (myNode.getIncomingEdges()[0]->isTurningDirectionAt(myNode.getOutgoingEdges()[0])) {
      87       13826 :             return computeNodeShapeSmall();
      88             :         }
      89             :     }
      90       79477 :     const bool geometryLike = myNode.isSimpleContinuation(true, true);
      91       79477 :     const PositionVector& ret = computeNodeShapeDefault(geometryLike);
      92             :     // fail fall-back: use "computeNodeShapeSmall"
      93       79477 :     if (ret.size() < 3) {
      94         105 :         return computeNodeShapeSmall();
      95             :     }
      96             :     return ret;
      97       79477 : }
      98             : 
      99             : 
     100             : void
     101       97189 : NBNodeShapeComputer::computeSameEnd(PositionVector& l1, PositionVector& l2) {
     102             :     assert(l1[0].distanceTo2D(l1[1]) >= EXT);
     103             :     assert(l2[0].distanceTo2D(l2[1]) >= EXT);
     104       97189 :     PositionVector tmp;
     105       97189 :     tmp.push_back(PositionVector::positionAtOffset2D(l1[0], l1[1], EXT));
     106       97189 :     tmp.push_back(l1[1]);
     107       97189 :     tmp[1].sub(tmp[0]);
     108       97189 :     tmp[1].set(-tmp[1].y(), tmp[1].x());
     109       97189 :     tmp[1].add(tmp[0]);
     110       97189 :     tmp.extrapolate2D(EXT);
     111       97189 :     if (l2.intersects(tmp[0], tmp[1])) {
     112       97102 :         const double offset = l2.intersectsAtLengths2D(tmp)[0];
     113       97102 :         if (l2.length2D() - offset > POSITION_EPS) {
     114       97102 :             PositionVector tl2 = l2.getSubpart2D(offset, l2.length2D());
     115       97102 :             tl2.extrapolate2D(EXT);
     116       97102 :             l2.erase(l2.begin(), l2.begin() + (l2.size() - tl2.size()));
     117       97102 :             l2[0] = tl2[0];
     118       97102 :         }
     119             :     }
     120       97189 : }
     121             : 
     122             : 
     123             : const PositionVector
     124       79477 : NBNodeShapeComputer::computeNodeShapeDefault(bool simpleContinuation) {
     125             :     // if we have less than two edges, we can not compute the node's shape this way
     126       79477 :     if (myNode.getEdges().size() < 2) {
     127           0 :         return PositionVector();
     128             :     }
     129             :     // magic values
     130       79477 :     const OptionsCont& oc = OptionsCont::getOptions();
     131       79477 :     const double defaultRadius = getDefaultRadius(oc);
     132       79477 :     const bool useDefaultRadius = myNode.getRadius() == NBNode::UNSPECIFIED_RADIUS || myNode.getRadius() == defaultRadius;
     133       79477 :     myRadius = (useDefaultRadius ? defaultRadius : myNode.getRadius());
     134       79477 :     double smallRadius = useDefaultRadius ? oc.getFloat("junctions.small-radius") : myRadius;
     135       79477 :     const int cornerDetail = oc.getInt("junctions.corner-detail");
     136       79477 :     const double sCurveStretch = oc.getFloat("junctions.scurve-stretch");
     137       79477 :     const bool useEndpoints = oc.getBool("junctions.endpoint-shape");
     138       79477 :     const bool rectangularCut = oc.getBool("rectangular-lane-cut");
     139       79477 :     const bool openDriveOutput = oc.isSet("opendrive-output");
     140             : 
     141             :     // Extend geometries to move the stop line forward.
     142             :     // In OpenDrive the junction starts whenever the geometry changes. Stop
     143             :     // line information is not given or ambiguous (sign positions at most)
     144             :     // In SUMO, stop lines are where the junction starts. This is computed
     145             :     // heuristically from intersecting the junctions roads geometries.
     146      154163 :     const double advanceStopLine = oc.exists("opendrive-files") && oc.isSet("opendrive-files") ? oc.getFloat("opendrive.advance-stopline") : 0;
     147             : 
     148             : 
     149             : #ifdef DEBUG_NODE_SHAPE
     150             :     if (DEBUGCOND) {
     151             :         std::cout << "\ncomputeNodeShapeDefault node " << myNode.getID() << " simple=" << simpleContinuation << " useDefaultRadius=" << useDefaultRadius << " radius=" << myRadius << "\n";
     152             :     }
     153             : #endif
     154             : 
     155             :     // initialise
     156       79477 :     EdgeVector::const_iterator i;
     157             :     // edges located in the value-vector have the same direction as the key edge
     158             :     std::map<NBEdge*, std::set<NBEdge*> > same;
     159             :     // the counter-clockwise boundary of the edge regarding possible same-direction edges
     160             :     GeomsMap geomsCCW;
     161             :     // the clockwise boundary of the edge regarding possible same-direction edges
     162             :     GeomsMap geomsCW;
     163       79477 :     EdgeVector usedEdges = myNode.getEdges();
     164       79477 :     computeEdgeBoundaries(usedEdges, geomsCCW, geomsCW);
     165             : 
     166             :     // check which edges are parallel
     167       79477 :     joinSameDirectionEdges(usedEdges, same, useEndpoints);
     168             :     // compute unique direction list
     169       79477 :     EdgeVector newAll = computeUniqueDirectionList(usedEdges, same, geomsCCW, geomsCW);
     170             :     // if we have only two "directions", let's not compute the geometry using this method
     171       79477 :     if (newAll.size() < 2) {
     172         105 :         return PositionVector();
     173             :     }
     174             : 
     175             :     // All geoms are outgoing from myNode.
     176             :     // for every direction in newAll we compute the offset at which the
     177             :     // intersection ends and the edge starts. This value is saved in 'distances'
     178             :     // If the geometries need to be extended to get an intersection, this is
     179             :     // recorded in 'myExtended'
     180             :     std::map<NBEdge*, double> distances;
     181             :     std::map<NBEdge*, double> distances2;
     182             :     std::map<NBEdge*, bool> myExtended;
     183             : 
     184      307344 :     for (i = newAll.begin(); i != newAll.end(); ++i) {
     185      227972 :         EdgeVector::const_iterator cwi = i;
     186      227972 :         EdgeVector::const_iterator ccwi = i;
     187             :         double ccad;
     188             :         double cad;
     189      227972 :         initNeighbors(newAll, i, geomsCW, geomsCCW, cwi, ccwi, cad, ccad);
     190             :         assert(geomsCCW.find(*i) != geomsCCW.end());
     191             :         assert(geomsCW.find(*ccwi) != geomsCW.end());
     192             :         assert(geomsCW.find(*cwi) != geomsCW.end());
     193             : 
     194             :         // there are only 2 directions and they are almost parallel
     195      227972 :         if (*cwi == *ccwi &&
     196             :                 (
     197             :                     // no change in lane numbers, even low angles still give a good intersection
     198       39864 :                     (simpleContinuation && fabs(ccad - cad) < (double) 0.1)
     199             :                     // lane numbers change, a direct intersection could be far away from the node position
     200             :                     // so we use a larger threshold
     201       32402 :                     || (!simpleContinuation && fabs(ccad - cad) < DEG2RAD(22.5)))
     202             :            ) {
     203             :             // compute the mean position between both edges ends ...
     204             :             Position p;
     205       38110 :             if (myExtended.find(*ccwi) != myExtended.end()) {
     206           0 :                 p = geomsCCW[*ccwi][0];
     207           0 :                 p.add(geomsCW[*ccwi][0]);
     208             :                 p.mul(0.5);
     209             : #ifdef DEBUG_NODE_SHAPE
     210             :                 if (DEBUGCOND) {
     211             :                     std::cout << " extended: p=" << p << " angle=" << (ccad - cad) << "\n";
     212             :                 }
     213             : #endif
     214             :             } else {
     215       38110 :                 p = geomsCCW[*ccwi][0];
     216       38110 :                 p.add(geomsCW[*ccwi][0]);
     217       38110 :                 p.add(geomsCCW[*i][0]);
     218       38110 :                 p.add(geomsCW[*i][0]);
     219             :                 p.mul(0.25);
     220             : #ifdef DEBUG_NODE_SHAPE
     221             :                 if (DEBUGCOND) {
     222             :                     std::cout << " unextended: p=" << p << " angle=" << (ccad - cad) << "\n";
     223             :                 }
     224             : #endif
     225             :             }
     226             :             // ... compute the distance to this point ...
     227       76220 :             double dist = MAX2(
     228       38110 :                               geomsCCW[*i].nearest_offset_to_point2D(p),
     229       38110 :                               geomsCW[*i].nearest_offset_to_point2D(p));
     230       38110 :             if (dist < 0) {
     231           0 :                 if (isRailway((*i)->getPermissions())) {
     232             :                     // better not mess up bidi geometries
     233           0 :                     return PositionVector();
     234             :                 }
     235             :                 // ok, we have the problem that even the extrapolated geometry
     236             :                 //  does not reach the point
     237             :                 // in this case, the geometry has to be extenden... too bad ...
     238             :                 // ... let's append the mean position to the geometry
     239           0 :                 PositionVector g = (*i)->getGeometry();
     240           0 :                 if (myNode.hasIncoming(*i)) {
     241           0 :                     g.push_back_noDoublePos(p);
     242             :                 } else {
     243           0 :                     g.push_front_noDoublePos(p);
     244             :                 }
     245           0 :                 (*i)->setGeometry(g);
     246             :                 // and rebuild previous information
     247           0 :                 geomsCCW[*i] = (*i)->getCCWBoundaryLine(myNode);
     248           0 :                 geomsCCW[*i].extrapolate(EXT);
     249           0 :                 geomsCW[*i] = (*i)->getCWBoundaryLine(myNode);
     250           0 :                 geomsCW[*i].extrapolate(EXT);
     251             :                 // the distance is now = zero (the point we have appended)
     252           0 :                 distances[*i] = EXT;
     253           0 :                 myExtended[*i] = true;
     254             : #ifdef DEBUG_NODE_SHAPE
     255             :                 if (DEBUGCOND) {
     256             :                     std::cout << " extending (dist=" << dist << ")\n";
     257             :                 }
     258             : #endif
     259           0 :             } else {
     260       38110 :                 if (!simpleContinuation) {
     261       13246 :                     dist += myRadius;
     262             :                 } else {
     263             :                     // if the angles change, junction should have some size to avoid degenerate shape
     264       24864 :                     double radius2 = fabs(ccad - cad) * (*i)->getNumLanes();
     265       24864 :                     if (radius2 > NUMERICAL_EPS || openDriveOutput) {
     266             :                         radius2 = MAX2(0.15, radius2);
     267             :                     }
     268       49728 :                     if (myNode.getCrossings().size() > 0) {
     269          42 :                         double width = myNode.getCrossings()[0]->customWidth;
     270          42 :                         if (width == NBEdge::UNSPECIFIED_WIDTH) {
     271          76 :                             width = OptionsCont::getOptions().getFloat("default.crossing-width");
     272             :                         }
     273          42 :                         radius2 = MAX2(radius2, width / 2);
     274             :                     }
     275       24864 :                     if (!useDefaultRadius) {
     276           2 :                         radius2 = MAX2(radius2, myRadius);
     277             :                     }
     278       24864 :                     dist += radius2;
     279             : #ifdef DEBUG_NODE_SHAPE
     280             :                     if (DEBUGCOND) {
     281             :                         std::cout << " using radius=" << radius2 << " ccad=" << ccad << " cad=" << cad << "\n";
     282             :                     }
     283             : #endif
     284             :                 }
     285       38110 :                 distances[*i] = dist;
     286             :             }
     287             : 
     288             :         } else {
     289             :             // the angles are different enough to compute the intersection of
     290             :             // the outer boundaries directly (or there are more than 2 directions). The "nearer" neighbor causes the furthest distance
     291      189862 :             const bool ccwCloser = ccad < cad;
     292      189862 :             const bool cwLargeTurn = needsLargeTurn(*i, *cwi, same);
     293      189862 :             const bool ccwLargeTurn = needsLargeTurn(*i, *ccwi, same);
     294      189862 :             const bool neighLargeTurn = ccwCloser ? ccwLargeTurn : cwLargeTurn;
     295      189862 :             const bool neigh2LargeTurn =  ccwCloser ? cwLargeTurn : ccwLargeTurn;
     296             :             // the border facing the closer neighbor
     297      189862 :             const PositionVector& currGeom = ccwCloser ? geomsCCW[*i] : geomsCW[*i];
     298             :             // the border facing the far neighbor
     299      189862 :             const PositionVector& currGeom2 = ccwCloser ? geomsCW[*i] : geomsCCW[*i];
     300             :             // the border of the closer neighbor
     301      189862 :             const PositionVector& neighGeom = ccwCloser ? geomsCW[*ccwi] : geomsCCW[*cwi];
     302             :             // the border of the far neighbor
     303      189862 :             const PositionVector& neighGeom2 = ccwCloser ? geomsCCW[*cwi] : geomsCW[*ccwi];
     304             :             // whether the current edge/direction spans a divided road
     305      379724 :             const bool keepBothDistances = isDivided(*i, same[*i], geomsCCW[*i], geomsCW[*i]);
     306             : #ifdef DEBUG_NODE_SHAPE
     307             :             if (DEBUGCOND) {
     308             :                 std::cout << " i=" << (*i)->getID() << " neigh=" << (*ccwi)->getID() << " neigh2=" << (*cwi)->getID() << "\n";
     309             :                 std::cout << "    ccwCloser=" << ccwCloser << " divided=" << keepBothDistances
     310             :                           << "\n      currGeom=" << currGeom << " neighGeom=" << neighGeom
     311             :                           << "\n      currGeom2=" << currGeom2 << " neighGeom2=" << neighGeom2
     312             :                           << "\n";
     313             :             }
     314             : #endif
     315      189862 :             if (!simpleContinuation) {
     316      174763 :                 if (useEndpoints && !(*i)->hasDefaultGeometryEndpointAtNode(&myNode)) {
     317          16 :                     distances[*i] = EXT;
     318      174747 :                 } else if (currGeom.intersects(neighGeom)) {
     319      173933 :                     distances[*i] = (neighLargeTurn ? myRadius : smallRadius) + closestIntersection(currGeom, neighGeom, EXT);
     320             : #ifdef DEBUG_NODE_SHAPE
     321             :                     if (DEBUGCOND) {
     322             :                         std::cout << "   neigh intersects dist=" << distances[*i] << " currGeom=" << currGeom << " neighGeom=" << neighGeom << "\n";
     323             :                     }
     324             : #endif
     325      173933 :                     if (*cwi != *ccwi && currGeom2.intersects(neighGeom2)) {
     326             :                         // also use the second intersection point
     327             :                         // but prevent very large node shapes
     328      166240 :                         const double farAngleDist = ccwCloser ? cad : ccad;
     329      166240 :                         double a1 = distances[*i];
     330      166240 :                         double a2 = (neigh2LargeTurn ? myRadius : smallRadius) + closestIntersection(currGeom2, neighGeom2, EXT);
     331             : #ifdef DEBUG_NODE_SHAPE
     332             :                         if (DEBUGCOND) {
     333             :                             std::cout << "      neigh2 also intersects a1=" << a1 << " a2=" << a2 << " ccad=" << RAD2DEG(ccad) << " cad=" << RAD2DEG(cad) << " dist[cwi]=" << distances[*cwi] << " dist[ccwi]=" << distances[*ccwi] << " farAngleDist=" << RAD2DEG(farAngleDist) << " currGeom2=" << currGeom2 << " neighGeom2=" << neighGeom2 << "\n";
     334             :                         }
     335             : #endif
     336             :                         //if (RAD2DEG(farAngleDist) < 175) {
     337             :                         //    distances[*i] = MAX2(a1, MIN2(a2, a1 + 180 - RAD2DEG(farAngleDist)));
     338             :                         //}
     339      166240 :                         if (a2 <= EXT) {
     340        5851 :                             if (keepBothDistances) {
     341          73 :                                 if (ccwCloser) {
     342          19 :                                     distances2[*i] = a2;
     343             :                                 } else {
     344          54 :                                     distances[*i] = a2;
     345          54 :                                     distances2[*i] = a1;
     346             :                                 }
     347             :                             } else {
     348       11556 :                                 distances[*i] = MAX2(a1, a2);
     349             :                             }
     350      160389 :                         } else if (ccad > DEG2RAD(90. + 45.) && cad > DEG2RAD(90. + 45.)) {
     351             :                             // do nothing.
     352      155115 :                         } else if (farAngleDist < DEG2RAD(135) || (fabs(RAD2DEG(farAngleDist) - 180) > 1 && fabs(a2 - a1) < 10)) {
     353      120898 :                             if (keepBothDistances) {
     354         367 :                                 if (ccwCloser) {
     355         171 :                                     distances2[*i] = a2;
     356             :                                 } else {
     357         196 :                                     distances[*i] = a2;
     358         196 :                                     distances2[*i] = a1;
     359             :                                 }
     360             :                             } else {
     361      241062 :                                 distances[*i] = MAX2(a1, a2);
     362             :                             }
     363             :                         }
     364             : #ifdef DEBUG_NODE_SHAPE
     365             :                         if (DEBUGCOND) {
     366             :                             std::cout << "   a1=" << a1 << " a2=" << a2 << " keepBoth=" << keepBothDistances << " dist=" << distances[*i] << "\n";
     367             :                         }
     368             : #endif
     369             :                     }
     370             :                 } else {
     371         814 :                     if (*cwi != *ccwi && currGeom2.intersects(neighGeom2)) {
     372        1426 :                         distances[*i] = (neigh2LargeTurn ? myRadius : smallRadius) + currGeom2.intersectsAtLengths2D(neighGeom2)[0];
     373             : #ifdef DEBUG_NODE_SHAPE
     374             :                         if (DEBUGCOND) {
     375             :                             std::cout << "   neigh2 intersects dist=" << distances[*i] << " currGeom2=" << currGeom2 << " neighGeom2=" << neighGeom2 << "\n";
     376             :                         }
     377             : #endif
     378             :                     } else {
     379         101 :                         distances[*i] = EXT + myRadius;
     380             : #ifdef DEBUG_NODE_SHAPE
     381             :                         if (DEBUGCOND) {
     382             :                             std::cout << "   no intersects dist=" << distances[*i]  << " currGeom=" << currGeom << " neighGeom=" << neighGeom << " currGeom2=" << currGeom2 << " neighGeom2=" << neighGeom2 << "\n";
     383             :                         }
     384             : #endif
     385             :                     }
     386             :                 }
     387             :             } else {
     388       15099 :                 if (currGeom.intersects(neighGeom)) {
     389       30128 :                     distances[*i] = currGeom.intersectsAtLengths2D(neighGeom)[0];
     390             :                 } else {
     391          35 :                     distances[*i] = (double) EXT;
     392             :                 }
     393             :             }
     394             :         }
     395      227972 :         if (useDefaultRadius && sCurveStretch > 0) {
     396          32 :             double sCurveWidth = myNode.getDisplacementError();
     397          32 :             if (sCurveWidth > 0) {
     398           2 :                 const double sCurveRadius = myRadius + sCurveWidth / SUMO_const_laneWidth * sCurveStretch * pow((*i)->getSpeed(), 2 + sCurveStretch) / 1000;
     399           2 :                 const double stretch = EXT + sCurveRadius - distances[*i];
     400           2 :                 if (stretch > 0) {
     401           2 :                     distances[*i] += stretch;
     402             :                     // fixate extended geometry for repeated computation
     403           2 :                     const double shorten = distances[*i] - EXT;
     404           2 :                     (*i)->shortenGeometryAtNode(&myNode, shorten);
     405           2 :                     for (std::set<NBEdge*>::iterator k = same[*i].begin(); k != same[*i].end(); ++k) {
     406           0 :                         (*k)->shortenGeometryAtNode(&myNode, shorten);
     407             :                     }
     408             : #ifdef DEBUG_NODE_SHAPE
     409             :                     if (DEBUGCOND) {
     410             :                         std::cout << "   stretching junction: sCurveWidth=" << sCurveWidth << " sCurveRadius=" << sCurveRadius << " stretch=" << stretch << " dist=" << distances[*i]  << "\n";
     411             :                     }
     412             : #endif
     413             :                 }
     414             :             }
     415             :         }
     416             :     }
     417             : 
     418      307344 :     for (NBEdge* const edge : newAll) {
     419      227972 :         if (distances.find(edge) == distances.end()) {
     420             :             assert(false);
     421           0 :             distances[edge] = EXT;
     422             :         }
     423             :     }
     424             :     // because of lane spread right the crossing point may be identical to the junction center and thus the distance is exactly EXT
     425       79372 :     const double off = EXT - NUMERICAL_EPS;
     426             :     // prevent inverted node shapes
     427             :     // (may happen with near-parallel edges)
     428       79372 :     const double minDistSum = 2 * (EXT + myRadius);
     429      307344 :     for (NBEdge* const edge : newAll) {
     430      227972 :         if (distances[edge] < off && edge->hasDefaultGeometryEndpointAtNode(&myNode)) {
     431        8838 :             for (EdgeVector::const_iterator j = newAll.begin(); j != newAll.end(); ++j) {
     432        6602 :                 if (distances[*j] > off && (*j)->hasDefaultGeometryEndpointAtNode(&myNode) && distances[edge] + distances[*j] < minDistSum) {
     433         370 :                     const double angleDiff = fabs(NBHelpers::relAngle(edge->getAngleAtNode(&myNode), (*j)->getAngleAtNode(&myNode)));
     434         370 :                     if (angleDiff > 160 || angleDiff < 20) {
     435             : #ifdef DEBUG_NODE_SHAPE
     436             :                         if (DEBUGCOND) {
     437             :                             std::cout << "   increasing dist for i=" << edge->getID() << " because of j=" << (*j)->getID() << " jDist=" << distances[*j]
     438             :                                       << "  oldI=" << distances[edge] << " newI=" << minDistSum - distances[*j]
     439             :                                       << " angleDiff=" << angleDiff
     440             :                                       << " geomI=" << edge->getGeometry() << " geomJ=" << (*j)->getGeometry() << "\n";
     441             :                         }
     442             : #endif
     443         190 :                         distances[edge] = minDistSum - distances[*j];
     444             :                     }
     445             :                 }
     446             :             }
     447             :         }
     448             :     }
     449             : 
     450             : 
     451             :     // build
     452       79372 :     PositionVector ret;
     453      307344 :     for (i = newAll.begin(); i != newAll.end(); ++i) {
     454      227972 :         const PositionVector& ccwBound = geomsCCW[*i];
     455      227972 :         const PositionVector& cwBound = geomsCW[*i];
     456             :         //double offset = MIN3(distances[*i], cwBound.length2D() - POSITION_EPS, ccwBound.length2D() - POSITION_EPS);
     457      227972 :         double offset = distances[*i];
     458         440 :         double offset2 = distances2.count(*i) != 0 ? distances2[*i] : offset;
     459      227972 :         if (offset != offset2) {
     460             :             // keep rectangular cuts if the difference is small or the roads aren't
     461             :             // really divided by much (unless the angle is very different)
     462        1320 :             const double dWidth = divisionWidth(*i, same[*i],
     463         440 :                                                 ccwBound.positionAtOffset2D(offset),
     464         440 :                                                 cwBound.positionAtOffset2D(offset2));
     465         440 :             const double angle = RAD2DEG(GeomHelper::angleDiff(ccwBound.angleAt2D(0), cwBound.angleAt2D(0)));
     466         440 :             const double oDelta = fabs(offset - offset2);
     467             :             //std::cout << " i=" << (*i)->getID() << " offset=" << offset << " offset2=" << offset2 << " dWidth=" << dWidth << " angle=" << angle << " same=" << joinNamedToStringSorting(same[*i], ",") << "\n";
     468         440 :             if ((((oDelta < 5 || dWidth < 10) && fabs(angle) < 30)) || (fabs(angle) < 5 && myNode.getType() != SumoXMLNodeType::RAIL_CROSSING)) {
     469             : #ifdef DEBUG_NODE_SHAPE
     470             :                 std::cout << " i=" << (*i)->getID() << " offset=" << offset << " offset2=" << offset2 << " dWidth=" << dWidth << " angle=" << angle << " same=" << joinNamedToStringSorting(same[*i], ",") << "\n";
     471             : #endif
     472             :                 offset = MAX2(offset, offset2);
     473             :                 offset2 = offset;
     474             :             }
     475             :         }
     476      227972 :         if (!(*i)->hasDefaultGeometryEndpointAtNode(&myNode)) {
     477             :             // for non geometry-endpoints, only shorten but never extend the geometry
     478       23702 :             if (advanceStopLine > 0 && offset < EXT) {
     479             : #ifdef DEBUG_NODE_SHAPE
     480             :                 std::cout << " i=" << (*i)->getID() << " offset=" << offset << " advanceStopLine=" << advanceStopLine << "\n";
     481             : #endif
     482             :                 // fixate extended geometry for repeated computation
     483           0 :                 (*i)->extendGeometryAtNode(&myNode, advanceStopLine);
     484           0 :                 for (std::set<NBEdge*>::iterator k = same[*i].begin(); k != same[*i].end(); ++k) {
     485           0 :                     (*k)->extendGeometryAtNode(&myNode, advanceStopLine);
     486             :                 }
     487             :             }
     488       23702 :             offset = MAX2(EXT - advanceStopLine, offset);
     489             :             offset2 = MAX2(EXT - advanceStopLine, offset2);
     490             :         }
     491      227972 :         if (offset == -1) {
     492           0 :             WRITE_WARNINGF(TL("Fixing offset for edge '%' at node '%."), (*i)->getID(), myNode.getID());
     493             :             offset = -.1;
     494             :             offset2 = -.1;
     495             :         }
     496      227972 :         Position p = ccwBound.positionAtOffset2D(offset);
     497      227972 :         p.setz(myNode.getPosition().z());
     498      227972 :         if (i != newAll.begin()) {
     499      445800 :             ret.append(getSmoothCorner(geomsCW[*(i - 1)], ccwBound, ret[-1], p, cornerDetail));
     500             :         }
     501      227972 :         Position p2 = cwBound.positionAtOffset2D(offset2);
     502      227972 :         p2.setz(myNode.getPosition().z());
     503             :         //ret.append(getEdgeCuts(*i, geomsCCW, geomsCW, offset, offset2, same));
     504      227972 :         ret.push_back_noDoublePos(p);
     505      227972 :         ret.push_back_noDoublePos(p2);
     506             : #ifdef DEBUG_NODE_SHAPE
     507             :         if (DEBUGCOND) {
     508             :             std::cout << "   build stopLine for i=" << (*i)->getID() << " offset=" << offset << " offset2=" << offset2 << " dist=" << distances[*i] << " cwLength=" << cwBound.length2D() << " ccwLength=" << ccwBound.length2D() << " p=" << p << " p2=" << p2 << " ccwBound=" <<  ccwBound << " cwBound=" << cwBound << "\n";
     509             :         }
     510             : #endif
     511      227972 :         (*i)->setNodeBorder(&myNode, p, p2, rectangularCut);
     512      327954 :         for (std::set<NBEdge*>::iterator k = same[*i].begin(); k != same[*i].end(); ++k) {
     513       99982 :             (*k)->setNodeBorder(&myNode, p, p2, rectangularCut);
     514             :         }
     515             :     }
     516             :     // final curve segment
     517      238116 :     ret.append(getSmoothCorner(geomsCW[*(newAll.end() - 1)], geomsCCW[*newAll.begin()], ret[-1], ret[0], cornerDetail));
     518             : #ifdef DEBUG_NODE_SHAPE
     519             :     if (DEBUGCOND) {
     520             :         std::cout << " final shape=" << ret << "\n";
     521             :     }
     522             : #endif
     523             :     return ret;
     524       79372 : }
     525             : 
     526             : 
     527             : double
     528      340173 : NBNodeShapeComputer::closestIntersection(const PositionVector& geom1, const PositionVector& geom2, double offset) {
     529      340173 :     std::vector<double> intersections = geom1.intersectsAtLengths2D(geom2);
     530      340173 :     double result = intersections[0];
     531      360977 :     for (std::vector<double>::iterator it = intersections.begin() + 1; it != intersections.end(); ++it) {
     532       20804 :         if (fabs(*it - offset) < fabs(result - offset)) {
     533             :             result = *it;
     534             :         }
     535             :     }
     536      340173 :     return result;
     537             : }
     538             : 
     539             : bool
     540      379724 : NBNodeShapeComputer::needsLargeTurn(NBEdge* e1, NBEdge* e2,
     541             :                                     std::map<NBEdge*, std::set<NBEdge*> >& same) const {
     542      379724 :     const SVCPermissions p1 = e1->getPermissions();
     543      379724 :     const SVCPermissions p2 = e2->getPermissions();
     544      379724 :     if ((p1 & p2 & SVC_LARGE_TURN) != 0) {
     545             :         // note: would could also check whether there is actually a connection
     546             :         // between those edges
     547             :         return true;
     548             :     }
     549             :     // maybe edges in the same direction need a large turn
     550      254839 :     for (NBEdge* e2s : same[e2]) {
     551       59515 :         if ((p1 & e2s->getPermissions() & SVC_LARGE_TURN) != 0
     552       59515 :                 && (e1->getToNode() == e2s->getFromNode() || e2s->getToNode() == e1->getFromNode())) {
     553             :             return true;
     554             :         }
     555       99035 :         for (NBEdge* e1s : same[e1]) {
     556       41048 :             if ((e2s->getPermissions() & e1s->getPermissions() & SVC_LARGE_TURN) != 0
     557       41048 :                     && (e2s->getToNode() == e1s->getFromNode() || e1s->getToNode() == e2s->getFromNode())) {
     558             :                 return true;
     559             :             }
     560             :         }
     561             :     }
     562      252476 :     for (NBEdge* e1s : same[e1]) {
     563       57922 :         if ((p2 & e1s->getPermissions() & SVC_LARGE_TURN) != 0
     564       57922 :                 && (e2->getToNode() == e1s->getFromNode() || e1s->getToNode() == e2->getFromNode())) {
     565             :             return true;
     566             :         }
     567             :     }
     568             :     //std::cout << " e1=" << e1->getID() << " e2=" << e2->getID() << " sameE1=" << toString(same[e1]) << " sameE2=" << toString(same[e2]) << "\n";
     569             :     return false;
     570             : }
     571             : 
     572             : PositionVector
     573      227972 : NBNodeShapeComputer::getSmoothCorner(PositionVector begShape, PositionVector endShape,
     574             :                                      const Position& begPoint, const Position& endPoint, int cornerDetail) {
     575      227972 :     PositionVector ret;
     576      227972 :     if (cornerDetail > 0) {
     577      185409 :         PositionVector begShape2 = begShape.reverse().getSubpart2D(EXT2, begShape.length());
     578      185409 :         const double begSplit = begShape2.nearest_offset_to_point2D(begPoint, false);
     579             : #ifdef DEBUG_SMOOTH_CORNERS
     580             :         if (DEBUGCOND) {
     581             :             std::cout << " begLength=" << begShape2.length2D() << " begSplit=" << begSplit << "\n";
     582             :         }
     583             : #endif
     584      185409 :         if (begSplit > POSITION_EPS && begSplit < begShape2.length2D() - POSITION_EPS) {
     585      168738 :             begShape2 = begShape2.splitAt(begSplit, true).first;
     586             :         } else {
     587             :             return ret;
     588             :         }
     589      168738 :         PositionVector endShape2 = endShape.getSubpart(0, endShape.length() - EXT2);
     590      168738 :         const double endSplit = endShape2.nearest_offset_to_point2D(endPoint, false);
     591             : #ifdef DEBUG_SMOOTH_CORNERS
     592             :         if (DEBUGCOND) {
     593             :             std::cout << " endLength=" << endShape2.length2D() << " endSplit=" << endSplit << "\n";
     594             :         }
     595             : #endif
     596      168738 :         if (endSplit > POSITION_EPS && endSplit < endShape2.length2D() - POSITION_EPS) {
     597      160192 :             endShape2 = endShape2.splitAt(endSplit, true).second;
     598             :         } else {
     599             :             return ret;
     600             :         }
     601             :         // flatten z to junction z level
     602      320384 :         begShape2 = begShape2.interpolateZ(myNode.getPosition().z(), myNode.getPosition().z());
     603      320384 :         endShape2 = endShape2.interpolateZ(myNode.getPosition().z(), myNode.getPosition().z());
     604             : #ifdef DEBUG_SMOOTH_CORNERS
     605             :         if (DEBUGCOND) {
     606             :             std::cout << "getSmoothCorner begPoint=" << begPoint << " endPoint=" << endPoint
     607             :                       << " begShape=" << begShape << " endShape=" << endShape
     608             :                       << " begShape2=" << begShape2 << " endShape2=" << endShape2
     609             :                       << "\n";
     610             :         }
     611             : #endif
     612      160192 :         if (begShape2.size() < 2 || endShape2.size() < 2) {
     613             :             return ret;
     614             :         }
     615      160192 :         const double angle = GeomHelper::angleDiff(begShape2.angleAt2D(-2), endShape2.angleAt2D(0));
     616             :         NBNode* recordError = nullptr;
     617             : #ifdef DEBUG_SMOOTH_CORNERS
     618             :         if (DEBUGCOND) {
     619             :             std::cout << "   angle=" << RAD2DEG(angle) << "\n";
     620             :         }
     621             :         recordError = const_cast<NBNode*>(&myNode);
     622             : #endif
     623             :         // fill highly acute corners
     624             :         //if (fabs(angle) > DEG2RAD(135)) {
     625             :         //    return ret;
     626             :         //}
     627      160192 :         PositionVector curve = myNode.computeSmoothShape(begShape2, endShape2, cornerDetail + 2, false, 25, 25, recordError, NBNode::AVOID_WIDE_LEFT_TURN);
     628             :         //PositionVector curve = myNode.computeSmoothShape(begShape2, endShape2, cornerDetail + 2, false, 25, 25, recordError, 0);
     629      160192 :         const double curvature = curve.length2D() / MAX2(NUMERICAL_EPS, begPoint.distanceTo2D(endPoint));
     630             : #ifdef DEBUG_SMOOTH_CORNERS
     631             :         if (DEBUGCOND) {
     632             :             std::cout << "   curve=" << curve << " curveLength=" << curve.length2D() << " dist=" << begPoint.distanceTo2D(endPoint) << " curvature=" << curvature << "\n";
     633             :         }
     634             : #endif
     635      160192 :         if (curvature > 2 && angle > DEG2RAD(85)) {
     636             :             // simplify dubious inside corner shape
     637             :             return ret;
     638             :         }
     639      158067 :         if (curve.size() > 2) {
     640             :             curve.erase(curve.begin());
     641             :             curve.pop_back();
     642             :             ret = curve;
     643             :         }
     644      185409 :     }
     645             :     return ret;
     646           0 : }
     647             : 
     648             : void
     649       79477 : NBNodeShapeComputer::computeEdgeBoundaries(const EdgeVector& edges,
     650             :         GeomsMap& geomsCCW,
     651             :         GeomsMap& geomsCW) {
     652             :     // compute boundary lines and extend it by EXT m
     653      407716 :     for (NBEdge* const edge : edges) {
     654             :         // store current edge's boundary as current ccw/cw boundary
     655             :         try {
     656      656478 :             geomsCCW[edge] = edge->getCCWBoundaryLine(myNode);
     657           0 :         } catch (InvalidArgument& e) {
     658           0 :             WRITE_WARNING("While computing intersection geometry at junction '" + myNode.getID() + "': " + std::string(e.what()));
     659           0 :             geomsCCW[edge] = edge->getGeometry();
     660           0 :         }
     661             :         try {
     662      656478 :             geomsCW[edge] = edge->getCWBoundaryLine(myNode);
     663           0 :         } catch (InvalidArgument& e) {
     664           0 :             WRITE_WARNING("While computing intersection geometry at junction '" + myNode.getID() + "': " + std::string(e.what()));
     665           0 :             geomsCW[edge] = edge->getGeometry();
     666           0 :         }
     667             :         // ensure the boundary is valid
     668      328239 :         if (geomsCCW[edge].length2D() < NUMERICAL_EPS) {
     669           0 :             geomsCCW[edge] = edge->getGeometry();
     670             :         }
     671      328239 :         if (geomsCW[edge].length2D() < NUMERICAL_EPS) {
     672           0 :             geomsCW[edge] = edge->getGeometry();
     673             :         }
     674             :         // cut off all parts beyond EXT to avoid issues with curved-back roads
     675      656490 :         geomsCCW[edge] = geomsCCW[edge].getSubpart2D(0, MAX2(EXT, edge->getTotalWidth()));
     676      656490 :         geomsCW[edge] = geomsCW[edge].getSubpart2D(0, MAX2(EXT, edge->getTotalWidth()));
     677             :         // extend the boundary by extrapolating it by EXT m towards the junction
     678      328239 :         geomsCCW[edge].extrapolate2D(EXT, true);
     679      328239 :         geomsCW[edge].extrapolate2D(EXT, true);
     680             :         // ensure minimum length by extending it away from the junction
     681      328239 :         geomsCCW[edge].extrapolate(EXT2, false, true);
     682      328239 :         geomsCW[edge].extrapolate(EXT2, false, true);
     683             :     }
     684       79477 : }
     685             : 
     686             : void
     687       79477 : NBNodeShapeComputer::joinSameDirectionEdges(const EdgeVector& edges, std::map<NBEdge*, std::set<NBEdge*> >& same, bool useEndpoints) {
     688             :     // compute same (edges where an intersection doesn't work well
     689             :     // (always check an edge and its cw neighbor)
     690             :     const double angleChangeLookahead = 35; // distance to look ahead for a misleading angle
     691      228324 :     const bool isXodr = OptionsCont::getOptions().exists("opendrive-files") && OptionsCont::getOptions().isSet("opendrive-files");
     692             :     EdgeSet foundOpposite;
     693      407716 :     for (EdgeVector::const_iterator i = edges.begin(); i != edges.end(); i++) {
     694             :         EdgeVector::const_iterator j;
     695      328239 :         if (i == edges.end() - 1) {
     696             :             j = edges.begin();
     697             :         } else {
     698             :             j = i + 1;
     699             :         }
     700          16 :         if (useEndpoints
     701          16 :                 && !(*i)->hasDefaultGeometryEndpointAtNode(&myNode)
     702      328255 :                 && !(*j)->hasDefaultGeometryEndpointAtNode(&myNode)) {
     703          16 :             continue;
     704             :         }
     705      328223 :         const bool incoming = (*i)->getToNode() == &myNode;
     706      328223 :         const bool incoming2 = (*j)->getToNode() == &myNode;
     707             :         const bool differentDirs = (incoming != incoming2);
     708      382178 :         const bool sameGeom = (*i)->getGeometry() == (differentDirs ? (*j)->getGeometry().reverse() : (*j)->getGeometry());
     709      328223 :         const PositionVector g1 = incoming ? (*i)->getCCWBoundaryLine(myNode) : (*i)->getCWBoundaryLine(myNode);
     710      328223 :         const PositionVector g2 = incoming ? (*j)->getCCWBoundaryLine(myNode) : (*j)->getCWBoundaryLine(myNode);
     711      328223 :         const double angle1further = (g1.size() > 2 && g1[0].distanceTo2D(g1[1]) < angleChangeLookahead ?
     712      328223 :                                       g1.angleAt2D(1) : g1.angleAt2D(0));
     713      328223 :         const double angle2further = (g2.size() > 2 && g2[0].distanceTo2D(g2[1]) < angleChangeLookahead ?
     714      328223 :                                       g2.angleAt2D(1) : g2.angleAt2D(0));
     715      328223 :         const double angleDiff = GeomHelper::angleDiff(g1.angleAt2D(0), g2.angleAt2D(0));
     716      328223 :         const double angleDiffFurther = GeomHelper::angleDiff(angle1further, angle2further);
     717      328223 :         const bool ambiguousGeometry = ((angleDiff > 0 && angleDiffFurther < 0) || (angleDiff < 0 && angleDiffFurther > 0));
     718             :         //if (ambiguousGeometry) {
     719             :         //    @todo: this warning would be helpful in many cases. However, if angle and angleFurther jump between 179 and -179 it is misleading
     720             :         //    WRITE_WARNINGF(TL("Ambiguous angles at junction '%' for edges '%' and '%'."), myNode.getID(), (*i)->getID(), (*j)->getID());
     721             :         //}
     722             : #ifdef DEBUG_NODE_SHAPE
     723             :         if (DEBUGCOND) {
     724             :             std::cout << "   checkSameDirection " << (*i)->getID() << " " << (*j)->getID()
     725             :                       << " diffDirs=" << differentDirs
     726             :                       << " isOpposite=" << (differentDirs && foundOpposite.count(*i) == 0)
     727             :                       << " angleDiff=" << angleDiff
     728             :                       << " ambiguousGeometry=" << ambiguousGeometry
     729             :                       << " badInsersection=" << badIntersection(*i, *j, EXT)
     730             :                       << "\n";
     731             : 
     732             :         }
     733             : #endif
     734      328223 :         if (sameGeom || fabs(angleDiff) < DEG2RAD(20)) {
     735      102670 :             const bool isOpposite = differentDirs && foundOpposite.count(*i) == 0;
     736             :             if (isOpposite) {
     737             :                 foundOpposite.insert(*i);
     738             :                 foundOpposite.insert(*j);
     739             :             }
     740      102670 :             if (isOpposite || ambiguousGeometry || (!isXodr && badIntersection(*i, *j, EXT))) {
     741             :                 // maintain equivalence relation for all members of the equivalence class
     742      106415 :                 for (std::set<NBEdge*>::iterator k = same[*i].begin(); k != same[*i].end(); ++k) {
     743        6193 :                     if (*j != *k) {
     744        6133 :                         same[*k].insert(*j);
     745        6133 :                         same[*j].insert(*k);
     746             :                     }
     747             :                 }
     748      106700 :                 for (std::set<NBEdge*>::iterator k = same[*j].begin(); k != same[*j].end(); ++k) {
     749        6478 :                     if (*i != *k) {
     750        6418 :                         same[*k].insert(*i);
     751        6418 :                         same[*i].insert(*k);
     752             :                     }
     753             :                 }
     754      100222 :                 same[*i].insert(*j);
     755      100222 :                 same[*j].insert(*i);
     756             : #ifdef DEBUG_NODE_SHAPE
     757             :                 if (DEBUGCOND) {
     758             :                     std::cout << "   joinedSameDirectionEdges " << (*i)->getID() << "   " << (*j)->getID() << " isOpposite=" << isOpposite << " ambiguousGeometry=" << ambiguousGeometry << "\n";
     759             :                 }
     760             : #endif
     761             :             }
     762             :         }
     763      328223 :     }
     764       79477 : }
     765             : 
     766             : 
     767             : bool
     768        6776 : NBNodeShapeComputer::badIntersection(const NBEdge* e1, const NBEdge* e2, double distance) {
     769             :     // check whether the two edges are on top of each other. In that case they should be joined
     770             :     // also, if they never touch along their common length
     771        6776 :     const double commonLength = MIN3(distance, e1->getGeometry().length(), e2->getGeometry().length());
     772             :     PositionVector geom1 = e1->getGeometry();
     773             :     PositionVector geom2 = e2->getGeometry();
     774             :     // shift to make geom the centerline of the edge regardless of spreadtype
     775        6776 :     if (e1->getLaneSpreadFunction() == LaneSpreadFunction::RIGHT) {
     776        1251 :         geom1.move2side(e1->getTotalWidth() / 2);
     777             :     }
     778        6776 :     if (e2->getLaneSpreadFunction() == LaneSpreadFunction::RIGHT) {
     779        1214 :         geom2.move2side(e2->getTotalWidth() / 2);
     780             :     }
     781             :     // always let geometry start at myNode
     782        6776 :     if (e1->getToNode() == &myNode) {
     783        5588 :         geom1 = geom1.reverse();
     784             :     }
     785        6776 :     if (e2->getToNode() == &myNode) {
     786        7082 :         geom2 = geom2.reverse();
     787             :     }
     788       13552 :     geom1 = geom1.getSubpart2D(0, commonLength);
     789       13552 :     geom2 = geom2.getSubpart2D(0, commonLength);
     790             :     double endAngleDiff = 0;
     791        6776 :     if (geom1.size() >= 2 && geom2.size() >= 2) {
     792        6776 :         endAngleDiff = fabs(RAD2DEG(GeomHelper::angleDiff(
     793             :                                         geom1.angleAt2D((int)geom1.size() - 2),
     794             :                                         geom2.angleAt2D((int)geom2.size() - 2))));
     795             :     }
     796        6776 :     const double minDistanceThreshold = (e1->getTotalWidth() + e2->getTotalWidth()) / 2 + POSITION_EPS;
     797        6776 :     std::vector<double> distances = geom1.distances(geom2, true);
     798        6776 :     std::vector<double> distances2 = geom1.distances(geom2);
     799             :     const double minDist = VectorHelper<double>::minValue(distances2);
     800             :     const double maxDist = VectorHelper<double>::maxValue(distances);
     801        6776 :     const bool curvingTowards = geom1[0].distanceTo2D(geom2[0]) > minDistanceThreshold && minDist < minDistanceThreshold;
     802        6776 :     const bool onTop = (maxDist - POSITION_EPS < minDistanceThreshold) && endAngleDiff < 30;
     803        6776 :     const bool bothDefault = e1->hasDefaultGeometryEndpointAtNode(&myNode) && e2->hasDefaultGeometryEndpointAtNode(&myNode);
     804        6776 :     const bool neverTouch = minDist > minDistanceThreshold * 2 && !bothDefault;
     805        6776 :     geom1.extrapolate2D(EXT);
     806        6776 :     geom2.extrapolate2D(EXT);
     807        6776 :     Position intersect = geom1.intersectionPosition2D(geom2);
     808        5012 :     const bool intersects = intersect != Position::INVALID && geom1.distance2D(intersect) < POSITION_EPS;
     809             : #ifdef DEBUG_NODE_SHAPE
     810             :     if (DEBUGCOND) {
     811             :         std::cout << "    badIntersect: onTop=" << onTop << " curveTo=" << curvingTowards << " intersects=" << intersects
     812             :                   << " endAngleDiff=" << endAngleDiff
     813             :                   << " geom1=" << geom1 << " geom2=" << geom2
     814             :                   << " distances=" << toString(distances) << " minDist=" << minDist << " maxDist=" << maxDist << " thresh=" << minDistanceThreshold
     815             :                   << " neverTouch=" << neverTouch
     816             :                   << " intersectPos=" << intersect
     817             :                   << "\n";
     818             :     }
     819             : #endif
     820       13552 :     return onTop || curvingTowards || !intersects || neverTouch;
     821        6776 : }
     822             : 
     823             : 
     824             : EdgeVector
     825       79477 : NBNodeShapeComputer::computeUniqueDirectionList(
     826             :     const EdgeVector& all,
     827             :     std::map<NBEdge*, std::set<NBEdge*> >& same,
     828             :     GeomsMap& geomsCCW,
     829             :     GeomsMap& geomsCW) {
     830             :     // store relationships
     831       79477 :     EdgeVector newAll = all;
     832      407716 :     for (NBEdge* e1 : all) {
     833             :         // determine which of the edges marks the outer boundary
     834      328239 :         auto e2NewAll = std::find(newAll.begin(), newAll.end(), e1);
     835             : #ifdef DEBUG_NODE_SHAPE
     836             :         if (DEBUGCOND) std::cout << "computeUniqueDirectionList e1=" << e1->getID()
     837             :                                      << " deleted=" << (e2NewAll == newAll.end())
     838             :                                      << " same=" << joinNamedToStringSorting(same[e1], ',') << "\n";
     839             : #endif
     840      328239 :         if (e2NewAll == newAll.end()) {
     841      100162 :             continue;
     842             :         }
     843      228077 :         auto e1It = std::find(all.begin(), all.end(), e1);
     844             :         auto bestCCW = e1It;
     845             :         auto bestCW = e1It;
     846             :         bool changed = true;
     847      554506 :         while (changed) {
     848             :             changed = false;
     849      532705 :             for (NBEdge* e2 : same[e1]) {
     850             : #ifdef DEBUG_NODE_SHAPE
     851             :                 if (DEBUGCOND) {
     852             :                     std::cout << "  e2=" << e2->getID() << "\n";
     853             :                 }
     854             : #endif
     855      206276 :                 auto e2It = std::find(all.begin(), all.end(), e2);
     856      206276 :                 if (e2It + 1 == bestCCW || (e2It == (all.end() - 1) && bestCCW == all.begin())) {
     857             :                     bestCCW = e2It;
     858             :                     changed = true;
     859             : #ifdef DEBUG_NODE_SHAPE
     860             :                     if (DEBUGCOND) {
     861             :                         std::cout << "    bestCCW=" << e2->getID() << "\n";
     862             :                     }
     863             : #endif
     864      203973 :                 } else if (bestCW + 1 == e2It || (bestCW == (all.end() - 1) && e2It == all.begin())) {
     865             :                     bestCW = e2It;
     866             :                     changed = true;
     867             : #ifdef DEBUG_NODE_SHAPE
     868             :                     if (DEBUGCOND) {
     869             :                         std::cout << "    bestCW=" << e2->getID() << "\n";
     870             :                     }
     871             : #endif
     872             :                 }
     873             :             }
     874             :         }
     875      228077 :         if (bestCW != e1It) {
     876       95115 :             geomsCW[e1] = geomsCW[*bestCW];
     877       95115 :             computeSameEnd(geomsCW[e1], geomsCCW[e1]);
     878             :         }
     879      228077 :         if (bestCCW != e1It) {
     880        2074 :             geomsCCW[e1] = geomsCCW[*bestCCW];
     881        2074 :             computeSameEnd(geomsCW[e1], geomsCCW[e1]);
     882             :         }
     883             :         // clean up
     884      328239 :         for (NBEdge* e2 : same[e1]) {
     885      100162 :             auto e2NewAllIt = std::find(newAll.begin(), newAll.end(), e2);
     886      100162 :             if (e2NewAllIt != newAll.end()) {
     887             :                 newAll.erase(e2NewAllIt);
     888             :             }
     889             :         }
     890             :     }
     891             : #ifdef DEBUG_NODE_SHAPE
     892             :     if (DEBUGCOND) {
     893             :         std::cout << "  newAll:\n";
     894             :         for (NBEdge* e : newAll) {
     895             :             std::cout << "    " << e->getID() << " geomCCW=" << geomsCCW[e] << " geomsCW=" << geomsCW[e] << "\n";
     896             :         }
     897             :     }
     898             : #endif
     899       79477 :     return newAll;
     900             : }
     901             : 
     902             : 
     903             : void
     904      227972 : NBNodeShapeComputer::initNeighbors(const EdgeVector& edges, const EdgeVector::const_iterator& current,
     905             :                                    GeomsMap& geomsCW,
     906             :                                    GeomsMap& geomsCCW,
     907             :                                    EdgeVector::const_iterator& cwi,
     908             :                                    EdgeVector::const_iterator& ccwi,
     909             :                                    double& cad,
     910             :                                    double& ccad) {
     911             :     const double twoPI = (double)(2 * M_PI);
     912      227972 :     cwi = current;
     913             :     cwi++;
     914      227972 :     if (cwi == edges.end()) {
     915       79372 :         std::advance(cwi, -((int)edges.size())); // set to edges.begin();
     916             :     }
     917      227972 :     ccwi = current;
     918      227972 :     if (ccwi == edges.begin()) {
     919             :         std::advance(ccwi, edges.size() - 1); // set to edges.end() - 1;
     920             :     } else {
     921             :         ccwi--;
     922             :     }
     923             : 
     924      227972 :     const double angleCurCCW = geomsCCW[*current].angleAt2D(0);
     925      227972 :     const double angleCurCW = geomsCW[*current].angleAt2D(0);
     926      227972 :     const double angleCCW = geomsCW[*ccwi].angleAt2D(0);
     927      227972 :     const double angleCW = geomsCCW[*cwi].angleAt2D(0);
     928      227972 :     ccad = angleCCW - angleCurCCW;
     929      309813 :     while (ccad < 0.) {
     930       81841 :         ccad += twoPI;
     931             :     }
     932      227972 :     cad = angleCurCW - angleCW;
     933      309813 :     while (cad < 0.) {
     934       81841 :         cad += twoPI;
     935             :     }
     936      227972 : }
     937             : 
     938             : 
     939             : 
     940             : const PositionVector
     941       27226 : NBNodeShapeComputer::computeNodeShapeSmall() {
     942             : #ifdef DEBUG_NODE_SHAPE
     943             :     if (DEBUGCOND) {
     944             :         std::cout << "computeNodeShapeSmall node=" << myNode.getID() << "\n";
     945             :     }
     946             : #endif
     947       27226 :     PositionVector ret;
     948       68490 :     for (NBEdge* e : myNode.getEdges()) {
     949             :         // compute crossing with normal
     950       41264 :         PositionVector edgebound1 = e->getCCWBoundaryLine(myNode).getSubpartByIndex(0, 2);
     951       41264 :         PositionVector edgebound2 = e->getCWBoundaryLine(myNode).getSubpartByIndex(0, 2);
     952       41264 :         Position delta = edgebound1[1] - edgebound1[0];
     953       41264 :         delta.set(-delta.y(), delta.x()); // rotate 90 degrees
     954       41264 :         PositionVector cross(myNode.getPosition(), myNode.getPosition() + delta);
     955       41264 :         cross.extrapolate2D(500);
     956       41264 :         edgebound1.extrapolate2D(500);
     957       41264 :         edgebound2.extrapolate2D(500);
     958       41264 :         if (cross.intersects(edgebound1)) {
     959       41264 :             Position np = cross.intersectionPosition2D(edgebound1);
     960       41264 :             np.set(np.x(), np.y(), myNode.getPosition().z());
     961       41264 :             ret.push_back_noDoublePos(np);
     962             :         }
     963       41264 :         if (cross.intersects(edgebound2)) {
     964       41260 :             Position np = cross.intersectionPosition2D(edgebound2);
     965       41260 :             np.set(np.x(), np.y(), myNode.getPosition().z());
     966       41260 :             ret.push_back_noDoublePos(np);
     967             :         }
     968       41264 :         e->resetNodeBorder(&myNode);
     969       41264 :     }
     970       27226 :     return ret;
     971           0 : }
     972             : 
     973             : 
     974             : double
     975       79477 : NBNodeShapeComputer::getDefaultRadius(const OptionsCont& oc) {
     976             :     // look for incoming/outgoing edge pairs that do not go straight and allow wide vehicles
     977             :     // (connection information is not available yet)
     978             :     // @TODO compute the radius for each pair of neighboring edge intersections in computeNodeShapeDefault rather than use the maximum
     979       79477 :     const double radius = oc.getFloat("default.junctions.radius");
     980       79477 :     const double smallRadius = oc.getFloat("junctions.small-radius");
     981             :     double maxRightAngle = 0; // rad
     982             :     double extraWidthRight = 0; // m
     983             :     double maxLeftAngle = 0; // rad
     984             :     double extraWidthLeft = 0; // m
     985             :     int laneDelta = 0;
     986             :     int totalWideLanesIn = 0;
     987      242902 :     for (NBEdge* in : myNode.getIncomingEdges()) {
     988             :         int wideLanesIn = 0;
     989      372963 :         for (int i = 0; i < in->getNumLanes(); i++) {
     990      209538 :             if ((in->getPermissions(i) & SVC_LARGE_TURN) != 0) {
     991      129860 :                 wideLanesIn++;
     992             :             }
     993             :         }
     994      163425 :         totalWideLanesIn += wideLanesIn;
     995      569125 :         for (NBEdge* out : myNode.getOutgoingEdges()) {
     996      405700 :             if ((in->getPermissions() & out->getPermissions() & SVC_LARGE_TURN) != 0) {
     997      244480 :                 if (myNode.getDirection(in, out) == LinkDirection::TURN) {
     998       68152 :                     continue;
     999             :                 };
    1000      176328 :                 const double angle = GeomHelper::angleDiff(
    1001      176328 :                                          in->getGeometry().angleAt2D(-2),
    1002             :                                          out->getGeometry().angleAt2D(0));
    1003      176328 :                 if (angle < 0) {
    1004       82001 :                     if (maxRightAngle < -angle) {
    1005             :                         maxRightAngle = -angle;
    1006       51710 :                         extraWidthRight = MAX2(getExtraWidth(in, SVC_LARGE_TURN), getExtraWidth(out, SVC_LARGE_TURN));
    1007             :                     }
    1008             :                 } else {
    1009       94327 :                     if (maxLeftAngle < angle) {
    1010             :                         maxLeftAngle = angle;
    1011             :                         // all edges clockwise between in and out count as extra width
    1012             :                         extraWidthLeft = 0;
    1013       45210 :                         EdgeVector::const_iterator pIn = std::find(myNode.getEdges().begin(), myNode.getEdges().end(), in);
    1014       45210 :                         NBContHelper::nextCW(myNode.getEdges(), pIn);
    1015      118886 :                         while (*pIn != out) {
    1016       73676 :                             extraWidthLeft += (*pIn)->getTotalWidth();
    1017             : #ifdef DEBUG_RADIUS
    1018             :                             if (DEBUGCOND) {
    1019             :                                 std::cout << "   in=" << in->getID() << " out=" << out->getID() << " extra=" << (*pIn)->getID() << " extraWidthLeft=" << extraWidthLeft << "\n";
    1020             :                             }
    1021             : #endif
    1022       73676 :                             NBContHelper::nextCW(myNode.getEdges(), pIn);
    1023             :                         }
    1024             :                     }
    1025             :                 }
    1026             :                 int wideLanesOut = 0;
    1027      426033 :                 for (int i = 0; i < out->getNumLanes(); i++) {
    1028      249705 :                     if ((out->getPermissions(i) & SVC_LARGE_TURN) != 0) {
    1029      232473 :                         wideLanesOut++;
    1030             :                     }
    1031             :                 }
    1032             : #ifdef DEBUG_RADIUS
    1033             :                 if (DEBUGCOND) {
    1034             :                     std::cout << "   in=" << in->getID() << " out=" << out->getID() << " wideLanesIn=" << wideLanesIn << " wideLanesOut=" << wideLanesOut << "\n";
    1035             :                 }
    1036             : #endif
    1037      176328 :                 laneDelta = MAX2(laneDelta, abs(wideLanesOut - wideLanesIn));
    1038             :             }
    1039             :         }
    1040             :     }
    1041             :     // special case: on/off-ramp
    1042       79477 :     if (myNode.getOutgoingEdges().size() == 1 || myNode.getIncomingEdges().size() == 1) {
    1043             :         int totalWideLanesOut = 0;
    1044       78621 :         for (NBEdge* out : myNode.getOutgoingEdges()) {
    1045      106129 :             for (int i = 0; i < out->getNumLanes(); i++) {
    1046       61176 :                 if ((out->getPermissions(i) & SVC_LARGE_TURN) != 0) {
    1047       31382 :                     totalWideLanesOut++;
    1048             :                 }
    1049             :             }
    1050             :         }
    1051       33668 :         if (totalWideLanesIn == totalWideLanesOut) {
    1052             :             // use total laneDelta instead of individual edge lane delta
    1053             :             laneDelta = 0;
    1054             :         }
    1055             :     }
    1056             :     // changing the number of wide-vehicle lanes on a straight segment requires a larger junction to allow for smooth driving
    1057             :     // otherwise we can reduce the radius according to the angle
    1058             :     double result = radius;
    1059             :     // left turns are assumed to cross additional edges and thus du not determine the required radius in most cases
    1060             :     double maxTurnAngle = maxRightAngle;
    1061             :     double extraWidth = extraWidthRight;
    1062       79477 :     if (maxRightAngle < DEG2RAD(5)) {
    1063             :         maxTurnAngle = maxLeftAngle;
    1064             :         extraWidth = extraWidthLeft;
    1065             :     }
    1066       79477 :     const double minRadius = maxTurnAngle >= DEG2RAD(30) ? MIN2(smallRadius, radius) : smallRadius;
    1067       79477 :     if (laneDelta == 0 || maxTurnAngle >= DEG2RAD(30) || myNode.isConstantWidthTransition()) {
    1068             :         // subtract radius gained from extra lanes
    1069             :         // do not increase radius for turns that are sharper than a right angle
    1070       75672 :         result = radius * tan(0.5 * MIN2(0.5 * M_PI, maxTurnAngle)) - extraWidth;
    1071             :     }
    1072             :     result = MAX2(minRadius, result);
    1073             : #ifdef DEBUG_RADIUS
    1074             :     if (DEBUGCOND) {
    1075             :         std::cout << "getDefaultRadius n=" << myNode.getID()
    1076             :                   << " r=" << radius << " sr=" << smallRadius
    1077             :                   << " mr=" << minRadius
    1078             :                   << " laneDelta=" << laneDelta
    1079             :                   << " rightA=" << RAD2DEG(maxRightAngle)
    1080             :                   << " leftA=" << RAD2DEG(maxLeftAngle)
    1081             :                   << " maxA=" << RAD2DEG(maxTurnAngle)
    1082             :                   << " extraWidth=" << extraWidth
    1083             :                   << " result=" << result << "\n";
    1084             :     }
    1085             : #endif
    1086       79477 :     return result;
    1087             : }
    1088             : 
    1089             : 
    1090             : bool
    1091      189862 : NBNodeShapeComputer::isDivided(const NBEdge* e, std::set<NBEdge*> same, const PositionVector& ccw, const PositionVector& cw) const {
    1092      189862 :     if (same.size() < 2) {
    1093             :         return false;
    1094             :     }
    1095             :     std::set<Position> endPoints;
    1096        2606 :     endPoints.insert(e->getEndpointAtNode(&myNode));
    1097        5128 :     for (NBEdge* s : same) {
    1098        7650 :         endPoints.insert(s->getEndpointAtNode(&myNode));
    1099             :     }
    1100        1303 :     if (endPoints.size() > 1) {
    1101        1082 :         std::vector<double> distances = ccw.distances(cw, true);
    1102        1082 :         double width = e->getTotalWidth();
    1103        4380 :         for (const NBEdge* e2 : same) {
    1104        3298 :             width += e2->getTotalWidth();
    1105             :         }
    1106             :         const double maxDist = VectorHelper<double>::maxValue(distances);
    1107        1082 :         const double maxDivider = maxDist - width;
    1108        1082 :         return maxDivider >= 5;
    1109             :     }
    1110             :     return false;
    1111             : }
    1112             : 
    1113             : 
    1114             : double
    1115      103420 : NBNodeShapeComputer::getExtraWidth(const NBEdge* e, SVCPermissions exclude) {
    1116             :     double result = 0;
    1117             :     int lane = 0;
    1118      103420 :     while (lane < e->getNumLanes() && e->getPermissions(lane) == 0) {
    1119             :         // ignore forbidden lanes out the outside
    1120           0 :         lane++;
    1121             :     }
    1122      115380 :     while (lane < e->getNumLanes() && (e->getPermissions(lane) & exclude) == 0) {
    1123       11960 :         result += e->getLaneWidth(lane);
    1124       11960 :         lane++;
    1125             :     }
    1126      103420 :     return result;
    1127             : }
    1128             : 
    1129             : 
    1130             : double
    1131         440 : NBNodeShapeComputer::divisionWidth(const NBEdge* e, std::set<NBEdge*> same, const Position& p, const Position& p2) {
    1132             :     double result = p.distanceTo2D(p2);
    1133         440 :     result -= e->getTotalWidth();
    1134        1962 :     for (NBEdge* e2 : same) {
    1135        1522 :         result -= e2->getTotalWidth();
    1136             :     }
    1137         440 :     return MAX2(0.0, result);
    1138             : }
    1139             : 
    1140             : /****************************************************************************/

Generated by: LCOV version 1.14