add fallback strategy for tinshift transform to use closest triangle for points not in any

[josch]

Hi,

this is a WIP pull request for #2903 and I'm asking for guidance to finish it up. Open questions are:

• is "fallback_strategy" a good name? The default does not change the current behaviour.
• should "fallback_strategy" be part of the json or the proj string?
• currently "fallback_strategy" is a string -- should it be an enum? Maybe there are other strategies that could be useful in future implementations?
• should the json version be bumped because of the new field?
• is there a way to compute distance to a triangle with less computations? I can only come up with heuristics like the commented out piece of code which is not precise but only checks the vertices of the triangles.
• the quadtree does not seem to be ready for performing a nearest-neighbor kind of search. It would need to be extended to support a distance function and a shared variable to track the currently found closest feature. For my own use-case, the current implementation that does not use the quad tree is sufficient.

TODO:

• Tests added
• Added clear title that can be used to generate release notes
• Fully documented, including updating docs/source/*.rst for new API

[rouault]

• is "fallback_strategy" a good name? The default does not change the current behaviour.

Sounds good to me

• should "fallback_strategy" be part of the json or the proj string?

I'm fine with it in the json.

• currently "fallback_strategy" is a string -- should it be an enum? Maybe there are other strategies that could be useful in future implementations?

In the C++ code you mean I guess ?(as you already made it an enum in the json schema) . Yes that would be fine

• should the json version be bumped because of the new field?

Yes, adding 1.1 in the format_version enum of the json schema

• is there a way to compute distance to a triangle with less computations? I can only come up with heuristics like the commented out piece of code which is not precise but only checks the vertices of the triangles.

I don't really have ready-made suggestions here. This is certainly a well-known problem, and must be implemented in GEOS for the Distance() function.

• the quadtree does not seem to be ready for performing a nearest-neighbor kind of search. It would need to be extended to support a distance function and a shared variable to track the currently found closest feature. For my own use-case, the current implementation that does not use the quad tree is sufficient.

That would probably be fine to use that "naive" approach for now.

[josch]

Thanks a lot for your input! I will work on these points and push a new version once I'm ready.

• is there a way to compute distance to a triangle with less computations? I can only come up with heuristics like the commented out piece of code which is not precise but only checks the vertices of the triangles.

I don't really have ready-made suggestions here. This is certainly a well-known problem, and must be implemented in GEOS for the Distance() function.

Can you elaborate on this point? This is my first contact with PROJ and I don't know what you mean with GEOS and the Distance() function.

Thanks!

[rouault]

Can you elaborate on this point? This is my first contact with PROJ and I don't know what you mean with GEOS and the Distance() function.

I was thinking to the GEOS library that has gazillons of geometry computation algorithms : https://github.com/libgeos/geos/blob/main/capi/geos_c.h.in#L1520 . But we don't want to depend or draw GEOS in PROJ. Maybe https://stackoverflow.com/questions/849211/shortest-distance-between-a-point-and-a-line-segment could be of help ?

[josch]

I was thinking to the GEOS library that has gazillons of geometry computation algorithms : https://github.com/libgeos/geos/blob/main/capi/geos_c.h.in#L1520 . But we don't want to depend or draw GEOS in PROJ.

Ah, now I understand. You were suggesting that maybe GEOS has a better implementation. I just checked and no it does not: https://github.com/libgeos/geos/blob/5a3cebe6e9b1b1d612638a07fce6d81ac04cd684/src/operation/distance/DistanceOp.cpp#L601 Just as my code above, GEOS just computes the distance between line segments and points. Its documentation even warns that its approach is O(n²): http://geos.refractions.net/ro/doxygen_docs/html/classgeos_1_1operation_1_1distance_1_1DistanceOp.html#_details

Maybe https://stackoverflow.com/questions/849211/shortest-distance-between-a-point-and-a-line-segment could be of help ?

The answers there do the same thing that my code does: build a line from two points, obtain the line parameter t by projecting the point onto the line and check whether t is less than 0, larger than 1 or in the middle to compute the distance.

[rouault]

Beyond the few details pointed, overall looks good to me !

[busstoptaktik]

currently "fallback_strategy" is a string -- should it be an enum? Maybe there are other strategies that could be useful in future implementations?

Sorry for being late to the party, but I think you are right here: Actually I think, the obvious thing to do is to select the triangle for which the external point is closest to its centroid (i.e. the intersection of the medians).

In many cases, the triangle selected will not differ from "nearest side", since geodetic TINs tend to be conscoiusly and sensibly constructed. But in the case of a long narrow triangle with the short side facing the outer edge of the TIN, its neighbouring more harmoniously shaped triangles may be more fit for extrapolation?

And if having the centroids of all triangles precomputed at instantiation, finding the nearest (in the centroid sense) is simply a matter of comparing squared point-to -point distances (x - x_n) * (x - x_n) + (y - y_n) * (y - y_n)

[josch]

Sorry for being late to the party, but I think you are right here: Actually I think, the obvious thing to do is to select the triangle for which the external point is closest to its centroid (i.e. the intersection of the medians).

This is a good point and it also means that "nearest" is a bad name for the format. I renamed it to "nearest_side" and added a new method called "nearest_centroid" which finds the triangle with the closest centroid to the query point.