LCOV - code coverage report
Current view: top level - src/netbuild - NBNodeShapeComputer.cpp (source / functions) Coverage Total Hit
Test: lcov.info Lines: 91.3 % 461 421
Test Date: 2025-11-13 15:38:19 Functions: 100.0 % 18 18

            Line data    Source code
       1              : /****************************************************************************/
       2              : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
       3              : // Copyright (C) 2001-2025 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        78199 : NBNodeShapeComputer::NBNodeShapeComputer(const NBNode& node) :
      55        78199 :     myNode(node),
      56        78199 :     myRadius(node.getRadius()) {
      57        78199 :     if (node.getEdges().size() > 4 && !NBNodeTypeComputer::isRailwayNode(&node)) {
      58        19547 :         EXT = 50;
      59              :     } else {
      60        58652 :         EXT = 100;
      61              :     }
      62        78199 : }
      63              : 
      64              : 
      65        78199 : NBNodeShapeComputer::~NBNodeShapeComputer() {}
      66              : 
      67              : 
      68              : const PositionVector
      69        78199 : 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        78199 :     if (myNode.getEdges().size() == 1 || forceSmall) {
      83         9940 :         return computeNodeShapeSmall();
      84              :     }
      85        68259 :     if (myNode.getEdges().size() == 2 && myNode.getIncomingEdges().size() == 1) {
      86        21827 :         if (myNode.getIncomingEdges()[0]->isTurningDirectionAt(myNode.getOutgoingEdges()[0])) {
      87        10971 :             return computeNodeShapeSmall();
      88              :         }
      89              :     }
      90        57288 :     const bool geometryLike = myNode.isSimpleContinuation(true, true);
      91        57288 :     const PositionVector& ret = computeNodeShapeDefault(geometryLike);
      92              :     // fail fall-back: use "computeNodeShapeSmall"
      93        57288 :     if (ret.size() < 3) {
      94           83 :         return computeNodeShapeSmall();
      95              :     }
      96              :     return ret;
      97        57288 : }
      98              : 
      99              : 
     100              : void
     101        75552 : NBNodeShapeComputer::computeSameEnd(PositionVector& l1, PositionVector& l2) {
     102              :     assert(l1[0].distanceTo2D(l1[1]) >= EXT);
     103              :     assert(l2[0].distanceTo2D(l2[1]) >= EXT);
     104        75552 :     PositionVector tmp;
     105        75552 :     tmp.push_back(PositionVector::positionAtOffset2D(l1[0], l1[1], EXT));
     106        75552 :     tmp.push_back(l1[1]);
     107        75552 :     tmp[1].sub(tmp[0]);
     108        75552 :     tmp[1].set(-tmp[1].y(), tmp[1].x());
     109        75552 :     tmp[1].add(tmp[0]);
     110        75552 :     tmp.extrapolate2D(EXT);
     111        75552 :     if (l2.intersects(tmp[0], tmp[1])) {
     112        75489 :         const double offset = l2.intersectsAtLengths2D(tmp)[0];
     113        75489 :         if (l2.length2D() - offset > POSITION_EPS) {
     114        75489 :             PositionVector tl2 = l2.getSubpart2D(offset, l2.length2D());
     115        75489 :             tl2.extrapolate2D(EXT);
     116        75489 :             l2.erase(l2.begin(), l2.begin() + (l2.size() - tl2.size()));
     117        75489 :             l2[0] = tl2[0];
     118        75489 :         }
     119              :     }
     120        75552 : }
     121              : 
     122              : 
     123              : const PositionVector
     124        57288 : NBNodeShapeComputer::computeNodeShapeDefault(bool simpleContinuation) {
     125              :     // if we have less than two edges, we can not compute the node's shape this way
     126        57288 :     if (myNode.getEdges().size() < 2) {
     127            0 :         return PositionVector();
     128              :     }
     129              :     // magic values
     130        57288 :     const OptionsCont& oc = OptionsCont::getOptions();
     131        57288 :     const double defaultRadius = getDefaultRadius(oc);
     132        57288 :     const bool useDefaultRadius = myNode.getRadius() == NBNode::UNSPECIFIED_RADIUS || myNode.getRadius() == defaultRadius;
     133        57288 :     myRadius = (useDefaultRadius ? defaultRadius : myNode.getRadius());
     134        57288 :     double smallRadius = useDefaultRadius ? oc.getFloat("junctions.small-radius") : myRadius;
     135        57288 :     const int cornerDetail = oc.getInt("junctions.corner-detail");
     136        57288 :     const double sCurveStretch = oc.getFloat("junctions.scurve-stretch");
     137        57288 :     const bool useEndpoints = oc.getBool("junctions.endpoint-shape");
     138        57288 :     const bool rectangularCut = oc.getBool("rectangular-lane-cut");
     139        57288 :     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       109903 :     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        57288 :     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*, ComparatorIdLess> > 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        57288 :     EdgeVector usedEdges = myNode.getEdges();
     164        57288 :     computeEdgeBoundaries(usedEdges, geomsCCW, geomsCW);
     165              : 
     166              :     // check which edges are parallel
     167        57288 :     joinSameDirectionEdges(usedEdges, same, useEndpoints);
     168              :     // compute unique direction list
     169        57288 :     EdgeVector newAll = computeUniqueDirectionList(usedEdges, same, geomsCCW, geomsCW);
     170              :     // if we have only two "directions", let's not compute the geometry using this method
     171        57288 :     if (newAll.size() < 2) {
     172           83 :         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       221739 :     for (i = newAll.begin(); i != newAll.end(); ++i) {
     185       164534 :         EdgeVector::const_iterator cwi = i;
     186       164534 :         EdgeVector::const_iterator ccwi = i;
     187              :         double ccad;
     188              :         double cad;
     189       164534 :         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       164534 :         if (*cwi == *ccwi &&
     196              :                 (
     197              :                     // no change in lane numbers, even low angles still give a good intersection
     198        28990 :                     (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        12210 :                     || (!simpleContinuation && fabs(ccad - cad) < DEG2RAD(22.5)))
     202              :            ) {
     203              :             // compute the mean position between both edges ends ...
     204              :             Position p;
     205        27634 :             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        27634 :                 p = geomsCCW[*ccwi][0];
     216        27634 :                 p.add(geomsCW[*ccwi][0]);
     217        27634 :                 p.add(geomsCCW[*i][0]);
     218        27634 :                 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        55268 :             double dist = MAX2(
     228        27634 :                               geomsCCW[*i].nearest_offset_to_point2D(p),
     229        27634 :                               geomsCW[*i].nearest_offset_to_point2D(p));
     230        27634 :             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        27634 :                 if (!simpleContinuation) {
     261         9484 :                     dist += myRadius;
     262              :                 } else {
     263              :                     // if the angles change, junction should have some size to avoid degenerate shape
     264        18150 :                     double radius2 = fabs(ccad - cad) * (*i)->getNumLanes();
     265        18150 :                     if (radius2 > NUMERICAL_EPS || openDriveOutput) {
     266              :                         radius2 = MAX2(0.15, radius2);
     267              :                     }
     268        18150 :                     if (myNode.getCrossings().size() > 0) {
     269           50 :                         double width = myNode.getCrossings()[0]->customWidth;
     270           50 :                         if (width == NBEdge::UNSPECIFIED_WIDTH) {
     271           56 :                             width = OptionsCont::getOptions().getFloat("default.crossing-width");
     272              :                         }
     273           50 :                         radius2 = MAX2(radius2, width / 2);
     274              :                     }
     275        18150 :                     if (!useDefaultRadius) {
     276            4 :                         radius2 = MAX2(radius2, myRadius);
     277              :                     }
     278        18150 :                     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        27634 :                 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       136900 :             const bool ccwCloser = ccad < cad;
     292       136900 :             const bool cwLargeTurn = needsLargeTurn(*i, *cwi, same);
     293       136900 :             const bool ccwLargeTurn = needsLargeTurn(*i, *ccwi, same);
     294       136900 :             const bool neighLargeTurn = ccwCloser ? ccwLargeTurn : cwLargeTurn;
     295       136900 :             const bool neigh2LargeTurn =  ccwCloser ? cwLargeTurn : ccwLargeTurn;
     296              :             // the border facing the closer neighbor
     297       136900 :             const PositionVector& currGeom = ccwCloser ? geomsCCW[*i] : geomsCW[*i];
     298              :             // the border facing the far neighbor
     299       136900 :             const PositionVector& currGeom2 = ccwCloser ? geomsCW[*i] : geomsCCW[*i];
     300              :             // the border of the closer neighbor
     301       136900 :             const PositionVector& neighGeom = ccwCloser ? geomsCW[*ccwi] : geomsCCW[*cwi];
     302              :             // the border of the far neighbor
     303       136900 :             const PositionVector& neighGeom2 = ccwCloser ? geomsCCW[*cwi] : geomsCW[*ccwi];
     304              :             // whether the current edge/direction spans a divided road
     305       273800 :             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       136900 :             if (!simpleContinuation) {
     316       125994 :                 if (useEndpoints && !(*i)->hasDefaultGeometryEndpointAtNode(&myNode)) {
     317           16 :                     distances[*i] = EXT;
     318       125978 :                 } else if (currGeom.intersects(neighGeom)) {
     319       125425 :                     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       125425 :                     if (*cwi != *ccwi && currGeom2.intersects(neighGeom2)) {
     326              :                         // also use the second intersection point
     327              :                         // but prevent very large node shapes
     328       120062 :                         const double farAngleDist = ccwCloser ? cad : ccad;
     329       120062 :                         double a1 = distances[*i];
     330       120062 :                         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       120062 :                         if (a2 <= EXT) {
     340         4662 :                             if (keepBothDistances) {
     341           51 :                                 if (ccwCloser) {
     342           14 :                                     distances2[*i] = a2;
     343              :                                 } else {
     344           37 :                                     distances[*i] = a2;
     345           37 :                                     distances2[*i] = a1;
     346              :                                 }
     347              :                             } else {
     348         4611 :                                 distances[*i] = MAX2(a1, a2);
     349              :                             }
     350        28506 :                         } else if (ccad > DEG2RAD(90. + 45.) && cad > DEG2RAD(90. + 45.)
     351       119426 :                                    && (fabs(ccad - cad) > DEG2RAD(10)
     352         1174 :                                        || MAX2(ccad, cad) > DEG2RAD(160)
     353          176 :                                        || (a2 - a1) > 7
     354              :                                        // keep roundabouts nodes small
     355          169 :                                        || myNode.isRoundabout())) {
     356              : #ifdef DEBUG_NODE_SHAPE
     357              :                             if (DEBUGCOND) {
     358              :                                 std::cout << "     ignore a2\n";
     359              :                             }
     360              : #endif
     361              :                             // do nothing.
     362       111533 :                         } else if (farAngleDist < DEG2RAD(135) || (fabs(RAD2DEG(farAngleDist) - 180) > 1 && fabs(a2 - a1) < 10)) {
     363        87377 :                             if (keepBothDistances) {
     364          262 :                                 if (ccwCloser) {
     365          129 :                                     distances2[*i] = a2;
     366              :                                 } else {
     367          133 :                                     distances[*i] = a2;
     368          133 :                                     distances2[*i] = a1;
     369              :                                 }
     370              :                             } else {
     371        87115 :                                 distances[*i] = MAX2(a1, a2);
     372              :                             }
     373              :                         }
     374              : #ifdef DEBUG_NODE_SHAPE
     375              :                         if (DEBUGCOND) {
     376              :                             std::cout << "   a1=" << a1 << " a2=" << a2 << " keepBoth=" << keepBothDistances << " dist=" << distances[*i] << "\n";
     377              :                         }
     378              : #endif
     379              :                     }
     380              :                 } else {
     381          553 :                     if (*cwi != *ccwi && currGeom2.intersects(neighGeom2)) {
     382          485 :                         distances[*i] = (neigh2LargeTurn ? myRadius : smallRadius) + currGeom2.intersectsAtLengths2D(neighGeom2)[0];
     383              : #ifdef DEBUG_NODE_SHAPE
     384              :                         if (DEBUGCOND) {
     385              :                             std::cout << "   neigh2 intersects dist=" << distances[*i] << " currGeom2=" << currGeom2 << " neighGeom2=" << neighGeom2 << "\n";
     386              :                         }
     387              : #endif
     388              :                     } else {
     389           68 :                         distances[*i] = EXT + myRadius;
     390              : #ifdef DEBUG_NODE_SHAPE
     391              :                         if (DEBUGCOND) {
     392              :                             std::cout << "   no intersects dist=" << distances[*i]  << " currGeom=" << currGeom << " neighGeom=" << neighGeom << " currGeom2=" << currGeom2 << " neighGeom2=" << neighGeom2 << "\n";
     393              :                         }
     394              : #endif
     395              :                     }
     396              :                 }
     397              :             } else {
     398        10906 :                 if (currGeom.intersects(neighGeom)) {
     399        10878 :                     distances[*i] = currGeom.intersectsAtLengths2D(neighGeom)[0];
     400              :                 } else {
     401           28 :                     distances[*i] = (double) EXT;
     402              :                 }
     403              :             }
     404              :         }
     405       164534 :         if (useDefaultRadius && sCurveStretch > 0) {
     406           32 :             double sCurveWidth = myNode.getDisplacementError();
     407           32 :             if (sCurveWidth > 0) {
     408            2 :                 const double sCurveRadius = myRadius + sCurveWidth / SUMO_const_laneWidth * sCurveStretch * pow((*i)->getSpeed(), 2 + sCurveStretch) / 1000;
     409            2 :                 const double stretch = EXT + sCurveRadius - distances[*i];
     410            2 :                 if (stretch > 0) {
     411            2 :                     distances[*i] += stretch;
     412              :                     // fixate extended geometry for repeated computation
     413            2 :                     const double shorten = distances[*i] - EXT;
     414            2 :                     (*i)->shortenGeometryAtNode(&myNode, shorten);
     415            2 :                     for (std::set<NBEdge*>::iterator k = same[*i].begin(); k != same[*i].end(); ++k) {
     416            0 :                         (*k)->shortenGeometryAtNode(&myNode, shorten);
     417              :                     }
     418              : #ifdef DEBUG_NODE_SHAPE
     419              :                     if (DEBUGCOND) {
     420              :                         std::cout << "   stretching junction: sCurveWidth=" << sCurveWidth << " sCurveRadius=" << sCurveRadius << " stretch=" << stretch << " dist=" << distances[*i]  << "\n";
     421              :                     }
     422              : #endif
     423              :                 }
     424              :             }
     425              :         }
     426              :     }
     427              : 
     428       221739 :     for (NBEdge* const edge : newAll) {
     429       164534 :         if (distances.find(edge) == distances.end()) {
     430              :             assert(false);
     431            0 :             distances[edge] = EXT;
     432              :         }
     433              :     }
     434              :     // because of lane spread right the crossing point may be identical to the junction center and thus the distance is exactly EXT
     435        57205 :     const double off = EXT - NUMERICAL_EPS;
     436              :     // prevent inverted node shapes
     437              :     // (may happen with near-parallel edges)
     438        57205 :     const double minDistSum = 2 * (EXT + myRadius);
     439       221739 :     for (NBEdge* const edge : newAll) {
     440       164534 :         if (distances[edge] < off && edge->hasDefaultGeometryEndpointAtNode(&myNode)) {
     441         7696 :             for (EdgeVector::const_iterator j = newAll.begin(); j != newAll.end(); ++j) {
     442         5756 :                 if (distances[*j] > off && (*j)->hasDefaultGeometryEndpointAtNode(&myNode) && distances[edge] + distances[*j] < minDistSum) {
     443          395 :                     const double angleDiff = fabs(NBHelpers::relAngle(edge->getAngleAtNode(&myNode), (*j)->getAngleAtNode(&myNode)));
     444          395 :                     if (angleDiff > 160 || angleDiff < 20) {
     445              : #ifdef DEBUG_NODE_SHAPE
     446              :                         if (DEBUGCOND) {
     447              :                             std::cout << "   increasing dist for i=" << edge->getID() << " because of j=" << (*j)->getID() << " jDist=" << distances[*j]
     448              :                                       << "  oldI=" << distances[edge] << " newI=" << minDistSum - distances[*j]
     449              :                                       << " angleDiff=" << angleDiff
     450              :                                       << " geomI=" << edge->getGeometry() << " geomJ=" << (*j)->getGeometry() << "\n";
     451              :                         }
     452              : #endif
     453          166 :                         distances[edge] = minDistSum - distances[*j];
     454              :                     }
     455              :                 }
     456              :             }
     457              :         }
     458              :     }
     459              : 
     460              : 
     461              :     // build
     462        57205 :     PositionVector ret;
     463       221739 :     for (i = newAll.begin(); i != newAll.end(); ++i) {
     464       164534 :         const PositionVector& ccwBound = geomsCCW[*i];
     465       164534 :         const PositionVector& cwBound = geomsCW[*i];
     466              :         //double offset = MIN3(distances[*i], cwBound.length2D() - POSITION_EPS, ccwBound.length2D() - POSITION_EPS);
     467       164534 :         double offset = distances[*i];
     468          313 :         double offset2 = distances2.count(*i) != 0 ? distances2[*i] : offset;
     469       164534 :         if (offset != offset2) {
     470              :             // keep rectangular cuts if the difference is small or the roads aren't
     471              :             // really divided by much (unless the angle is very different)
     472          626 :             const double dWidth = divisionWidth(*i, same[*i],
     473          313 :                                                 ccwBound.positionAtOffset2D(offset),
     474          313 :                                                 cwBound.positionAtOffset2D(offset2));
     475          313 :             const double angle = RAD2DEG(GeomHelper::angleDiff(ccwBound.angleAt2D(0), cwBound.angleAt2D(0)));
     476          313 :             const double oDelta = fabs(offset - offset2);
     477              :             //std::cout << " i=" << (*i)->getID() << " offset=" << offset << " offset2=" << offset2 << " dWidth=" << dWidth << " angle=" << angle << " same=" << joinNamedToStringSorting(same[*i], ",") << "\n";
     478          313 :             if ((((oDelta < 5 || dWidth < 10) && fabs(angle) < 30)) || (fabs(angle) < 5 && myNode.getType() != SumoXMLNodeType::RAIL_CROSSING)) {
     479              : #ifdef DEBUG_NODE_SHAPE
     480              :                 std::cout << " i=" << (*i)->getID() << " offset=" << offset << " offset2=" << offset2 << " dWidth=" << dWidth << " angle=" << angle << " same=" << joinNamedToStringSorting(same[*i], ",") << "\n";
     481              : #endif
     482              :                 offset = MAX2(offset, offset2);
     483              :                 offset2 = offset;
     484              :             }
     485              :         }
     486       164534 :         if (!(*i)->hasDefaultGeometryEndpointAtNode(&myNode)) {
     487              :             // for non geometry-endpoints, only shorten but never extend the geometry
     488        16378 :             if (advanceStopLine > 0 && offset < EXT) {
     489              : #ifdef DEBUG_NODE_SHAPE
     490              :                 std::cout << " i=" << (*i)->getID() << " offset=" << offset << " advanceStopLine=" << advanceStopLine << "\n";
     491              : #endif
     492              :                 // fixate extended geometry for repeated computation
     493            0 :                 (*i)->extendGeometryAtNode(&myNode, advanceStopLine);
     494            0 :                 for (std::set<NBEdge*>::iterator k = same[*i].begin(); k != same[*i].end(); ++k) {
     495            0 :                     (*k)->extendGeometryAtNode(&myNode, advanceStopLine);
     496              :                 }
     497              :             }
     498        16378 :             offset = MAX2(EXT - advanceStopLine, offset);
     499              :             offset2 = MAX2(EXT - advanceStopLine, offset2);
     500              :         }
     501       164534 :         if (offset == -1) {
     502            0 :             WRITE_WARNINGF(TL("Fixing offset for edge '%' at node '%."), (*i)->getID(), myNode.getID());
     503              :             offset = -.1;
     504              :             offset2 = -.1;
     505              :         }
     506       164534 :         Position p = ccwBound.positionAtOffset2D(offset);
     507       164534 :         p.setz(myNode.getPosition().z());
     508       164534 :         if (i != newAll.begin()) {
     509       321987 :             ret.append(getSmoothCorner(geomsCW[*(i - 1)], ccwBound, ret[-1], p, cornerDetail));
     510              :         }
     511       164534 :         Position p2 = cwBound.positionAtOffset2D(offset2);
     512       164534 :         p2.setz(myNode.getPosition().z());
     513              :         //ret.append(getEdgeCuts(*i, geomsCCW, geomsCW, offset, offset2, same));
     514       164534 :         ret.push_back_noDoublePos(p);
     515       164534 :         ret.push_back_noDoublePos(p2);
     516              : #ifdef DEBUG_NODE_SHAPE
     517              :         if (DEBUGCOND) {
     518              :             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";
     519              :         }
     520              : #endif
     521       164534 :         (*i)->setNodeBorder(&myNode, p, p2, rectangularCut);
     522       242012 :         for (std::set<NBEdge*>::iterator k = same[*i].begin(); k != same[*i].end(); ++k) {
     523        77478 :             (*k)->setNodeBorder(&myNode, p, p2, rectangularCut);
     524              :         }
     525              :     }
     526              :     // final curve segment
     527       171615 :     ret.append(getSmoothCorner(geomsCW[*(newAll.end() - 1)], geomsCCW[*newAll.begin()], ret[-1], ret[0], cornerDetail));
     528              : #ifdef DEBUG_NODE_SHAPE
     529              :     if (DEBUGCOND) {
     530              :         std::cout << " final shape=" << ret << "\n";
     531              :     }
     532              : #endif
     533              :     return ret;
     534       114493 : }
     535              : 
     536              : 
     537              : double
     538       245487 : NBNodeShapeComputer::closestIntersection(const PositionVector& geom1, const PositionVector& geom2, double offset) {
     539       245487 :     std::vector<double> intersections = geom1.intersectsAtLengths2D(geom2);
     540       245487 :     double result = intersections[0];
     541       260254 :     for (std::vector<double>::iterator it = intersections.begin() + 1; it != intersections.end(); ++it) {
     542        14767 :         if (fabs(*it - offset) < fabs(result - offset)) {
     543              :             result = *it;
     544              :         }
     545              :     }
     546       245487 :     return result;
     547       245487 : }
     548              : 
     549              : bool
     550       273800 : NBNodeShapeComputer::needsLargeTurn(NBEdge* e1, NBEdge* e2,
     551              :                                     std::map<NBEdge*, std::set<NBEdge*, ComparatorIdLess> >& same) const {
     552       273800 :     const SVCPermissions p1 = e1->getPermissions();
     553       273800 :     const SVCPermissions p2 = e2->getPermissions();
     554       273800 :     if ((p1 & p2 & SVC_LARGE_TURN) != 0) {
     555              :         // note: would could also check whether there is actually a connection
     556              :         // between those edges
     557              :         return true;
     558              :     }
     559              :     // maybe edges in the same direction need a large turn
     560       167607 :     for (NBEdge* e2s : same[e2]) {
     561        43648 :         if ((p1 & e2s->getPermissions() & SVC_LARGE_TURN) != 0
     562        43648 :                 && (e1->getToNode() == e2s->getFromNode() || e2s->getToNode() == e1->getFromNode())) {
     563              :             return true;
     564              :         }
     565        74424 :         for (NBEdge* e1s : same[e1]) {
     566        31775 :             if ((e2s->getPermissions() & e1s->getPermissions() & SVC_LARGE_TURN) != 0
     567        31775 :                     && (e2s->getToNode() == e1s->getFromNode() || e1s->getToNode() == e2s->getFromNode())) {
     568              :                 return true;
     569              :             }
     570              :         }
     571              :     }
     572       166065 :     for (NBEdge* e1s : same[e1]) {
     573        42571 :         if ((p2 & e1s->getPermissions() & SVC_LARGE_TURN) != 0
     574        42571 :                 && (e2->getToNode() == e1s->getFromNode() || e1s->getToNode() == e2->getFromNode())) {
     575              :             return true;
     576              :         }
     577              :     }
     578              :     //std::cout << " e1=" << e1->getID() << " e2=" << e2->getID() << " sameE1=" << toString(same[e1]) << " sameE2=" << toString(same[e2]) << "\n";
     579              :     return false;
     580              : }
     581              : 
     582              : PositionVector
     583       164534 : NBNodeShapeComputer::getSmoothCorner(PositionVector begShape, PositionVector endShape,
     584              :                                      const Position& begPoint, const Position& endPoint, int cornerDetail) {
     585       164534 :     PositionVector ret;
     586       164534 :     if (cornerDetail > 0) {
     587       127797 :         PositionVector begShape2 = begShape.reverse().getSubpart2D(EXT2, begShape.length());
     588       127797 :         const double begSplit = begShape2.nearest_offset_to_point2D(begPoint, false);
     589              : #ifdef DEBUG_SMOOTH_CORNERS
     590              :         if (DEBUGCOND) {
     591              :             std::cout << " begLength=" << begShape2.length2D() << " begSplit=" << begSplit << "\n";
     592              :         }
     593              : #endif
     594       127797 :         if (begSplit > POSITION_EPS && begSplit < begShape2.length2D() - POSITION_EPS) {
     595       118221 :             begShape2 = begShape2.splitAt(begSplit, true).first;
     596              :         } else {
     597              :             return ret;
     598              :         }
     599       118221 :         PositionVector endShape2 = endShape.getSubpart(0, endShape.length() - EXT2);
     600       118221 :         const double endSplit = endShape2.nearest_offset_to_point2D(endPoint, false);
     601              : #ifdef DEBUG_SMOOTH_CORNERS
     602              :         if (DEBUGCOND) {
     603              :             std::cout << " endLength=" << endShape2.length2D() << " endSplit=" << endSplit << "\n";
     604              :         }
     605              : #endif
     606       118221 :         if (endSplit > POSITION_EPS && endSplit < endShape2.length2D() - POSITION_EPS) {
     607       112986 :             endShape2 = endShape2.splitAt(endSplit, true).second;
     608              :         } else {
     609              :             return ret;
     610              :         }
     611              :         // flatten z to junction z level
     612       225972 :         begShape2 = begShape2.interpolateZ(myNode.getPosition().z(), myNode.getPosition().z());
     613       225972 :         endShape2 = endShape2.interpolateZ(myNode.getPosition().z(), myNode.getPosition().z());
     614              : #ifdef DEBUG_SMOOTH_CORNERS
     615              :         if (DEBUGCOND) {
     616              :             std::cout << "getSmoothCorner begPoint=" << begPoint << " endPoint=" << endPoint
     617              :                       << " begShape=" << begShape << " endShape=" << endShape
     618              :                       << " begShape2=" << begShape2 << " endShape2=" << endShape2
     619              :                       << "\n";
     620              :         }
     621              : #endif
     622       112986 :         if (begShape2.size() < 2 || endShape2.size() < 2) {
     623              :             return ret;
     624              :         }
     625       112986 :         const double angle = GeomHelper::angleDiff(begShape2.angleAt2D(-2), endShape2.angleAt2D(0));
     626              :         NBNode* recordError = nullptr;
     627              : #ifdef DEBUG_SMOOTH_CORNERS
     628              :         if (DEBUGCOND) {
     629              :             std::cout << "   angle=" << RAD2DEG(angle) << "\n";
     630              :         }
     631              :         recordError = const_cast<NBNode*>(&myNode);
     632              : #endif
     633              :         // fill highly acute corners
     634              :         //if (fabs(angle) > DEG2RAD(135)) {
     635              :         //    return ret;
     636              :         //}
     637       112986 :         PositionVector curve = myNode.computeSmoothShape(begShape2, endShape2, cornerDetail + 2, false, 25, 25, recordError, NBNode::AVOID_WIDE_LEFT_TURN);
     638              :         //PositionVector curve = myNode.computeSmoothShape(begShape2, endShape2, cornerDetail + 2, false, 25, 25, recordError, 0);
     639       112986 :         const double curvature = curve.length2D() / MAX2(NUMERICAL_EPS, begPoint.distanceTo2D(endPoint));
     640              : #ifdef DEBUG_SMOOTH_CORNERS
     641              :         if (DEBUGCOND) {
     642              :             std::cout << "   curve=" << curve << " curveLength=" << curve.length2D() << " dist=" << begPoint.distanceTo2D(endPoint) << " curvature=" << curvature << "\n";
     643              :         }
     644              : #endif
     645       112986 :         if (curvature > 2 && angle > DEG2RAD(85)) {
     646              :             // simplify dubious inside corner shape
     647              :             return ret;
     648              :         }
     649       111514 :         if (curve.size() > 2) {
     650              :             curve.erase(curve.begin());
     651              :             curve.pop_back();
     652              :             ret = curve;
     653              :         }
     654       127797 :     }
     655              :     return ret;
     656            0 : }
     657              : 
     658              : void
     659        57288 : NBNodeShapeComputer::computeEdgeBoundaries(const EdgeVector& edges,
     660              :         GeomsMap& geomsCCW,
     661              :         GeomsMap& geomsCW) {
     662              :     // compute boundary lines and extend it by EXT m
     663       299525 :     for (NBEdge* const edge : edges) {
     664              :         // store current edge's boundary as current ccw/cw boundary
     665              :         try {
     666       484474 :             geomsCCW[edge] = edge->getCCWBoundaryLine(myNode);
     667            0 :         } catch (InvalidArgument& e) {
     668            0 :             WRITE_WARNING("While computing intersection geometry at junction '" + myNode.getID() + "': " + std::string(e.what()));
     669            0 :             geomsCCW[edge] = edge->getGeometry();
     670            0 :         }
     671              :         try {
     672       484474 :             geomsCW[edge] = edge->getCWBoundaryLine(myNode);
     673            0 :         } catch (InvalidArgument& e) {
     674            0 :             WRITE_WARNING("While computing intersection geometry at junction '" + myNode.getID() + "': " + std::string(e.what()));
     675            0 :             geomsCW[edge] = edge->getGeometry();
     676            0 :         }
     677              :         // ensure the boundary is valid
     678       242237 :         if (geomsCCW[edge].length2D() < NUMERICAL_EPS) {
     679            0 :             geomsCCW[edge] = edge->getGeometry();
     680              :         }
     681       242237 :         if (geomsCW[edge].length2D() < NUMERICAL_EPS) {
     682            0 :             geomsCW[edge] = edge->getGeometry();
     683              :         }
     684              :         // cut off all parts beyond EXT to avoid issues with curved-back roads
     685       484486 :         geomsCCW[edge] = geomsCCW[edge].getSubpart2D(0, MAX2(EXT, edge->getTotalWidth()));
     686       484486 :         geomsCW[edge] = geomsCW[edge].getSubpart2D(0, MAX2(EXT, edge->getTotalWidth()));
     687              :         // extend the boundary by extrapolating it by EXT m towards the junction
     688       242237 :         geomsCCW[edge].extrapolate2D(EXT, true);
     689       242237 :         geomsCW[edge].extrapolate2D(EXT, true);
     690              :         // ensure minimum length by extending it away from the junction
     691       242237 :         geomsCCW[edge].extrapolate(EXT2, false, true);
     692       242237 :         geomsCW[edge].extrapolate(EXT2, false, true);
     693       242237 :         if (geomsCCW[edge].isNAN() || geomsCW[edge].isNAN()) {
     694            0 :             WRITE_WARNING("While computing intersection geometry at junction '" + myNode.getID() + "': found invalid boundary line for edge '" + edge->getID() + "'.")
     695              :         }
     696       242237 :         if (geomsCCW[edge].isNAN()) {
     697            0 :             geomsCCW[edge] = edge->getGeometry();
     698              :         }
     699       242237 :         if (geomsCW[edge].isNAN()) {
     700            0 :             geomsCW[edge] = edge->getGeometry();
     701              :         }
     702              :     }
     703        57288 : }
     704              : 
     705              : void
     706        57288 : NBNodeShapeComputer::joinSameDirectionEdges(const EdgeVector& edges, std::map<NBEdge*, std::set<NBEdge*, ComparatorIdLess> >& same, bool useEndpoints) {
     707              :     // compute same (edges where an intersection doesn't work well
     708              :     // (always check an edge and its cw neighbor)
     709              :     const double angleChangeLookahead = 35; // distance to look ahead for a misleading angle
     710       109784 :     const bool isXodr = OptionsCont::getOptions().exists("opendrive-files") && OptionsCont::getOptions().isSet("opendrive-files");
     711              :     EdgeSet foundOpposite;
     712       299525 :     for (EdgeVector::const_iterator i = edges.begin(); i != edges.end(); i++) {
     713              :         EdgeVector::const_iterator j;
     714       242237 :         if (i == edges.end() - 1) {
     715              :             j = edges.begin();
     716              :         } else {
     717              :             j = i + 1;
     718              :         }
     719           16 :         if (useEndpoints
     720           16 :                 && !(*i)->hasDefaultGeometryEndpointAtNode(&myNode)
     721       242253 :                 && !(*j)->hasDefaultGeometryEndpointAtNode(&myNode)) {
     722           16 :             continue;
     723              :         }
     724       242221 :         const bool incoming = (*i)->getToNode() == &myNode;
     725       242221 :         const bool incoming2 = (*j)->getToNode() == &myNode;
     726              :         const bool differentDirs = (incoming != incoming2);
     727       278124 :         const bool sameGeom = (*i)->getGeometry() == (differentDirs ? (*j)->getGeometry().reverse() : (*j)->getGeometry());
     728       242221 :         const PositionVector g1 = incoming ? (*i)->getCCWBoundaryLine(myNode) : (*i)->getCWBoundaryLine(myNode);
     729       242221 :         const PositionVector g2 = incoming ? (*j)->getCCWBoundaryLine(myNode) : (*j)->getCWBoundaryLine(myNode);
     730       242221 :         const double angle1further = (g1.size() > 2 && g1[0].distanceTo2D(g1[1]) < angleChangeLookahead ?
     731       242221 :                                       g1.angleAt2D(1) : g1.angleAt2D(0));
     732       242221 :         const double angle2further = (g2.size() > 2 && g2[0].distanceTo2D(g2[1]) < angleChangeLookahead ?
     733       242221 :                                       g2.angleAt2D(1) : g2.angleAt2D(0));
     734       242221 :         const double angleDiff = GeomHelper::angleDiff(g1.angleAt2D(0), g2.angleAt2D(0));
     735       242221 :         const double angleDiffFurther = GeomHelper::angleDiff(angle1further, angle2further);
     736       242221 :         const bool ambiguousGeometry = ((angleDiff > 0 && angleDiffFurther < 0) || (angleDiff < 0 && angleDiffFurther > 0));
     737              :         //if (ambiguousGeometry) {
     738              :         //    @todo: this warning would be helpful in many cases. However, if angle and angleFurther jump between 179 and -179 it is misleading
     739              :         //    WRITE_WARNINGF(TL("Ambiguous angles at junction '%' for edges '%' and '%'."), myNode.getID(), (*i)->getID(), (*j)->getID());
     740              :         //}
     741              : #ifdef DEBUG_NODE_SHAPE
     742              :         if (DEBUGCOND) {
     743              :             std::cout << "   checkSameDirection " << (*i)->getID() << " " << (*j)->getID()
     744              :                       << " sameGeom=" << sameGeom
     745              :                       << " diffDirs=" << differentDirs
     746              :                       << " isOpposite=" << (differentDirs && foundOpposite.count(*i) == 0)
     747              :                       << " angleDiff=" << angleDiff
     748              :                       << " ambiguousGeometry=" << ambiguousGeometry
     749              :                       << " badInsersection=" << badIntersection(*i, *j, EXT)
     750              :                       << "\n";
     751              : 
     752              :         }
     753              : #endif
     754       242221 :         if (sameGeom || fabs(angleDiff) < DEG2RAD(20)) {
     755        79658 :             const bool isOpposite = differentDirs && foundOpposite.count(*i) == 0;
     756              :             if (isOpposite) {
     757              :                 foundOpposite.insert(*i);
     758              :                 foundOpposite.insert(*j);
     759              :             }
     760        79658 :             if (isOpposite || ambiguousGeometry || (!isXodr && badIntersection(*i, *j, EXT))) {
     761              :                 // maintain equivalence relation for all members of the equivalence class
     762        82241 :                 for (std::set<NBEdge*>::iterator k = same[*i].begin(); k != same[*i].end(); ++k) {
     763         4563 :                     if (*j != *k) {
     764         4505 :                         same[*k].insert(*j);
     765         4505 :                         same[*j].insert(*k);
     766              :                     }
     767              :                 }
     768        82434 :                 for (std::set<NBEdge*>::iterator k = same[*j].begin(); k != same[*j].end(); ++k) {
     769         4756 :                     if (*i != *k) {
     770         4698 :                         same[*k].insert(*i);
     771         4698 :                         same[*i].insert(*k);
     772              :                     }
     773              :                 }
     774        77678 :                 same[*i].insert(*j);
     775        77678 :                 same[*j].insert(*i);
     776              : #ifdef DEBUG_NODE_SHAPE
     777              :                 if (DEBUGCOND) {
     778              :                     std::cout << "   joinedSameDirectionEdges " << (*i)->getID() << "   " << (*j)->getID() << " isOpposite=" << isOpposite << " ambiguousGeometry=" << ambiguousGeometry << "\n";
     779              :                 }
     780              : #endif
     781              :             }
     782              :         }
     783       242221 :     }
     784        57288 : }
     785              : 
     786              : 
     787              : bool
     788         4795 : NBNodeShapeComputer::badIntersection(const NBEdge* e1, const NBEdge* e2, double distance) {
     789              :     // check whether the two edges are on top of each other. In that case they should be joined
     790              :     // also, if they never touch along their common length
     791         4795 :     const double commonLength = MIN3(distance, e1->getGeometry().length(), e2->getGeometry().length());
     792              :     PositionVector geom1 = e1->getGeometry();
     793              :     PositionVector geom2 = e2->getGeometry();
     794              :     // shift to make geom the centerline of the edge regardless of spreadtype
     795         4795 :     if (e1->getLaneSpreadFunction() == LaneSpreadFunction::RIGHT) {
     796         1012 :         geom1.move2side(e1->getTotalWidth() / 2);
     797              :     }
     798         4795 :     if (e2->getLaneSpreadFunction() == LaneSpreadFunction::RIGHT) {
     799          967 :         geom2.move2side(e2->getTotalWidth() / 2);
     800              :     }
     801              :     // always let geometry start at myNode
     802         4795 :     if (e1->getToNode() == &myNode) {
     803         3748 :         geom1 = geom1.reverse();
     804              :     }
     805         4795 :     if (e2->getToNode() == &myNode) {
     806         5146 :         geom2 = geom2.reverse();
     807              :     }
     808         9590 :     geom1 = geom1.getSubpart2D(0, commonLength);
     809         9590 :     geom2 = geom2.getSubpart2D(0, commonLength);
     810              :     double endAngleDiff = 0;
     811         4795 :     if (geom1.size() >= 2 && geom2.size() >= 2) {
     812         4795 :         endAngleDiff = fabs(RAD2DEG(GeomHelper::angleDiff(
     813              :                                         geom1.angleAt2D((int)geom1.size() - 2),
     814              :                                         geom2.angleAt2D((int)geom2.size() - 2))));
     815              :     }
     816         4795 :     const double minDistanceThreshold = (e1->getTotalWidth() + e2->getTotalWidth()) / 2 + POSITION_EPS;
     817         4795 :     std::vector<double> distances = geom1.distances(geom2, true);
     818         4795 :     std::vector<double> distances2 = geom1.distances(geom2);
     819              :     const double minDist = VectorHelper<double>::minValue(distances2);
     820              :     const double maxDist = VectorHelper<double>::maxValue(distances);
     821         4795 :     const bool curvingTowards = geom1[0].distanceTo2D(geom2[0]) > minDistanceThreshold && minDist < minDistanceThreshold;
     822         4795 :     const bool onTop = (maxDist - POSITION_EPS < minDistanceThreshold) && endAngleDiff < 30;
     823         4795 :     const bool bothDefault = e1->hasDefaultGeometryEndpointAtNode(&myNode) && e2->hasDefaultGeometryEndpointAtNode(&myNode);
     824         4795 :     const bool neverTouch = minDist > minDistanceThreshold * 2 && !bothDefault;
     825         4795 :     geom1.extrapolate2D(EXT);
     826         4795 :     geom2.extrapolate2D(EXT);
     827         4795 :     Position intersect = geom1.intersectionPosition2D(geom2);
     828         3540 :     const bool intersects = intersect != Position::INVALID && geom1.distance2D(intersect) < POSITION_EPS;
     829              : #ifdef DEBUG_NODE_SHAPE
     830              :     if (DEBUGCOND) {
     831              :         std::cout << "    badIntersect: onTop=" << onTop << " curveTo=" << curvingTowards << " intersects=" << intersects
     832              :                   << " endAngleDiff=" << endAngleDiff
     833              :                   << " geom1=" << geom1 << " geom2=" << geom2
     834              :                   << " distances=" << toString(distances) << " minDist=" << minDist << " maxDist=" << maxDist << " thresh=" << minDistanceThreshold
     835              :                   << " neverTouch=" << neverTouch
     836              :                   << " intersectPos=" << intersect
     837              :                   << "\n";
     838              :     }
     839              : #endif
     840         9590 :     return onTop || curvingTowards || !intersects || neverTouch;
     841         4795 : }
     842              : 
     843              : 
     844              : EdgeVector
     845        57288 : NBNodeShapeComputer::computeUniqueDirectionList(
     846              :     const EdgeVector& all,
     847              :     std::map<NBEdge*, std::set<NBEdge*, ComparatorIdLess> >& same,
     848              :     GeomsMap& geomsCCW,
     849              :     GeomsMap& geomsCW) {
     850              :     // store relationships
     851        57288 :     EdgeVector newAll = all;
     852       299525 :     for (NBEdge* e1 : all) {
     853              :         // determine which of the edges marks the outer boundary
     854       242237 :         auto e2NewAll = std::find(newAll.begin(), newAll.end(), e1);
     855              : #ifdef DEBUG_NODE_SHAPE
     856              :         if (DEBUGCOND) std::cout << "computeUniqueDirectionList e1=" << e1->getID()
     857              :                                      << " deleted=" << (e2NewAll == newAll.end())
     858              :                                      << " same=" << joinNamedToStringSorting(same[e1], ',') << "\n";
     859              : #endif
     860       242237 :         if (e2NewAll == newAll.end()) {
     861        77620 :             continue;
     862              :         }
     863       164617 :         auto e1It = std::find(all.begin(), all.end(), e1);
     864              :         auto bestCCW = e1It;
     865              :         auto bestCW = e1It;
     866              :         bool changed = true;
     867       405669 :         while (changed) {
     868              :             changed = false;
     869       400990 :             for (NBEdge* e2 : same[e1]) {
     870              : #ifdef DEBUG_NODE_SHAPE
     871              :                 if (DEBUGCOND) {
     872              :                     std::cout << "  e2=" << e2->getID() << "\n";
     873              :                 }
     874              : #endif
     875       159938 :                 auto e2It = std::find(all.begin(), all.end(), e2);
     876       159938 :                 if (e2It + 1 == bestCCW || (e2It == (all.end() - 1) && bestCCW == all.begin())) {
     877              :                     bestCCW = e2It;
     878              :                     changed = true;
     879              : #ifdef DEBUG_NODE_SHAPE
     880              :                     if (DEBUGCOND) {
     881              :                         std::cout << "    bestCCW=" << e2->getID() << "\n";
     882              :                     }
     883              : #endif
     884       158911 :                 } else if (bestCW + 1 == e2It || (bestCW == (all.end() - 1) && e2It == all.begin())) {
     885              :                     bestCW = e2It;
     886              :                     changed = true;
     887              : #ifdef DEBUG_NODE_SHAPE
     888              :                     if (DEBUGCOND) {
     889              :                         std::cout << "    bestCW=" << e2->getID() << "\n";
     890              :                     }
     891              : #endif
     892              :                 }
     893              :             }
     894              :         }
     895       164617 :         if (bestCW != e1It) {
     896        74703 :             geomsCW[e1] = geomsCW[*bestCW];
     897        74703 :             computeSameEnd(geomsCW[e1], geomsCCW[e1]);
     898              :         }
     899       164617 :         if (bestCCW != e1It) {
     900          849 :             geomsCCW[e1] = geomsCCW[*bestCCW];
     901          849 :             computeSameEnd(geomsCW[e1], geomsCCW[e1]);
     902              :         }
     903              :         // clean up
     904       242237 :         for (NBEdge* e2 : same[e1]) {
     905        77620 :             auto e2NewAllIt = std::find(newAll.begin(), newAll.end(), e2);
     906        77620 :             if (e2NewAllIt != newAll.end()) {
     907              :                 newAll.erase(e2NewAllIt);
     908              :             }
     909              :         }
     910              :     }
     911              : #ifdef DEBUG_NODE_SHAPE
     912              :     if (DEBUGCOND) {
     913              :         std::cout << "  newAll:\n";
     914              :         for (NBEdge* e : newAll) {
     915              :             std::cout << "    " << e->getID() << " geomCCW=" << geomsCCW[e] << " geomsCW=" << geomsCW[e] << "\n";
     916              :         }
     917              :     }
     918              : #endif
     919        57288 :     return newAll;
     920            0 : }
     921              : 
     922              : 
     923              : void
     924       164534 : NBNodeShapeComputer::initNeighbors(const EdgeVector& edges, const EdgeVector::const_iterator& current,
     925              :                                    GeomsMap& geomsCW,
     926              :                                    GeomsMap& geomsCCW,
     927              :                                    EdgeVector::const_iterator& cwi,
     928              :                                    EdgeVector::const_iterator& ccwi,
     929              :                                    double& cad,
     930              :                                    double& ccad) {
     931              :     const double twoPI = (double)(2 * M_PI);
     932       164534 :     cwi = current;
     933              :     cwi++;
     934       164534 :     if (cwi == edges.end()) {
     935        57205 :         std::advance(cwi, -((int)edges.size())); // set to edges.begin();
     936              :     }
     937       164534 :     ccwi = current;
     938       164534 :     if (ccwi == edges.begin()) {
     939        57205 :         std::advance(ccwi, edges.size() - 1); // set to edges.end() - 1;
     940              :     } else {
     941              :         ccwi--;
     942              :     }
     943              : 
     944       164534 :     const double angleCurCCW = geomsCCW[*current].angleAt2D(0);
     945       164534 :     const double angleCurCW = geomsCW[*current].angleAt2D(0);
     946       164534 :     const double angleCCW = geomsCW[*ccwi].angleAt2D(0);
     947       164534 :     const double angleCW = geomsCCW[*cwi].angleAt2D(0);
     948       164534 :     ccad = angleCCW - angleCurCCW;
     949       223294 :     while (ccad < 0.) {
     950        58760 :         ccad += twoPI;
     951              :     }
     952       164534 :     cad = angleCurCW - angleCW;
     953       223294 :     while (cad < 0.) {
     954        58760 :         cad += twoPI;
     955              :     }
     956       164534 : }
     957              : 
     958              : 
     959              : 
     960              : const PositionVector
     961        20994 : NBNodeShapeComputer::computeNodeShapeSmall() {
     962              : #ifdef DEBUG_NODE_SHAPE
     963              :     if (DEBUGCOND) {
     964              :         std::cout << "computeNodeShapeSmall node=" << myNode.getID() << "\n";
     965              :     }
     966              : #endif
     967        20994 :     PositionVector ret;
     968        53133 :     for (NBEdge* e : myNode.getEdges()) {
     969              :         // compute crossing with normal
     970        32139 :         PositionVector edgebound1 = e->getCCWBoundaryLine(myNode).getSubpartByIndex(0, 2);
     971        32139 :         PositionVector edgebound2 = e->getCWBoundaryLine(myNode).getSubpartByIndex(0, 2);
     972        32139 :         Position delta = edgebound1[1] - edgebound1[0];
     973              :         delta.set(-delta.y(), delta.x()); // rotate 90 degrees
     974        32139 :         PositionVector cross(myNode.getPosition(), myNode.getPosition() + delta);
     975        32139 :         cross.extrapolate2D(500);
     976        32139 :         edgebound1.extrapolate2D(500);
     977        32139 :         edgebound2.extrapolate2D(500);
     978        32139 :         if (cross.intersects(edgebound1)) {
     979        32139 :             Position np = cross.intersectionPosition2D(edgebound1);
     980        32139 :             np.set(np.x(), np.y(), myNode.getPosition().z());
     981        32139 :             ret.push_back_noDoublePos(np);
     982              :         }
     983        32139 :         if (cross.intersects(edgebound2)) {
     984        32135 :             Position np = cross.intersectionPosition2D(edgebound2);
     985        32135 :             np.set(np.x(), np.y(), myNode.getPosition().z());
     986        32135 :             ret.push_back_noDoublePos(np);
     987              :         }
     988        32139 :         e->resetNodeBorder(&myNode);
     989        32139 :     }
     990        20994 :     return ret;
     991            0 : }
     992              : 
     993              : 
     994              : double
     995        57288 : NBNodeShapeComputer::getDefaultRadius(const OptionsCont& oc) {
     996              :     // look for incoming/outgoing edge pairs that do not go straight and allow wide vehicles
     997              :     // (connection information is not available yet)
     998              :     // @TODO compute the radius for each pair of neighboring edge intersections in computeNodeShapeDefault rather than use the maximum
     999        57288 :     const double radius = oc.getFloat("default.junctions.radius");
    1000        57288 :     const double smallRadius = oc.getFloat("junctions.small-radius");
    1001              :     double maxRightAngle = 0; // rad
    1002              :     double extraWidthRight = 0; // m
    1003              :     double maxLeftAngle = 0; // rad
    1004              :     double extraWidthLeft = 0; // m
    1005              :     int laneDelta = 0;
    1006              :     int totalWideLanesIn = 0;
    1007       177953 :     for (NBEdge* in : myNode.getIncomingEdges()) {
    1008              :         int wideLanesIn = 0;
    1009       277869 :         for (int i = 0; i < in->getNumLanes(); i++) {
    1010       157204 :             if ((in->getPermissions(i) & SVC_LARGE_TURN) != 0) {
    1011       103078 :                 wideLanesIn++;
    1012              :             }
    1013              :         }
    1014       120665 :         totalWideLanesIn += wideLanesIn;
    1015       426829 :         for (NBEdge* out : myNode.getOutgoingEdges()) {
    1016       306164 :             if ((in->getPermissions() & out->getPermissions() & SVC_LARGE_TURN) != 0) {
    1017       194304 :                 if (myNode.getDirection(in, out) == LinkDirection::TURN) {
    1018        53102 :                     continue;
    1019              :                 };
    1020       141202 :                 const double angle = GeomHelper::angleDiff(
    1021       141202 :                                          in->getGeometry().angleAt2D(-2),
    1022              :                                          out->getGeometry().angleAt2D(0));
    1023       141202 :                 if (angle < 0) {
    1024        65781 :                     if (maxRightAngle < -angle) {
    1025              :                         maxRightAngle = -angle;
    1026        40664 :                         extraWidthRight = MAX2(getExtraWidth(in, SVC_LARGE_TURN), getExtraWidth(out, SVC_LARGE_TURN));
    1027              :                     }
    1028              :                 } else {
    1029        75421 :                     if (maxLeftAngle < angle) {
    1030              :                         maxLeftAngle = angle;
    1031              :                         // all edges clockwise between in and out count as extra width
    1032              :                         extraWidthLeft = 0;
    1033        35960 :                         EdgeVector::const_iterator pIn = std::find(myNode.getEdges().begin(), myNode.getEdges().end(), in);
    1034        35960 :                         NBContHelper::nextCW(myNode.getEdges(), pIn);
    1035        93314 :                         while (*pIn != out) {
    1036        57354 :                             extraWidthLeft += (*pIn)->getTotalWidth();
    1037              : #ifdef DEBUG_RADIUS
    1038              :                             if (DEBUGCOND) {
    1039              :                                 std::cout << "   in=" << in->getID() << " out=" << out->getID() << " extra=" << (*pIn)->getID() << " extraWidthLeft=" << extraWidthLeft << "\n";
    1040              :                             }
    1041              : #endif
    1042        57354 :                             NBContHelper::nextCW(myNode.getEdges(), pIn);
    1043              :                         }
    1044              :                     }
    1045              :                 }
    1046              :                 int wideLanesOut = 0;
    1047       341113 :                 for (int i = 0; i < out->getNumLanes(); i++) {
    1048       199911 :                     if ((out->getPermissions(i) & SVC_LARGE_TURN) != 0) {
    1049       186483 :                         wideLanesOut++;
    1050              :                     }
    1051              :                 }
    1052              : #ifdef DEBUG_RADIUS
    1053              :                 if (DEBUGCOND) {
    1054              :                     std::cout << "   in=" << in->getID() << " out=" << out->getID() << " wideLanesIn=" << wideLanesIn << " wideLanesOut=" << wideLanesOut << "\n";
    1055              :                 }
    1056              : #endif
    1057       141202 :                 laneDelta = MAX2(laneDelta, abs(wideLanesOut - wideLanesIn));
    1058              :             }
    1059              :         }
    1060              :     }
    1061              :     // special case: on/off-ramp
    1062        57288 :     if (myNode.getOutgoingEdges().size() == 1 || myNode.getIncomingEdges().size() == 1) {
    1063              :         int totalWideLanesOut = 0;
    1064        53870 :         for (NBEdge* out : myNode.getOutgoingEdges()) {
    1065        74349 :             for (int i = 0; i < out->getNumLanes(); i++) {
    1066        43645 :                 if ((out->getPermissions(i) & SVC_LARGE_TURN) != 0) {
    1067        25836 :                     totalWideLanesOut++;
    1068              :                 }
    1069              :             }
    1070              :         }
    1071        23166 :         if (totalWideLanesIn == totalWideLanesOut) {
    1072              :             // use total laneDelta instead of individual edge lane delta
    1073              :             laneDelta = 0;
    1074              :         }
    1075              :     }
    1076              :     // changing the number of wide-vehicle lanes on a straight segment requires a larger junction to allow for smooth driving
    1077              :     // otherwise we can reduce the radius according to the angle
    1078              :     double result = radius;
    1079              :     // left turns are assumed to cross additional edges and thus du not determine the required radius in most cases
    1080              :     double maxTurnAngle = maxRightAngle;
    1081              :     double extraWidth = extraWidthRight;
    1082        57288 :     if (maxRightAngle < DEG2RAD(5)) {
    1083              :         maxTurnAngle = maxLeftAngle;
    1084              :         extraWidth = extraWidthLeft;
    1085              :     }
    1086        57288 :     const double minRadius = maxTurnAngle >= DEG2RAD(30) ? MIN2(smallRadius, radius) : smallRadius;
    1087        57288 :     if (laneDelta == 0 || maxTurnAngle >= DEG2RAD(30) || myNode.isConstantWidthTransition()) {
    1088              :         // subtract radius gained from extra lanes
    1089              :         // do not increase radius for turns that are sharper than a right angle
    1090        54028 :         result = radius * tan(0.5 * MIN2(0.5 * M_PI, maxTurnAngle)) - extraWidth;
    1091              :     }
    1092              :     result = MAX2(minRadius, result);
    1093              : #ifdef DEBUG_RADIUS
    1094              :     if (DEBUGCOND) {
    1095              :         std::cout << "getDefaultRadius n=" << myNode.getID()
    1096              :                   << " r=" << radius << " sr=" << smallRadius
    1097              :                   << " mr=" << minRadius
    1098              :                   << " laneDelta=" << laneDelta
    1099              :                   << " rightA=" << RAD2DEG(maxRightAngle)
    1100              :                   << " leftA=" << RAD2DEG(maxLeftAngle)
    1101              :                   << " maxA=" << RAD2DEG(maxTurnAngle)
    1102              :                   << " extraWidth=" << extraWidth
    1103              :                   << " result=" << result << "\n";
    1104              :     }
    1105              : #endif
    1106        57288 :     return result;
    1107              : }
    1108              : 
    1109              : 
    1110              : bool
    1111       136900 : NBNodeShapeComputer::isDivided(const NBEdge* e, std::set<NBEdge*, ComparatorIdLess> same, const PositionVector& ccw, const PositionVector& cw) const {
    1112       136900 :     if (same.size() < 2) {
    1113              :         return false;
    1114              :     }
    1115              :     std::set<Position> endPoints;
    1116         1778 :     endPoints.insert(e->getEndpointAtNode(&myNode));
    1117         3547 :     for (NBEdge* s : same) {
    1118         5316 :         endPoints.insert(s->getEndpointAtNode(&myNode));
    1119              :     }
    1120          889 :     if (endPoints.size() > 1) {
    1121          743 :         std::vector<double> distances = ccw.distances(cw, true);
    1122          743 :         double width = e->getTotalWidth();
    1123         3047 :         for (const NBEdge* e2 : same) {
    1124         2304 :             width += e2->getTotalWidth();
    1125              :         }
    1126              :         const double maxDist = VectorHelper<double>::maxValue(distances);
    1127          743 :         const double maxDivider = maxDist - width;
    1128          743 :         return maxDivider >= 5;
    1129          743 :     }
    1130              :     return false;
    1131              : }
    1132              : 
    1133              : 
    1134              : double
    1135        81328 : NBNodeShapeComputer::getExtraWidth(const NBEdge* e, SVCPermissions exclude) {
    1136              :     double result = 0;
    1137              :     int lane = 0;
    1138        81328 :     while (lane < e->getNumLanes() && e->getPermissions(lane) == 0) {
    1139              :         // ignore forbidden lanes out the outside
    1140            0 :         lane++;
    1141              :     }
    1142        90249 :     while (lane < e->getNumLanes() && (e->getPermissions(lane) & exclude) == 0) {
    1143         8921 :         result += e->getLaneWidth(lane);
    1144         8921 :         lane++;
    1145              :     }
    1146        81328 :     return result;
    1147              : }
    1148              : 
    1149              : 
    1150              : double
    1151          313 : NBNodeShapeComputer::divisionWidth(const NBEdge* e, std::set<NBEdge*, ComparatorIdLess> same, const Position& p, const Position& p2) {
    1152              :     double result = p.distanceTo2D(p2);
    1153          313 :     result -= e->getTotalWidth();
    1154         1466 :     for (NBEdge* e2 : same) {
    1155         1153 :         result -= e2->getTotalWidth();
    1156              :     }
    1157          313 :     return MAX2(0.0, result);
    1158              : }
    1159              : 
    1160              : /****************************************************************************/
        

Generated by: LCOV version 2.0-1