Cutout bounds fix

cadc-data-ops-fits/README.md

@@ -5,12 +5,20 @@ The `cadc-data-ops-fits` library depends on the NASA led NOM TAM FITS library
 (https://github.com/nom-tam-fits/nom-tam-fits) version 1.15.3 (or newer).
 
 ## Building it
-You may use the provided Gradle Wrapper, or provide your own Gradle (< 7) installation.
+Use the repository Gradle Wrapper (8.x). From the repo root:
 
 ```sh
-$ ../gradlew -i clean build
+../gradlew -i :cadc-data-ops-fits:clean :cadc-data-ops-fits:build
 ```
 
+### WCS native dependency (`cadc-wcs`)
+
+World-coordinate cutouts use **`cadc-wcs`**, which loads a JNI library bundled in the JAR plus the system **WCSLib** C library at runtime. **`mavenLocal()` is listed before `mavenCentral()`** in this module so you can **`publishToMavenLocal`** a locally built `cadc-wcs` (with a JNI binary for your OS) and override the artifact from Central.
+
+**Linux:** Install WCSLib from your distribution (e.g. `wcslib-dev` on Debian/Ubuntu) or a prefix build. **macOS:** Install via Homebrew, MacPorts, or a prefix build ([WCSLIB](https://www.atnf.csiro.au/computing/software/wcs/)).
+
+If tests fail with `WCSLibInitializationException` or `UnsatisfiedLinkError`, build JNI in the [`opencadc/wcs`](https://github.com/opencadc/wcs) `cadc-wcs` project (`./gradlew -c settings-jni.gradle copyJniToResources` after installing WCSLib), then `./gradlew :cadc-wcs:publishToMavenLocal` from that repo. Optional overrides when building JNI: environment variables **`WCSLIB_LIB`** (path to `libwcs.so` / `libwcs.dylib`) or **`WCSLIB_LIB_DIR`**, or Gradle **`-Pwcslib.lib=...` / `-Pwcslib.libDir=...`**.
+
 ## Cutout API
 This library supports the commonly used cutout syntax to extract a sub-image from an Image HDU.
 

cadc-data-ops-fits/build.gradle

@@ -11,21 +11,24 @@ repositories {
 
 apply from: '../opencadc.gradle'
 
-sourceCompatibility = 11
+java {
+    sourceCompatibility = JavaVersion.VERSION_11
+    targetCompatibility = JavaVersion.VERSION_11
+}
 
 group = 'org.opencadc'
-version = '0.4.1'
+version = '0.4.2'
 
 description = 'OpenCADC FITS cutout library'
 def git_url = 'https://github.com/opencadc/dal'
 
 dependencies {
     implementation 'org.opencadc:cadc-dali:[1.2.10,2.0.0)'
-    implementation 'org.opencadc:cadc-util:[1.6,2.0)'
+    implementation 'org.opencadc:cadc-util:[1.12.5,2.0)'
     implementation 'org.opencadc:cadc-soda-server:[1.2.1,2.0)'
-    implementation 'org.opencadc:cadc-wcs:[2.1.4,3.0)'
+    implementation 'org.opencadc:cadc-wcs:[2.2.0,2.4.0)'
     implementation 'org.opencadc:jsky:[1.0.0,2.0.0)'
-    implementation 'gov.nasa.gsfc.heasarc:nom-tam-fits:1.20.0'
+    implementation 'gov.nasa.gsfc.heasarc:nom-tam-fits:1.22.0'
 
     // Use JUnit test framework
     testImplementation 'junit:junit:[4.13,5.0)'

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

@@ -0,0 +1,42 @@
+package org.opencadc.fits.slice;
+
+import java.util.Arrays;
+
+import nom.tam.fits.header.Standard;
+import org.apache.log4j.Logger;
+
+
+/**
+ * Fills a flat {@code long[naxis*2]} bounds array: the cut axis uses clipped interval coordinates;
+ * every other axis spans {@code [1, NAXISn]} for that axis (not a single global length).
+ */
+public class AxisBoundsFiller {
+    private static final Logger LOGGER = Logger.getLogger(AxisBoundsFiller.class);
+
+    static long[] fill(final int axes, final long[] clippedBounds, final int clipAxis, final int[] naxisPerAxis) {
+        LOGGER.debug("Filling bounds for " + axes + " axes, clip axis " + clipAxis + ", clipped bounds "
+                     + Arrays.toString(clippedBounds));
+        final long[] bounds = new long[axes];
+        for (int i = 0; i < axes; i += 2) {
+            final int axis = (i + 2) / 2;
+            if (axis == clipAxis) {
+                bounds[i] = clippedBounds.length > 0 ? clippedBounds[0] : 1L;
+                bounds[i + 1] = clippedBounds.length > 1 ? clippedBounds[1] : naxisPerAxis[axis - 1];
+            } else {
+                bounds[i] = 1L;
+                bounds[i + 1] = naxisPerAxis[axis - 1];
+            }
+        }
+        LOGGER.debug("Filled bounds: " + Arrays.toString(bounds));
+
+        return bounds;
+    }
+
+    static int[] naxisSizes(final FITSHeaderWCSKeywords wcs, final int naxis) {
+        final int[] sizes = new int[naxis];
+        for (int a = 1; a <= naxis; a++) {
+            sizes[a - 1] = wcs.getIntValue(Standard.NAXISn.n(a).key());
+        }
+        return sizes;
+    }
+}

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

@@ -133,10 +133,20 @@ public long[] getBounds(final Interval<Number> bounds) throws NoSuchKeywordExcep
             return null;
         } else {
             final int naxis = spectralWCSKeywords.getIntValue(Standard.NAXIS.key());
-            final Interval<Double> boundsIntervalPixel = getCutoutPixelInterval(bounds, energyAxis, naxis);
-            final Interval<Double> nativePixelsInterval =
-                    new Interval<>(0.0D, (double) spectralWCSKeywords.getIntValue(
-                            Standard.NAXISn.n(energyAxis).key()));
+            final String ctype = spectralWCSKeywords.getStringValue(Standard.CTYPEn.n(energyAxis).key());
+            final boolean isVelocity = CoordTypeCode.fromCType(ctype).isVelocity();
+            // Intersect user-requested wavelength (m) with the spectral range covered by the data at pixels 1..N
+            // so that extended / infinite physical bounds (e.g. -Inf) map to a finite WCS call instead of
+            // sky2pix endpoints that miss the field of view in pixel space.
+            final Interval<Number> requestMetres = isVelocity ? bounds : clampWavelengthToSpectralFieldOfView(bounds, energyAxis);
+            if (requestMetres == null) {
+                LOGGER.warn("No overlap with spectral extent on data.");
+                return null;
+            }
+            final int nchan = spectralWCSKeywords.getIntValue(Standard.NAXISn.n(energyAxis).key());
+            final Interval<Double> boundsIntervalPixel = getCutoutPixelInterval(requestMetres, energyAxis, naxis);
+            // FITS pixel indices 1..N; clip() in FITSCutout is 1-based to len inclusive.
+            final Interval<Double> nativePixelsInterval = new Interval<>(1.0D, (double) nchan);
             final Interval<Double> intersectionPixels = getOverlap(nativePixelsInterval, boundsIntervalPixel);
 
             if (intersectionPixels == null) {
@@ -152,39 +162,65 @@ public long[] getBounds(final Interval<Number> bounds) throws NoSuchKeywordExcep
                         clip(maxSpectralLength, (long) Math.floor(Math.min(low, up) + 0.5D),
                              (long) Math.ceil(Math.max(low, up) - 0.5D));
 
-                final long[] entireBounds = clippedSpectralBounds == null ? null : new long[naxis * 2];
-
-                if (entireBounds != null) {
-                    for (int i = 0; i < entireBounds.length; i += 2) {
-                        final int axis = (i + 2) / 2;
-                        if (axis == energyAxis) {
-                            entireBounds[i] = clippedSpectralBounds[0];
-                            entireBounds[i + 1] = clippedSpectralBounds[1];
-                        } else {
-                            entireBounds[i] = 1L;
-                            entireBounds[i + 1] = (long) this.fitsHeaderWCSKeywords.getDoubleValue(
-                                    Standard.NAXISn.n(axis).key());
-                        }
-                    }
-                }
-
-                return entireBounds;
+                final int axes = clippedSpectralBounds == null ? 0 : naxis * 2;
+                return AxisBoundsFiller.fill(axes, clippedSpectralBounds, energyAxis,
+                                             AxisBoundsFiller.naxisSizes(spectralWCSKeywords, naxis));
             }
         }
     }
 
-    private Interval<Double> getOverlap(final Interval<Double> headerWCSInterval, final Interval<Double> cutoutBounds) {
-        LOGGER.debug("Checking overlap between header pixels ("
-                     + headerWCSInterval.getLower() + ", " + headerWCSInterval.getUpper()
-                     + ") and requested bounds pixels ("
-                     + cutoutBounds.getLower() + ", " + cutoutBounds.getUpper() + ")");
-        if (headerWCSInterval.getLower() > cutoutBounds.getUpper()
-            || headerWCSInterval.getUpper() < cutoutBounds.getLower()) {
+    private Interval<Double> getOverlap(final Interval<Double> a, final Interval<Double> b) {
+        final double lo = Math.max(a.getLower(), b.getLower());
+        final double hi = Math.min(a.getUpper(), b.getUpper());
+        LOGGER.debug("Pixel interval intersection (" + a.getLower() + ", " + a.getUpper() + ") with ("
+                     + b.getLower() + ", " + b.getUpper() + ") -> (" + lo + ", " + hi + ")");
+        if (lo > hi) {
+            return null;
+        }
+        return new Interval<>(lo, hi);
+    }
+
+    /**
+     * Wavelength in metres (barycentric) for channels 1 and nchan, using a linear WCS in the native
+     * spectral unit (CUNIT) at the reference pixel (same approximation as a simple 1D grid).
+     * Returns the ordered pair [min, max] in metres.
+     */
+    private double[] spectralWavelengthMetresAtBandEdges(final int energyAxis) {
+        final int nchan = this.fitsHeaderWCSKeywords.getIntValue(Standard.NAXISn.n(energyAxis).key());
+        final double crval = this.fitsHeaderWCSKeywords.getDoubleValue(Standard.CRVALn.n(energyAxis).key());
+        final double cdelt = this.fitsHeaderWCSKeywords.getDoubleValue(Standard.CDELTn.n(energyAxis).key(), 0.0D);
+        final double crpix = this.fitsHeaderWCSKeywords.getDoubleValue(Standard.CRPIXn.n(energyAxis).key());
+        final String cunit = this.fitsHeaderWCSKeywords.getStringValue(CADCExt.CUNITn.n(energyAxis).key());
+        final double world1 = crval + cdelt * (1.0D - crpix);
+        final double worldN = crval + cdelt * ((double) nchan - crpix);
+        final EnergyConverter energyConverter = new EnergyConverter();
+        final double m1 = energyConverter.toMeters(world1, cunit);
+        final double m2 = energyConverter.toMeters(worldN, cunit);
+        return new double[]{Math.min(m1, m2), Math.max(m1, m2)};
+    }
+
+    /**
+     * Intersect the requested barycentric wavelength range (m) with the [min,max] in metres of the
+     * spectral band on the data (linear in native unit at pixels 1 and nchan). Handles ±Infinity on
+     * the request via {@code max}/{@code min} with the band; returns null if there is no intersection.
+     */
+    private Interval<Number> clampWavelengthToSpectralFieldOfView(final Interval<Number> bounds, final int energyAxis) {
+        final double wmin = bounds.getLower().doubleValue();
+        final double wmax = bounds.getUpper().doubleValue();
+        final double rlo = Math.min(wmin, wmax);
+        final double rhi = Math.max(wmin, wmax);
+        final double[] m = spectralWavelengthMetresAtBandEdges(energyAxis);
+        final double mmin = m[0];
+        final double mmax = m[1];
+        // max(rlo, mmin) lets -Inf select the in-band start; min(rhi, mmax) clips a request that runs past the band.
+        final double cLo = Math.max(rlo, mmin);
+        final double cHi = Math.min(rhi, mmax);
+        LOGGER.debug("Spectral (m) request [" + rlo + ", " + rhi + "] with header band ~[" + mmin + ", " + mmax
+                     + "] -> [" + cLo + ", " + cHi + "]");
+        if (cLo > cHi) {
             return null;
-        } else {
-            return new Interval<>(Math.max(headerWCSInterval.getLower(), cutoutBounds.getLower()),
-                                  Math.min(headerWCSInterval.getUpper(), cutoutBounds.getUpper()));
         }
+        return new Interval<>(cLo, cHi);
     }
 
     /**