Merge pull request #104 from at88mph/incorrect-wcs

Incorrect wcs

cadc-data-ops-fits/build.gradle

@@ -12,7 +12,7 @@ repositories {
 
 sourceCompatibility = 11
 group = 'org.opencadc'
-version = '0.4.0'
+version = '0.4.1'
 
 description = 'OpenCADC FITS cutout library'
 def git_url = 'https://github.com/opencadc/dal'

cadc-data-ops-fits/src/main/java/org/opencadc/fits/slice/WCSCutoutUtil.java

@@ -81,7 +81,9 @@
 import java.util.List;
 
 import nom.tam.fits.Header;
+import nom.tam.fits.HeaderCard;
 import nom.tam.fits.HeaderCardException;
+import nom.tam.fits.header.Standard;
 import org.apache.log4j.Logger;
 import org.opencadc.soda.PixelRange;
 import org.opencadc.soda.server.Cutout;
@@ -138,6 +140,55 @@ public static PixelRange[] getBounds(final Header header, final Cutout cutout)
         return allPixelRanges.toArray(new PixelRange[0]);
     }
 
+    /**
+     * Post cutout routine to adjust necessary Header Card values, such as CRPIXn, to match the dimensions of the resulting cutout.
+     * @param header            The Header to adjust.
+     * @param dimensionLength   The length of the dimensions to iterate.
+     * @param corners           The corners (starting co-ordinates) of the resulting image.
+     * @param steps             The striding value to skip while reading.
+     */
+    static void adjustHeaders(final Header header, final int dimensionLength, final int[] corners, final int[] steps) {
+        // CRPIX values are not set automatically.  Adjust them here, if present.
+        for (int i = 0; i < dimensionLength; i++) {
+            final HeaderCard crPixCard = header.findCard(Standard.CRPIXn.n(i + 1));
+            final int stepValue = steps[steps.length - i - 1];
+            if (crPixCard != null) {
+                // Need to run backwards (reverse order) to match the dimensions.
+                final double nextValue = corners[corners.length - i - 1];
+                final double crPixValue = Double.parseDouble(crPixCard.getValue()) - nextValue;
+
+                if (stepValue > 1) {
+                    final double newValue = (crPixValue / stepValue) + (1.0 - (1.0 / stepValue));
+                    crPixCard.setValue(newValue);
+                    LOGGER.debug("Adjusted " + crPixCard.getKey() + " to " + newValue);
+                } else {
+                    crPixCard.setValue(crPixValue);
+                    LOGGER.debug("Set " + crPixCard.getKey() + " to " + crPixValue);
+                }
+            }
+
+            // CDELTn cards are typically present for PC matrices, but due to the nature of Archive data,
+            // they could be present even if a CD matrix is present.  Since PC matrices are the default in
+            // the absence of a CD matrix, it's not unusual for it to be absent.
+            // See https://www.aanda.org/articles/aa/pdf/2002/45/aah3859.pdf
+            //
+            final HeaderCard cDeltCard = header.findCard(Standard.CDELTn.n(i + 1));
+            if (cDeltCard != null) {
+                final double cDeltValue = Double.parseDouble(cDeltCard.getValue());
+                cDeltCard.setValue(cDeltValue * (double) stepValue);
+            }
+
+            // Handle the entire CD matrix.
+            for (int j = 0; j < dimensionLength; j++) {
+                final HeaderCard cdMatrixCard = header.findCard(String.format("CD%d_%d", i + 1, j + 1));
+                if (cdMatrixCard != null) {
+                    final double cdMatrixValue = Double.parseDouble(cdMatrixCard.getValue());
+                    cdMatrixCard.setValue(cdMatrixValue * (double) stepValue);
+                }
+            }
+        }
+    }
+
     static PixelRange[] getSpatialBounds(final Header header, final Shape shape)
             throws HeaderCardException, NoSuchKeywordException {
         final long[] bounds;

cadc-data-ops-fits/src/test/java/org/opencadc/fits/FitsTest.java

@@ -77,6 +77,7 @@
 import nom.tam.fits.HeaderCard;
 import nom.tam.fits.header.IFitsHeader;
 import nom.tam.fits.header.Standard;
+import nom.tam.fits.header.extra.NOAOExt;
 import org.apache.log4j.Level;
 import org.apache.log4j.Logger;
 import org.junit.Assert;
@@ -91,10 +92,10 @@ public class FitsTest {
 
     private static final IFitsHeader[] HEADER_CARD_KEYS_TO_CHECK = new IFitsHeader[]{
             Standard.BITPIX, Standard.NAXIS, Standard.EXTNAME, Standard.XTENSION, Standard.SIMPLE, Standard.EXTVER,
-            Standard.BSCALE, Standard.BUNIT
+            Standard.BSCALE, Standard.BUNIT, NOAOExt.CD1_1, NOAOExt.CD1_2, Standard.CDELTn.n(1), Standard.CRPIXn.n(1)
     };
 
-    public static void assertFitsEqual(final Fits expected, final Fits result) throws Exception {
+    public static void assertFitsEqual(final Fits expected, final Fits result) {
         final BasicHDU<?>[] expectedHDUList = expected.read();
         final BasicHDU<?>[] resultHDUList = result.read();
 
@@ -105,22 +106,17 @@ public static void assertFitsEqual(final Fits expected, final Fits result) throw
             final BasicHDU<?> nextResultHDU = resultHDUList[expectedIndex];
 
             LOGGER.debug("On Extension at index " + expectedIndex);
-            FitsTest.assertHDUEqual(nextExpectedHDU, nextResultHDU);
+            FitsTest.assertHeadersEqual(nextExpectedHDU.getHeader(), nextResultHDU.getHeader());
         }
     }
 
-    public static void assertHDUEqual(final BasicHDU<?> expectedHDU, final BasicHDU<?> resultHDU) throws Exception {
-        final Header expectedHeader = expectedHDU.getHeader();
-        final Header resultHeader = resultHDU.getHeader();
-
-        FitsTest.assertHeadersEqual(expectedHeader, resultHeader);
-    }
-
-    public static void assertHeadersEqual(final Header expectedHeader, final Header resultHeader) throws Exception {
+    public static void assertHeadersEqual(final Header expectedHeader, final Header resultHeader) {
         Arrays.stream(HEADER_CARD_KEYS_TO_CHECK).forEach(headerCardKey -> {
             final HeaderCard expectedCard = expectedHeader.findCard(headerCardKey);
             final HeaderCard resultCard = resultHeader.findCard(headerCardKey);
 
+            LOGGER.info("Checking " + headerCardKey.key());
+
             if (expectedCard == null) {
                 Assert.assertNull("Card " + headerCardKey.key() + " should be null.", resultCard);
             } else {
@@ -131,7 +127,17 @@ public static void assertHeadersEqual(final Header expectedHeader, final Header
                 if (valueType == Float.class) {
                     Assert.assertEquals("Header " + headerCardKey.key() + " has the wrong value.",
                                         Float.parseFloat(expectedCard.getValue()),
-                                        Float.parseFloat(resultCard.getValue()), 0.0F);
+                                        Float.parseFloat(resultCard.getValue()), 1.0e-5F);
+                } else if (valueType == Double.class) {
+                    Assert.assertEquals("Header " + headerCardKey.key() + " has the wrong value.",
+                                        Double.parseDouble(expectedCard.getValue()),
+                                        Double.parseDouble(resultCard.getValue()), 1.0e-5D);
+                } else if (valueType == Integer.class) {
+                    // Expected type has been declared as Integer, but result may have been converted to Float (i.e. 0 == 0e0), so
+                    // allow some robustness here.
+                    Assert.assertEquals("Header " + headerCardKey.key() + " has the wrong value.",
+                                        Integer.parseInt(expectedCard.getValue()),
+                                        Math.round(Float.parseFloat(resultCard.getValue())));
                 } else {
                     Assert.assertEquals("Header " + headerCardKey.key() + " has the wrong value.",
                                         expectedCard.getValue(), resultCard.getValue());

cadc-data-ops-fits/src/test/java/org/opencadc/fits/slice/NDimensionalSlicerTest.java

@@ -85,6 +85,7 @@
 import org.apache.log4j.Level;
 import org.apache.log4j.Logger;
 import org.junit.Assert;
+import org.junit.Ignore;
 import org.junit.Test;
 import org.opencadc.fits.FitsTest;
 import org.opencadc.fits.NoOverlapException;
@@ -101,6 +102,36 @@ public class NDimensionalSlicerTest {
         Log4jInit.setLevel("org.opencadc.fits", Level.DEBUG);
     }
 
+    @Test
+    @Ignore("Requires a larger file to cut from.  Here to illustrate test input.")
+    public void testIncorrectWCS() throws Exception {
+        ExtensionSliceFormat fmt = new ExtensionSliceFormat();
+        List<ExtensionSlice> slices = new ArrayList<>();
+        slices.add(fmt.parse("[*:4,*:6,*:6]"));
+        final Cutout cutout = new Cutout();
+        cutout.pixelCutouts = slices;
+
+        final NDimensionalSlicer slicer = new NDimensionalSlicer();
+        final File file = FileUtil.getFileFromResource("test-alma-cube.fits", NDimensionalSlicerTest.class);
+
+        final String configuredTestWriteDir = System.getenv("TEST_WRITE_DIR");
+        final String testWriteDir = configuredTestWriteDir == null ? "/tmp" : configuredTestWriteDir;
+        final File expectedFile = FileUtil.getFileFromResource("test-alma-cube-cutout.fits",
+                                                               NDimensionalSlicerTest.class);
+        final Path outputPath = Files.createTempFile(new File(testWriteDir).toPath(), "test-alma-cube-cutout", ".fits");
+        LOGGER.debug("Writing out to " + outputPath);
+
+        try (final OutputStream outputStream = Files.newOutputStream(outputPath.toFile().toPath())) {
+            slicer.slice(file, cutout, outputStream);
+        }
+
+        final Fits expectedFits = new Fits(expectedFile);
+        final Fits resultFits = new Fits(outputPath.toFile());
+
+        FitsTest.assertFitsEqual(expectedFits, resultFits);
+        Files.deleteIfExists(outputPath);
+    }
+
     @Test
     public void testMEFFileSlice() throws Exception {
         ExtensionSliceFormat fmt = new ExtensionSliceFormat();

cadc-data-ops-fits/src/test/java/org/opencadc/fits/slice/WCSCutoutUtilTest.java

@@ -75,6 +75,8 @@
 import ca.nrc.cadc.util.Log4jInit;
 import nom.tam.fits.Header;
 import nom.tam.fits.HeaderCard;
+import nom.tam.fits.header.Standard;
+import nom.tam.fits.header.extra.NOAOExt;
 import org.apache.log4j.Level;
 import org.apache.log4j.Logger;
 import org.junit.Assert;
@@ -129,4 +131,34 @@ public void testMultipleWCS() throws Exception {
         }
         LOGGER.debug("WCSCutoutUtilTest.testMultipleWCS OK: " + (System.currentTimeMillis() - startMillis) + " ms");
     }
+
+    @Test
+    public void testWCSAdjustment() throws Exception {
+        final String headerFileName = "test-blast-header-1.txt";
+        final File testFile = FileUtil.getFileFromResource(headerFileName, CircleCutoutTest.class);
+
+        // Corners and striding values MUST be in reverse order as per what nom-tam-fits provides.  This accurately
+        // represents the use case.
+        final int[] corners = new int[]{0, 0};
+        final int[] stridingValues = new int[]{5, 1};
+
+        try (final InputStream inputStream = Files.newInputStream(testFile.toPath());
+             final BufferedReader bufferedReader = new BufferedReader(new InputStreamReader(inputStream))) {
+            final Header testHeader = new Header();
+            final char[] buff = new char[80];
+            while (bufferedReader.read(buff) >= 0) {
+                testHeader.addLine(HeaderCard.create(new String(buff)));
+            }
+
+            final double originalCD11 = testHeader.getDoubleValue(NOAOExt.CD1_1);
+            final double originalCD21 = testHeader.getDoubleValue(NOAOExt.CD2_1);
+
+            WCSCutoutUtil.adjustHeaders(testHeader, testHeader.getIntValue(Standard.NAXIS), corners, stridingValues);
+
+            Assert.assertEquals("Should remain unchanged as only applies to second axis.", originalCD11, testHeader.getDoubleValue(NOAOExt.CD1_1),
+                                0.0D);
+            Assert.assertEquals("Should be adjusted.", originalCD21 * ((double) stridingValues[1]), testHeader.getDoubleValue(NOAOExt.CD2_1),
+                                0.0D);
+        }
+    }
 }