correct turn angles (and instructions) at high latitudes #596

Closed
wants to merge 5 commits into
from

Conversation

Projects
None yet
2 participants
Contributor

emiltin commented Feb 23, 2013

solves the problem of distorted angles, and thus incorrect turn instructions, at high latitudes.

keeps the original atan() formula, but scales lat deltas by the longitudinal shortening at the C lat. only adds a single cos() call.

Contributor

emiltin commented Feb 23, 2013

the longitude shortening here in copenhagen is about 0.57, so the distortion of angles is significant without correction.

Contributor

DennisOSRM commented Feb 23, 2013

It's just missing Mercator projection, isn't it?

Contributor

emiltin commented Feb 23, 2013

right. needs to scale the lat

Contractor/EdgeBasedGraphFactory.cpp
+ const double v2x = B.lon - C.lon;
+ const double v2y = B.lat - C.lat;
+ const double latC = (C.lat/100000.)*M_PI/180.;
+ const double scale = cos(latC); //scale by length of longitude at latitude C
@DennisOSRM

DennisOSRM Feb 23, 2013

Contributor

Not sure how geometrically sound your solution is.

Why not use Mercartor projection to fix this?

template<class CoordinateT>
double EdgeBasedGraphFactory::GetAngleBetweenTwoEdges(const CoordinateT& A, const CoordinateT& C, const CoordinateT& B) const {
    const double v1x = (A.lon - C.lon)/100000.;
    const double v1y = lat2y((A.lat - C.lat)/100000.);
    const double v2x = (B.lon - C.lon)/100000.;
    const double v2y = lat2y((B.lat - C.lat)/100000.);

    double angle = (atan2(v2y,v2x) - atan2(v1y,v1x) )*180/M_PI;
    while(angle < 0)
        angle += 360;
    return angle;
}

The lat2y() function is in DataStructures/MercatorUtil.h and one of benefits is that Mercator Projection preserves angles.

Contributor

emiltin commented Feb 23, 2013

ok that should be more precise then.

Contributor

emiltin commented Feb 24, 2013

ok updated it to use mercator projection as suggested, however you need to use lat2y on each lat, not on the difference.

the diff of two mercator conversions can be combined since log(a)-log(b)=log(a/b). the more straightforward approach would be:

    const double latA = lat2yrad(int2rad(A.lat));
    const double latB = lat2yrad(int2rad(B.lat));
    const double latC = lat2yrad(int2rad(C.lat));

    const double lonA = int2rad(A.lon);
    const double lonB = int2rad(B.lon);
    const double lonC = int2rad(C.lon);

    const double v1x = lonA - lonC;
    const double v2x = lonB - lonC;
    const double v1y = latA - latC;
    const double v2y = latB - latC;

    const double a2 = rad2deg(atan2(v2y,v2x));
    const double a1 = rad2deg(atan2(v1y,v1x));
    double angle = a2-a1;   

wouldn't it make sense to move all the geometry utilities to one place? at the moment, geometry utilities are found in MercatorUtil.h, DescriptionFactory.cpp and Coordinate.h. i also needed a few more, for now i placed them in a new file Geometry.h.

btw, why do you use a while loop? isn't it enough to add 360 if the angle is below zero?

Contributor

emiltin commented Feb 24, 2013

to make testing things like this easier, cuke test could accept coords in x,y format, which would then be transformed to lat,lon using inverse mercator before being written to the osm file.

Contractor/EdgeBasedGraphFactory.cpp
+ const double v1x = lonA - lonC;
+ const double v2x = lonB - lonC;
+
+ // mercator projection: y=log(tan( pi/4 + lat/2 ))
@DennisOSRM

DennisOSRM Feb 24, 2013

Contributor

There is already functionality for that implemented.

@emiltin

emiltin Feb 24, 2013

Contributor

yes, but it works with degrees. using it would mean you first convert to deg, and the function internally converts to rad, computes and converts back to deg. then finally you convert that back to rad.

plus you can get replace a log() call with a division by combining the two angles. might be faster, not sure.

Contractor/EdgeBasedGraphFactory.cpp
+ const double v2y = log(tan(M_PI/4 + latB/2)/tan(M_PI/4 + latC/2));
+
+ const double a2 = atan2(v2y,v2x);
+ const double a1 = atan2(v1y,v1x);
@DennisOSRM

DennisOSRM Feb 24, 2013

Contributor

check indentation

Contributor

DennisOSRM commented Feb 24, 2013

The while loop is a precaution.

Contributor

DennisOSRM commented Feb 28, 2013

Merged the test and fixed the issue by re-using the existing mercartor projection code.

@DennisOSRM DennisOSRM closed this Feb 28, 2013

Zhdanovich pushed a commit to Zhdanovich/Project-OSRM that referenced this pull request Jul 20, 2014

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment