fits_io.h

Go to the documentation of this file.

// -*- lsst-c++ -*-
/*
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */

#ifndef LSST_AFW_IMAGE_fits_io_h_INCLUDED
#define LSST_AFW_IMAGE_fits_io_h_INCLUDED

#include <string>

#include "boost/format.hpp"
#include "lsst/pex/exceptions.h"
#include "lsst/geom.h"
#include "lsst/afw/fits.h"
#include "lsst/afw/geom/wcsUtils.h"
#include "ndarray.h"

namespace lsst {
namespace afw {
namespace image {

template <typename PixelT>
inline void fits_read_array(fits::Fits& fitsfile, ndarray::Array<PixelT, 2, 2>& array,
                            lsst::geom::Point2I& xy0, lsst::daf::base::PropertySet& metadata,
                            lsst::geom::Box2I bbox = lsst::geom::Box2I(), ImageOrigin origin = PARENT) {
    fitsfile.checkCompressedImagePhu();
    if (!fitsfile.checkImageType<PixelT>()) {
        throw LSST_FITS_EXCEPT(fits::FitsTypeError, fitsfile, "Incorrect image type for FITS image");
    }
    int nAxis = fitsfile.getImageDim();
    ndarray::Vector<ndarray::Size, 2> shape;
    if (nAxis == 2) {
        shape = fitsfile.getImageShape<2>();
    } else if (nAxis == 3) {
        ndarray::Vector<ndarray::Size, 3> shape3 = fitsfile.getImageShape<3>();
        if (shape3[0] != 1) {
            throw LSST_EXCEPT(fits::FitsError,
                              boost::str(boost::format("3rd dimension %d is not 1") % shape3[0]));
        }
        shape = shape3.last<2>();
    }
    /*
     * TODO: DM-11235 We cannot call fitsfile.readMetadata(metadata, true)
     * as it doesn't handle INHERIT correctly
     */
    const int hdu = fitsfile.getHdu();
    std::shared_ptr<daf::base::PropertyList> header(readMetadata(fitsfile, true));
    metadata.combine(header);
    fitsfile.setHdu(hdu);

    // Origin of part of image to read
    xy0 = lsst::geom::Point2I();

    lsst::geom::Extent2I xyOffset(geom::getImageXY0FromMetadata(metadata, detail::wcsNameForXY0));
    lsst::geom::Extent2I dimensions = lsst::geom::Extent2I(shape[1], shape[0]);

    if (!bbox.isEmpty()) {
        if (origin == PARENT) {
            bbox.shift(-xyOffset);
        }
        xy0 = bbox.getMin();

        if (bbox.getMinX() < 0 || bbox.getMinY() < 0 || bbox.getWidth() > dimensions.getX() ||
            bbox.getHeight() > dimensions.getY()) {
            throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                              (boost::format("BBox (%d,%d) %dx%d doesn't fit in image %dx%d") %
                               bbox.getMinX() % bbox.getMinY() % bbox.getWidth() % bbox.getHeight() %
                               dimensions.getX() % dimensions.getY())
                                      .str());
        }
        dimensions = bbox.getDimensions();
    }

    array = ndarray::allocate(dimensions.getY(), dimensions.getX());

    fitsfile.readImage(array, ndarray::makeVector(xy0.getY(), xy0.getX()));

    xy0 += xyOffset;
}

}  // namespace image
}  // namespace afw
}  // namespace lsst

#endif  // !LSST_AFW_IMAGE_fits_io_h_INCLUDED