All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.opencadc.fits.slice.WCSCutoutUtil Maven / Gradle / Ivy

The newest version!
/*
 ************************************************************************
 *******************  CANADIAN ASTRONOMY DATA CENTRE  *******************
 **************  CENTRE CANADIEN DE DONNÉES ASTRONOMIQUES  **************
 *
 *  (c) 2021.                            (c) 2021.
 *  Government of Canada                 Gouvernement du Canada
 *  National Research Council            Conseil national de recherches
 *  Ottawa, Canada, K1A 0R6              Ottawa, Canada, K1A 0R6
 *  All rights reserved                  Tous droits réservés
 *
 *  NRC disclaims any warranties,        Le CNRC dénie toute garantie
 *  expressed, implied, or               énoncée, implicite ou légale,
 *  statutory, of any kind with          de quelque nature que ce
 *  respect to the software,             soit, concernant le logiciel,
 *  including without limitation         y compris sans restriction
 *  any warranty of merchantability      toute garantie de valeur
 *  or fitness for a particular          marchande ou de pertinence
 *  purpose. NRC shall not be            pour un usage particulier.
 *  liable in any event for any          Le CNRC ne pourra en aucun cas
 *  damages, whether direct or           être tenu responsable de tout
 *  indirect, special or general,        dommage, direct ou indirect,
 *  consequential or incidental,         particulier ou général,
 *  arising from the use of the          accessoire ou fortuit, résultant
 *  software.  Neither the name          de l'utilisation du logiciel. Ni
 *  of the National Research             le nom du Conseil National de
 *  Council of Canada nor the            Recherches du Canada ni les noms
 *  names of its contributors may        de ses  participants ne peuvent
 *  be used to endorse or promote        être utilisés pour approuver ou
 *  products derived from this           promouvoir les produits dérivés
 *  software without specific prior      de ce logiciel sans autorisation
 *  written permission.                  préalable et particulière
 *                                       par écrit.
 *
 *  This file is part of the             Ce fichier fait partie du projet
 *  OpenCADC project.                    OpenCADC.
 *
 *  OpenCADC is free software:           OpenCADC est un logiciel libre ;
 *  you can redistribute it and/or       vous pouvez le redistribuer ou le
 *  modify it under the terms of         modifier suivant les termes de
 *  the GNU Affero General Public        la “GNU Affero General Public
 *  License as published by the          License” telle que publiée
 *  Free Software Foundation,            par la Free Software Foundation
 *  either version 3 of the              : soit la version 3 de cette
 *  License, or (at your option)         licence, soit (à votre gré)
 *  any later version.                   toute version ultérieure.
 *
 *  OpenCADC is distributed in the       OpenCADC est distribué
 *  hope that it will be useful,         dans l’espoir qu’il vous
 *  but WITHOUT ANY WARRANTY;            sera utile, mais SANS AUCUNE
 *  without even the implied             GARANTIE : sans même la garantie
 *  warranty of MERCHANTABILITY          implicite de COMMERCIALISABILITÉ
 *  or FITNESS FOR A PARTICULAR          ni d’ADÉQUATION À UN OBJECTIF
 *  PURPOSE.  See the GNU Affero         PARTICULIER. Consultez la Licence
 *  General Public License for           Générale Publique GNU Affero
 *  more details.                        pour plus de détails.
 *
 *  You should have received             Vous devriez avoir reçu une
 *  a copy of the GNU Affero             copie de la Licence Générale
 *  General Public License along         Publique GNU Affero avec
 *  with OpenCADC.  If not, see          OpenCADC ; si ce n’est
 *  .      pas le cas, consultez :
 *                                       .
 *
 *
 ************************************************************************
 */

package org.opencadc.fits.slice;

import ca.nrc.cadc.dali.Circle;
import ca.nrc.cadc.dali.Interval;
import ca.nrc.cadc.dali.PolarizationState;
import ca.nrc.cadc.dali.Polygon;
import ca.nrc.cadc.dali.Range;
import ca.nrc.cadc.dali.Shape;
import ca.nrc.cadc.wcs.exceptions.NoSuchKeywordException;

import java.util.ArrayList;
import java.util.Arrays;
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;


/**
 * Utility class to provide bounds for the given WCS Cutout.  In the unlikely event that multiple cutouts are specified,
 * only one will be honoured, with priority given to POSITION, then BAND, then TIME, and finally POLARIZATION.
 * Pixel cutouts are handled outside of this class.
 */
public class WCSCutoutUtil {
    private static final Logger LOGGER = Logger.getLogger(WCSCutoutUtil.class);

