fix(pq): resolve Peirce Quincuncial inverse singularities using complex Jacobi elliptic functions - #4811
fix(pq): resolve Peirce Quincuncial inverse singularities using complex Jacobi elliptic functions#4811dulong-lab wants to merge 1 commit into
Conversation
f3d41f5 to
275d801
Compare
275d801 to
8dbc8bb
Compare
fallenmi
left a comment
There was a problem hiding this comment.
Disclosure: I used Codex to assist with this review and verified the finding against exact commit 8dbc8bbf2c1b471593ffe164947ad27e9f8094ad.
Please undo scrolly with the opposite sign in the vertical inverse. The forward vertical path adds Q->scrolly * (4 * K) to xy.y (modulo the 4K wrap), but peirce_q_inverse() adds that same offset again at line 415. This breaks round trips for most nonzero vertical scroll values.
For example, with +proj=peirce_q +shape=vertical +scrolly=0.25 +R=1, exact head maps (10°, 20°) forward to (0.220422825640, -1.357737865642) and then inversely to (170.000000327°, -20.000001266°) instead of the original point. In a 72-point by 13-scroll oracle, exact head/current merge pass only 360/936 vertical round trips; every case at scrolly=±0.125, ±0.25, ±0.375, and ±0.75 fails. The analogous horizontal-scroll control passes 936/936.
Changing only this expression to xy.y + K - Q->scrolly * (4 * K) makes the vertical oracle pass 936/936, with maximum error below 0.0000023 degrees. The focused Peirce GIE suite passes 592/592, but its scrolled-vertical section is forward-only. Please subtract the offset and add a nonzero-scrolly inverse regression.
docs/source/*.rstfor new APIThe current implementation of the pq projection relies on pj_generic_inverse_2d for its inverse transformation. This approach fails to converge near the singularities.
This PR replaces the iterative solver with an analytical solution using complex Jacobi elliptic functions.
Visual Verification
The comparison images were generated by performing an inverse projection on a coordinate grid and sampling colors from a global TIF file.
Notes for Maintainers
Domain Checking: I haven't implemented strict input coordinate validation. Since the Peirce Quincuncial projection can tile the entire plane infinitely, I'm unsure what the standard PROJ practice is for handling inputs outside the primary square.
Test accuracy adjustment: I have updated the .gie test files to use a tolerance 350 mm for both shape=square and shape=diamond operations.
Additional Test Cases: While horizontal, vertical, nhemisphere, and shemisphere shapes have been visually verified to work correctly, I have not added formal .gie test cases for them. I am not familiar with the process of generating reference data for PROJ tests, and it would be more efficient for maintainers to add them if needed.
AI Use: AI was used to understanding the math of Peirce Quincuncial projection, assist in translating the technical description into English and for discussing C++ best practices. The core implementation of ellpj_5 was adapted from SciPy's Cephes library.