From 6fd50229ff285160f00cfd6d223c8dced72f4b8c Mon Sep 17 00:00:00 2001 From: Andrew Kirmse Date: Tue, 12 Mar 2024 19:09:12 -0700 Subject: [PATCH 1/2] Fix rounding in coordinates in filenames. --- code/degree_coordinate_system.cpp | 4 ++++ code/flt_loader.cpp | 6 +++++- code/prominence_task.cpp | 18 +++++++++++++----- code/tile_loading_policy.cpp | 9 +++++++++ 4 files changed, 31 insertions(+), 6 deletions(-) diff --git a/code/degree_coordinate_system.cpp b/code/degree_coordinate_system.cpp index 27967e2..4557b45 100644 --- a/code/degree_coordinate_system.cpp +++ b/code/degree_coordinate_system.cpp @@ -82,6 +82,10 @@ Offsets DegreeCoordinateSystem::offsetsTo(const CoordinateSystem &that) { int dx = static_cast(std::round((mMinLongitude - other->mMinLongitude) * mSamplesPerDegreeLongitude)); int dy = static_cast(std::round((other->mMaxLatitude - mMaxLatitude) * mSamplesPerDegreeLatitude)); + // XXX + VLOG(1) << "Ours " << toString(); + VLOG(1) << "Theirs " << that.toString(); + VLOG(1) << "dx " << dx << " dy " << dy; return Offsets(dx, dy); } diff --git a/code/flt_loader.cpp b/code/flt_loader.cpp index a360a02..86c7798 100644 --- a/code/flt_loader.cpp +++ b/code/flt_loader.cpp @@ -62,9 +62,13 @@ Tile *FltLoader::loadFromFltFile(const string &directory, float minLat, float mi filename = directory + "/" + filename; } + // XXX + VLOG(1) << "Loading FLT file " << filename; + FILE *infile = fopen(filename.c_str(), "rb"); if (infile == nullptr) { - VLOG(3) << "Failed to open file " << filename; +// XXX VLOG(3) << "Failed to open file " << filename; + VLOG(1) << "Failed to open file " << filename; return nullptr; } diff --git a/code/prominence_task.cpp b/code/prominence_task.cpp index fbae74e..8e2fdd3 100644 --- a/code/prominence_task.cpp +++ b/code/prominence_task.cpp @@ -55,7 +55,7 @@ bool ProminenceTask::run(float lat, float lng, const CoordinateSystem &coordinat VLOG(2) << "Couldn't load tile for " << lat << " " << lng; return false; } - + // Flip tile upside down if we're computing anti-prominence if (mOptions.antiprominence) { tile->flipElevations(); @@ -115,11 +115,19 @@ bool ProminenceTask::writeStringToOutputFile(const string &filename, const strin string ProminenceTask::getFilenamePrefix() const { char filename[PATH_MAX]; - int latHundredths = fractionalDegree(mCurrentLatitude); - int lngHundredths = fractionalDegree(mCurrentLongitude); + + // XXX + const float epsilon = 0.001f; + float lat = mCurrentLatitude; + float lng = mCurrentLongitude; + lat += (lat >= 0) ? epsilon : -epsilon; + lng += (lng >= 0) ? epsilon : -epsilon; + + int latHundredths = fractionalDegree(lat); + int lngHundredths = fractionalDegree(lng); snprintf(filename, sizeof(filename), "prominence-%02dx%02d-%03dx%02d", - static_cast(mCurrentLatitude), latHundredths, - static_cast(mCurrentLongitude), lngHundredths); + static_cast(lat), latHundredths, + static_cast(lng), lngHundredths); return mOptions.outputDir + "/" + filename; } diff --git a/code/tile_loading_policy.cpp b/code/tile_loading_policy.cpp index c4c792c..cea3e4f 100644 --- a/code/tile_loading_policy.cpp +++ b/code/tile_loading_policy.cpp @@ -103,6 +103,15 @@ Tile *BasicTileLoadingPolicy::loadTile(float minLat, float minLng) const { Tile *BasicTileLoadingPolicy::loadInternal(float minLat, float minLng) const { TileLoader *loader = nullptr; + // If the lat/lng can't be represented exactly in floating point, + // there's sometimes roundoff error when converting to integers for + // tile filenames. Adding a little slop prevents truncation. A + // tile is not going to be smaller than 0.1 degrees or so, so this + // amount should be safe. + const float epsilon = 0.001f; + minLat += (minLat >= 0) ? epsilon : -epsilon; + minLng += (minLng >= 0) ? epsilon : -epsilon; + switch (mFileFormat.value()) { case FileFormat::Value::HGT: loader = new HgtLoader(); From 0222b441962f9d089d2810862263a00d6f74b27b Mon Sep 17 00:00:00 2001 From: Andrew Kirmse Date: Tue, 12 Mar 2024 20:10:36 -0700 Subject: [PATCH 2/2] Fix floating point error accumulation. --- code/degree_coordinate_system.cpp | 4 ---- code/flt_loader.cpp | 6 +----- code/prominence.cpp | 17 +++++++++++------ code/prominence_task.cpp | 10 ++++------ code/tile_loading_policy.cpp | 9 +++------ code/util.cpp | 7 +++++++ code/util.h | 4 ++++ 7 files changed, 30 insertions(+), 27 deletions(-) diff --git a/code/degree_coordinate_system.cpp b/code/degree_coordinate_system.cpp index 4557b45..27967e2 100644 --- a/code/degree_coordinate_system.cpp +++ b/code/degree_coordinate_system.cpp @@ -82,10 +82,6 @@ Offsets DegreeCoordinateSystem::offsetsTo(const CoordinateSystem &that) { int dx = static_cast(std::round((mMinLongitude - other->mMinLongitude) * mSamplesPerDegreeLongitude)); int dy = static_cast(std::round((other->mMaxLatitude - mMaxLatitude) * mSamplesPerDegreeLatitude)); - // XXX - VLOG(1) << "Ours " << toString(); - VLOG(1) << "Theirs " << that.toString(); - VLOG(1) << "dx " << dx << " dy " << dy; return Offsets(dx, dy); } diff --git a/code/flt_loader.cpp b/code/flt_loader.cpp index 86c7798..a360a02 100644 --- a/code/flt_loader.cpp +++ b/code/flt_loader.cpp @@ -62,13 +62,9 @@ Tile *FltLoader::loadFromFltFile(const string &directory, float minLat, float mi filename = directory + "/" + filename; } - // XXX - VLOG(1) << "Loading FLT file " << filename; - FILE *infile = fopen(filename.c_str(), "rb"); if (infile == nullptr) { -// XXX VLOG(3) << "Failed to open file " << filename; - VLOG(1) << "Failed to open file " << filename; + VLOG(3) << "Failed to open file " << filename; return nullptr; } diff --git a/code/prominence.cpp b/code/prominence.cpp index 64f38d9..f6d6f39 100644 --- a/code/prominence.cpp +++ b/code/prominence.cpp @@ -198,9 +198,10 @@ int main(int argc, char **argv) { auto threadPool = std::make_unique(numThreads); int num_tiles_processed = 0; vector> results; - float lat = bounds[0]; + // Use double precision to avoid accumulating floating-point error during loop + double lat = bounds[0]; while (lat < bounds[1]) { - float lng = bounds[2]; + double lng = bounds[2]; while (lng < bounds[3]) { // Allow specifying longitude ranges that span the antimeridian (lng > 180) auto wrappedLng = lng; @@ -209,17 +210,21 @@ int main(int argc, char **argv) { } std::shared_ptr coordinateSystem( - fileFormat.coordinateSystemForOrigin(lat, wrappedLng, utmZone)); + fileFormat.coordinateSystemForOrigin( + static_cast(lat), static_cast(wrappedLng), utmZone)); // Skip tiles that don't intersect filtering polygon - if (!filter.intersects(lat, lat + fileFormat.degreesAcross(), - lng, lng + fileFormat.degreesAcross())) { + if (!filter.intersects(static_cast(lat), + static_cast(lat + fileFormat.degreesAcross()), + static_cast(lng), + static_cast(lng + fileFormat.degreesAcross()))) { VLOG(3) << "Skipping tile that doesn't intersect polygon " << lat << " " << lng; } else { // Actually calculate prominence ProminenceTask *task = new ProminenceTask(cache.get(), options); results.push_back(threadPool->enqueue([=] { - return task->run(lat, wrappedLng, *coordinateSystem); + return task->run(static_cast(lat), static_cast(wrappedLng), + *coordinateSystem); })); } diff --git a/code/prominence_task.cpp b/code/prominence_task.cpp index 8e2fdd3..1068806 100644 --- a/code/prominence_task.cpp +++ b/code/prominence_task.cpp @@ -27,6 +27,7 @@ #include "divide_tree.h" #include "island_tree.h" #include "tree_builder.h" +#include "util.h" #include "easylogging++.h" @@ -116,12 +117,9 @@ bool ProminenceTask::writeStringToOutputFile(const string &filename, const strin string ProminenceTask::getFilenamePrefix() const { char filename[PATH_MAX]; - // XXX - const float epsilon = 0.001f; - float lat = mCurrentLatitude; - float lng = mCurrentLongitude; - lat += (lat >= 0) ? epsilon : -epsilon; - lng += (lng >= 0) ? epsilon : -epsilon; + // Deal with floating point imprecision + float lat = adjustCoordinate(mCurrentLatitude); + float lng = adjustCoordinate(mCurrentLongitude); int latHundredths = fractionalDegree(lat); int lngHundredths = fractionalDegree(lng); diff --git a/code/tile_loading_policy.cpp b/code/tile_loading_policy.cpp index cea3e4f..ab07f72 100644 --- a/code/tile_loading_policy.cpp +++ b/code/tile_loading_policy.cpp @@ -105,12 +105,9 @@ Tile *BasicTileLoadingPolicy::loadInternal(float minLat, float minLng) const { // If the lat/lng can't be represented exactly in floating point, // there's sometimes roundoff error when converting to integers for - // tile filenames. Adding a little slop prevents truncation. A - // tile is not going to be smaller than 0.1 degrees or so, so this - // amount should be safe. - const float epsilon = 0.001f; - minLat += (minLat >= 0) ? epsilon : -epsilon; - minLng += (minLng >= 0) ? epsilon : -epsilon; + // tile filenames. Adding a little slop prevents truncation. + minLat = adjustCoordinate(minLat); + minLng = adjustCoordinate(minLng); switch (mFileFormat.value()) { case FileFormat::Value::HGT: diff --git a/code/util.cpp b/code/util.cpp index 93dc8d0..aa4a852 100644 --- a/code/util.cpp +++ b/code/util.cpp @@ -44,6 +44,13 @@ float metersToFeet(float meters) { return meters / 0.3048f; } +float adjustCoordinate(float coordinate) { + // A tile is not going to be smaller than 0.1 degrees or so, so this + // amount should be safe. + const float epsilon = 0.001f; + return coordinate + ((coordinate >= 0) ? epsilon : -epsilon); +} + string trim(const string &s) { static const std::string whitespace = " \t\f\v\n\r"; std::size_t start = s.find_first_not_of(whitespace); diff --git a/code/util.h b/code/util.h index 201b35b..485a4cf 100644 --- a/code/util.h +++ b/code/util.h @@ -36,6 +36,10 @@ float metersToFeet(float meters); // Trim whitespace from the start and end of the string, return a new string std::string trim(const std::string &s); +// Adjust the given coordinate by an epsilon value away from 0, so that truncation +// to int doesn't give incorrect values to due floating-point imprecision. +float adjustCoordinate(float coordinate); + /* * Split the given string by the given delimiter, putting the * result in elems and returning it.