Skip to content

Commit

Permalink
use double precision fp math for mercartor projection and point-line …
Browse files Browse the repository at this point in the history
…projections (for now), fixes #1191
  • Loading branch information
DennisOSRM committed Oct 9, 2014
1 parent 2b9e253 commit 440244e
Show file tree
Hide file tree
Showing 2 changed files with 17 additions and 17 deletions.
26 changes: 13 additions & 13 deletions DataStructures/Coordinate.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -258,16 +258,16 @@ float FixedPointCoordinate::ComputePerpendicularDistance(const FixedPointCoordin
BOOST_ASSERT(query_location.isValid());

// initialize values
const float x = lat2y(query_location.lat / COORDINATE_PRECISION);
const float y = query_location.lon / COORDINATE_PRECISION;
const float a = lat2y(segment_source.lat / COORDINATE_PRECISION);
const float b = segment_source.lon / COORDINATE_PRECISION;
const float c = lat2y(segment_target.lat / COORDINATE_PRECISION);
const float d = segment_target.lon / COORDINATE_PRECISION;
float p, q /*,mX*/, nY;
if (std::abs(a - c) > std::numeric_limits<float>::epsilon())
const double x = lat2y(query_location.lat / COORDINATE_PRECISION);
const double y = query_location.lon / COORDINATE_PRECISION;
const double a = lat2y(segment_source.lat / COORDINATE_PRECISION);
const double b = segment_source.lon / COORDINATE_PRECISION;
const double c = lat2y(segment_target.lat / COORDINATE_PRECISION);
const double d = segment_target.lon / COORDINATE_PRECISION;
double p, q /*,mX*/, nY;
if (std::abs(a - c) > std::numeric_limits<double>::epsilon())
{
const float m = (d - b) / (c - a); // slope
const double m = (d - b) / (c - a); // slope
// Projection of (x,y) on line joining (a,b) and (c,d)
p = ((x + (m * y)) + (m * m * a - m * b)) / (1.f + m * m);
q = b + m * (p - a);
Expand All @@ -293,11 +293,11 @@ float FixedPointCoordinate::ComputePerpendicularDistance(const FixedPointCoordin
{
ratio = (segment_target == query_location ? 1.f : 0.f);
}
else if (std::abs(ratio) <= std::numeric_limits<float>::epsilon())
else if (std::abs(ratio) <= std::numeric_limits<double>::epsilon())
{
ratio = 0.;
ratio = 0.f;
}
else if (std::abs(ratio - 1.f) <= std::numeric_limits<float>::epsilon())
else if (std::abs(ratio - 1.f) <= std::numeric_limits<double>::epsilon())
{
ratio = 1.f;
}
Expand All @@ -308,7 +308,7 @@ float FixedPointCoordinate::ComputePerpendicularDistance(const FixedPointCoordin
{
nearest_location = segment_source;
}
else if (ratio >= 1.)
else if (ratio >= 1.f)
{
nearest_location = segment_target;
}
Expand Down
8 changes: 4 additions & 4 deletions Util/MercatorUtil.h
Original file line number Diff line number Diff line change
Expand Up @@ -30,14 +30,14 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#include <cmath>

inline float y2lat(const float a)
inline double y2lat(const double a)
{
return 180.f * static_cast<float>(M_1_PI) * (2.f * std::atan(std::exp(a * static_cast<float>(M_PI) / 180.f)) - static_cast<float>(M_PI_2));
return 180. * M_1_PI * (2. * std::atan(std::exp(a * M_PI / 180.)) - M_PI_2);
}

inline float lat2y(const float a)
inline double lat2y(const double a)
{
return 180.f * static_cast<float>(M_1_PI) * std::log(std::tan(static_cast<float>(M_PI_4) + a * (static_cast<float>(M_PI) / 180.f) / 2.f));
return 180. * M_1_PI * std::log(std::tan(M_PI_4 + a * (M_PI / 180.) / 2.));
}

#endif // MERCATOR_UTIL_H

0 comments on commit 440244e

Please sign in to comment.