    /**
     * Calculate the pixel ranges that the set of WCS cutouts will produce for the HDU that the given Header represents.
     * @param header        The Header to check.
     * @param cutout        The Cutout specifications.
     * @return  Array of PixelRange objects, or empty if no overlap.  Never null.
     * @throws NoSuchKeywordException Unknown keyword found in Header.
     * @throws HeaderCardException    If a FITS Header card couldn't be read.
     */
    public static PixelRange[] getBounds(final Header header, final Cutout cutout)
            throws HeaderCardException, NoSuchKeywordException {
        final List allPixelRanges = new ArrayList<>();

        if (cutout.pos != null) {
            final PixelRange[] spatialPixelRanges = WCSCutoutUtil.getSpatialBounds(header, cutout.pos);
            if (spatialPixelRanges != null) {
                WCSCutoutUtil.merge(spatialPixelRanges, allPixelRanges);
            }
        }

        if (cutout.band != null) {
            final PixelRange[] spectralPixelRanges = WCSCutoutUtil.getSpectralBounds(header, cutout.band);
            if (spectralPixelRanges != null) {
                WCSCutoutUtil.merge(spectralPixelRanges, allPixelRanges);
            }
        }

        if (cutout.time != null) {
            final PixelRange[] temporalPixelRanges = WCSCutoutUtil.getTemporalBounds(header, cutout.time);
            if (temporalPixelRanges != null) {
                WCSCutoutUtil.merge(temporalPixelRanges, allPixelRanges);
            }
        }

        if (cutout.pol != null) {
            final PixelRange[] polarizationPixelRanges = WCSCutoutUtil.getPolarizationBounds(header, cutout.pol);
            if (polarizationPixelRanges != null) {
                WCSCutoutUtil.merge(polarizationPixelRanges, allPixelRanges);
            }
        }

        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;
        if (shape instanceof Circle) {
            bounds = new CircleCutout(header).getBounds((Circle) shape);
        } else if (shape instanceof Polygon) {
            bounds = new PolygonCutout(header).getBounds((Polygon) shape);
        } else if (shape instanceof Range) {
            bounds = new RangeCutout(header).getBounds((Range) shape);
        } else {
            bounds = null;
        }

        return WCSCutoutUtil.toPixelRanges(bounds);
    }

    static PixelRange[] getSpectralBounds(final Header header, final Interval spectralInterval)
            throws HeaderCardException, NoSuchKeywordException {
        return WCSCutoutUtil.toPixelRanges(new EnergyCutout(header).getBounds(spectralInterval));
    }

    static PixelRange[] getTemporalBounds(final Header header, final Interval temporalInterval)
            throws HeaderCardException {
        return WCSCutoutUtil.toPixelRanges(new TimeCutout(header).getBounds(temporalInterval));
    }

    static PixelRange[] getPolarizationBounds(final Header header, final List polarizationStates)
            throws HeaderCardException {
        return WCSCutoutUtil.toPixelRanges(new PolarizationCutout(header).getBounds(
                polarizationStates.toArray(new PolarizationState[0])));
    }

    /**
     * Merge in the given pixelCutoutRanges array into the entire list of PixelRange objects.  This is to support
     * multiple cutouts in a single HDU (Header), but along different axes.
     * @param pixelCutoutRanges         The NAXIS-length array of PixelRange objects.
     * @param completePixelRangeList    The List of PixelRange objects so far.  Can be empty but not null.
     */
    static void merge(final PixelRange[] pixelCutoutRanges, final List completePixelRangeList) {
        if (completePixelRangeList.isEmpty()) {
            completePixelRangeList.addAll(Arrays.asList(pixelCutoutRanges));
        } else {
            final int pixelCutoutRangesLength = pixelCutoutRanges.length;
            for (int i = 0; i < pixelCutoutRangesLength; i++) {
                final PixelRange pixelCutoutRange = pixelCutoutRanges[i];
                final PixelRange pixelRange = completePixelRangeList.get(i);

                if (!pixelCutoutRange.equals(pixelRange)) {
                    final PixelRange mergedPixelRange =
                            new PixelRange(Math.max(pixelCutoutRange.lowerBound, pixelRange.lowerBound),
                                           Math.min(pixelCutoutRange.upperBound, pixelRange.upperBound));
                    completePixelRangeList.set(i, mergedPixelRange);
                }
            }
        }
    }

    static PixelRange[] toPixelRanges(final long[] pixelBounds) {
        LOGGER.debug("toPixelRanges from bounds " + Arrays.toString(pixelBounds));
        if (pixelBounds == null) {
            return null;
        } else {
            final PixelRange[] pixelRanges = new PixelRange[pixelBounds.length / 2];
            for (int i = 0; i < pixelBounds.length; i += 2) {
                pixelRanges[i / 2] = new PixelRange((int) pixelBounds[i], (int) pixelBounds[i + 1]);
            }
            return pixelRanges;
        }
    }
}




© 2015 - 2024 Weber Informatics LLC | Privacy Policy