Improve floating point issues with tile/lat conversion

[Helium314]

This greatly reduces the number of LatLons where

StreetComplete/app/src/test/java/de/westnordost/streetcomplete/data/download/tiles/TilesRectTest.kt

Line 13 in 5cc6506

fun `convert bbox to tiles rect and back results in same bbox`() {

fails

It's a draft because the change to tile2lat causes the test for p(0.0, 0.0) to fail. This is some sort of "special" point that's used in many other tests, and thus I think it really needs to work.
So far the only way I found it to work for p(0.0, 0.0) is treating y == numTiles(zoom) / 2 as a special case in tile2lat, though this looks a bit ugly.

Slightly modified convert bbox to tiles rect and back results in same bbox using 1 million random coordinates (with lat limited to range from -85 to 85) fail count:

• master: ca 80k
• with lat2tile nextUp: ca 53k
• with tile2lat nextDown: ca 9200
• with tile2lat nextDown and moving / PI: ca 6700
• with all improvements (this PR): ca 2900

[westnordost]

I lost track of the issue. What issues does it cause exactly if converting forth and back between bbox and tile position does not work in 100% of cases? Also, where (and why) is converting forth and back between bbox and tile position necessary?

[Helium314]

This mainly affects the cache, which uses tiles internally but bboxes otherwise
The issue is not horrible, but frequently leads to needless dropping and subsequent reloading of tiles in

StreetComplete/app/src/main/java/de/westnordost/streetcomplete/util/SpatialCache.kt

Lines 77 to 93 in 5cc6506

	/**
	* Replaces all tiles fully contained in the [bbox] with empty tiles.
	* Any tile only partially contained in the [bbox] is removed (because all tiles cached must be
	* complete).
	* The given [items] are added to cache if the containing tile is cached. Note that for putting
	* and item, it does not have to be inside the [bbox].
	* If the number of tiles exceeds maxTiles, only maxSize tiles are cached.
	*/
	fun replaceAllInBBox(items: Collection<T>, bbox: BoundingBox) { synchronized(this) {
	val tiles = bbox.asListOfEnclosingTilePos()
	val (completelyContainedTiles, incompleteTiles) = tiles.partition { it.asBoundingBox(tileZoom).isCompletelyInside(bbox) }
	if (incompleteTiles.isNotEmpty()) {
	Log.w(TAG, "bbox does not align with tiles, clearing incomplete tiles from cache")
	for (tile in incompleteTiles) {
	removeTile(tile)
	}
	}

So depending on your latitude, you may get worse performance than expected.

[westnordost]

Ah, I see.

Hm well, I do not understand how nextDown / nextUp solves things. This formula,

180.0 / PI * atan(sinh(PI * (1.0 - 2.0 * y / (1 shl zoom))))

should be 100% correct minus the floating point imprecision. So, subtracting or adding a close to infinite amount from or to it should actually make it less correct.

Actually, I didn't know that nextDown and nextUp existed. I would say, those belong in TilePos.asBoundingBox and TilesRect.asBoundingBox in place of the 1e-12s (and maybe on the other way around, bbox to tilerect).

[Helium314]

180.0 / PI * atan(sinh(PI * (1.0 - 2.0 * y / (1 shl zoom))))

should be 100% correct minus the floating point imprecision. So, subtracting or adding a close to infinite amount from or to it should actually make it less correct.

Some comparison with higher precision numbers (80 bit using python / numpy):
Of all z15 y-tiles, converting tile2lat and lat2tile results in a different / wrong tile in 5716 cases (out of 32768). In Python it's slightly different with 5583 failures.
Comparing the tile2lat conversion for these 5583 failures: in 3521 cases the current formula gives the result more close to the 80 bit precision number, in 1983 cases the changes in this PR are more precise. So interestingly adding this little bit makes it more correct in ~36% of those cases.
(tile2lat -> lat2tile fails in only 138 z15 cases with this PR)

Actually, I didn't know that nextDown and nextUp existed. I would say, those belong in TilePos.asBoundingBox and TilesRect.asBoundingBox in place of the 1e-12s (and maybe on the other way around, bbox to tilerect).

When doing this, with lat close to 180 there are issues in:

val p = LatLon(latitude=81.07546376846332, longitude=179.99321275770234)
val tile = p.enclosingTilePos(15)
val bbox = tile.asBoundingBox(15)
val r = bbox.enclosingTilesRect(15)

require(left <= right && top <= bottom) fails in TilesRect.init

[westnordost]

Hm, a little difficult when the communication happens with such delays. Then you have to "reinfuchsen" in the matter each time you respond anew.

Anyway, as you see I am reluctant to merge a deviation to the standard formula that everybody uses as this is obviously the mathematically correct one. You have run a lot of tests and you came to the conclusion that it is more precise in edge cases, but I am missing the "why" here.

When doing this, with lat close to 180 there are issues in:

val p = LatLon(latitude=81.07546376846332, longitude=179.99321275770234)
val tile = p.enclosingTilePos(15)
val bbox = tile.asBoundingBox(15)
val r = bbox.enclosingTilesRect(15)

How can there be issues with nextUp if all it does is to replace 1e-12 with basically the same but a more sophisticated version? Or has it been a problem with the current solution too?

So, I would really like to replace the 1e-12 with nextUp / nextDown, but you are saying that this makes the rounding imprecision worse?

[Helium314]

Hm, a little difficult when the communication happens with such delays. Then you have to "reinfuchsen" in the matter each time you respond anew.

Sorry, finding a good way to compare with higher precision numbers took a while, and I had some other things to do I gave higher priority.

I am missing the "why" here

I honestly don't know. My initial idea was that there is some precision issue amplifying towards one direction in atan(sinh(...)), but this doesn't really work out...

How can there be issues with nextUp if all it does is to replace 1e-12 with basically the same but a more sophisticated version? Or has it been a problem with the current solution too?

No, with 1e-12 it works normally

So, I would really like to replace the 1e-12 with nextUp / nextDown, but you are saying that this makes the rounding imprecision worse?

I didn't yet investigate what causes the require(left <= right && top <= bottom) to fail, will have a look

[Helium314]

will have a look

The reason is in LatLon.enclosingTilePos:
179.99999999999997 + 180 = 360.0 (floating point precision), and 360 % 360 - 180 = -180
So the TilePos unexpectedly is 0

[westnordost]

Where does 179.99999999999997 come from?

[Helium314]

That's one of the longitudes from val bbox = tile.asBoundingBox(15)

[Helium314]

Anyway, I'm considering whether properly investiating all this is really worth more time... maybe I'll just close the PR

[westnordost]

179.99999999999997 + 180 = 360.0 (floating point precision), and 360 % 360 - 180 = -180

In another test, nextDown turned 9.0087890625 into 9.008789062499998 (-1e-15). So, the answer is that the deviation from the next tile through nextDown / nextUp is so small that it gets lost in the calculation that follows.

Looking again into the topic and why it is necessary to calculate close to the floating point precision, the only point where this is important is for the SpatialCache to see if new tiles need to be put into the cache on a bbox query.

So I wonder if the more stable solution wouldn't be to drop the -1e12 logic in TilesRect and do something like

- val (completelyContainedTiles, incompleteTiles) = 
-        tiles.partition { it.asBoundingBox(tileZoom).isCompletelyInside(bbox) }
+ val (completelyContainedTiles, incompleteTiles) = 
+        tiles.partition { it.asBoundingBox(tileZoom).toOsmPrecision().isCompletelyInside(bbox) }

...where toOsmPrecision() shrinks the bounding box to have 7 decimal places (OSM precision) for each coordinate, i.e.

BoundingBox(
  min=LatLon(latitude=52.99495027026899, longitude=8.997802734375)
  max=LatLon(latitude=53.001562274591464, longitude=9.008789062499)
)

becomes

BoundingBox(
  min=LatLon(latitude=52.9949503, longitude=8.9978028)
  max=LatLon(latitude=53.0015622, longitude=9.0087890)
)

Consideration: The gaps between the bounding boxes (up to 10cm IIRC) shrunk this way do not matter for application in the OSM world, because 7 digits is already the maximum precision, i.e. no elements would be outside any bounding box. At the same time, the bounding boxes thus are so well within tiles (compared to what would be near floating point precision), that imprecision through calculations on floating points should not matter anymore.

[Helium314]

I agree, then we don't have to care about precision issues.
Only thing where this should be considered is when adding / moving nodes. As far as I know the new coordinates aren't truncated to OSM precision.