cds-astro / cds-healpix-java

The CDS HEALPix library in Java
BSD 3-Clause "New" or "Revised" License
9 stars 9 forks source link

Numerical error in HealpixNestedPolygonComputer.overlappingCells #18

Closed mbtaylor closed 4 months ago

mbtaylor commented 4 months ago

Hi @fxpineau. I see a problem if I run the following code against cds-healpix-java v0.30.2:

import cds.healpix.Healpix;
import cds.healpix.HealpixNestedPolygonComputer;

public class Overlap {
    public static void main(String[] args) {
        int order = 12;
        HealpixNestedPolygonComputer pc = Healpix.getNested( 12 ).newPolygonComputer();
        double[][] verts1 = new double[][] {
            { -1.5702949547333407, -0.7295093151415473 },
            { -1.5702171673769950, -0.7295093016804524 },
            { -1.5701393800214274, -0.7295092852142693 },
            { -1.5700615926667945, -0.7295092657429985 },
        };
        double[][] verts2 = new double[][] {
            { -1.5706045044233712, -0.7295105218483977 },
            { -1.5705168372776197, -0.7295105199399403 },
            { -1.5704291701320918, -0.7295105142145686 },
            { -1.5703415029870114, -0.7295105046722821 },
        };
        pc.overlappingCells( verts1 );  // OK
        pc.overlappingCells( verts2 );  // trouble
    }
}

The calculation with the second set of vertices gives me an AssertionError at cds.healpix.NewtonMethod.arcSpecialPointInEqr(NewtonMethod.java:234), or a NullPointerException later on if assertions are turned off. This is probably some kind of edge case since this code has been in topcat for several years and I haven't noticed it before now.

fxpineau commented 4 months ago

Hi @mbtaylor ,

Thank you for the report.

The good news is that the rust library seems to be working well for both polygons. See the screenshot made in Aladin: Poly1and2

I am going to check the differences between both libraries to fix the Java one (today or tomorrow).

fx

fxpineau commented 4 months ago

Fixed in commit b6880155961430bc2e97dd3d308f1d1f160703cc (if the Newton method fails, we are in a corner case which probably does not contains "special points"; if we found a corner case having a special point we could fallback with a dichotomic appoach).

mbtaylor commented 4 months ago

Thanks @fxpineau, I confirm that fixes my test case and the application issue I encountered that led to the report. Will you assign a new release number/tag to this version?

fxpineau commented 4 months ago

Sure. I just updated the CHANGES file, the version number and added a tag. Thank you @mbtaylor .

mbtaylor commented 4 months ago

Thanks!