Wcs.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_WCS_H
#define LSST_AFW_IMAGE_WCS_H

#include <limits>
#include "Eigen/Core"
#include "lsst/base.h"
#include "lsst/daf/base/Persistable.h"
#include "lsst/daf/base/Citizen.h"
#include "lsst/afw/image/Image.h"
#include "lsst/afw/coord/Coord.h"
#include "lsst/afw/geom/AffineTransform.h"
#include "lsst/afw/geom/Point.h"
#include "lsst/afw/geom/Extent.h"
#include "lsst/afw/geom/XYTransform.h"
#include "lsst/afw/table/io/Persistable.h"

struct wcsprm;  // defined in wcs.h

namespace lsst {
namespace daf {
namespace base {
class PropertySet;
}
}
namespace afw {
namespace formatters {
class WcsFormatter;
}
namespace table {
class BaseRecord;
}
namespace image {


class Wcs : public lsst::daf::base::Persistable,
            public lsst::daf::base::Citizen,
            public afw::table::io::PersistableFacade<Wcs>,
            public afw::table::io::Persistable {
public:
    friend std::shared_ptr<Wcs> makeWcs(std::shared_ptr<lsst::daf::base::PropertySet> const& fitsMetadata,
                                        bool stripMetadata);

    Wcs(lsst::afw::geom::Point2D const& crval, lsst::afw::geom::Point2D const& crpix,
        Eigen::Matrix2d const& CD, std::string const& ctype1 = "RA---TAN",
        std::string const& ctype2 = "DEC--TAN", double equinox = 2000, std::string const& raDecSys = "ICRS",
        std::string const& cunits1 = "deg", std::string const& cunits2 = "deg");

    virtual ~Wcs();
    virtual std::shared_ptr<Wcs> clone(void) const;

    bool operator==(Wcs const& other) const;
    bool operator!=(Wcs const& other) const { return !(*this == other); }

    std::shared_ptr<lsst::afw::coord::Coord> getSkyOrigin() const;

    lsst::afw::geom::Point2D getPixelOrigin() const;

    Eigen::Matrix2d getCDMatrix() const;

    virtual void flipImage(int flipLR, int flipTB, lsst::afw::geom::Extent2I dimensions) const;

    virtual void rotateImageBy90(int nQuarter, lsst::afw::geom::Extent2I dimensions) const;

    virtual std::shared_ptr<lsst::daf::base::PropertyList> getFitsMetadata() const;

    bool isFlipped() const;

    double pixArea(lsst::afw::geom::Point2D pix00) const;

    geom::Angle pixelScale() const;

    std::shared_ptr<coord::Coord> pixelToSky(double pix1, double pix2) const;

    std::shared_ptr<coord::Coord> pixelToSky(lsst::afw::geom::Point2D const& pixel) const;

    void pixelToSky(double pixel1, double pixel2, geom::Angle& sky1, geom::Angle& sky2) const;

    geom::Point2D skyToPixel(geom::Angle sky1, geom::Angle sky2) const;

    geom::Point2D skyToPixel(coord::Coord const& coord) const;

    geom::Point2D skyToIntermediateWorldCoord(coord::Coord const& coord) const;

    virtual bool hasDistortion() const { return false; };

    afw::coord::CoordSystem getCoordSystem() const { return _coordSystem; };

    double getEquinox() const;

    bool isSameSkySystem(Wcs const& wcs) const;

    geom::LinearTransform getLinearTransform() const;

    geom::AffineTransform linearizePixelToSky(coord::Coord const& coord,
                                              geom::AngleUnit skyUnit = geom::degrees) const;

    geom::AffineTransform linearizePixelToSky(geom::Point2D const& pix,
                                              geom::AngleUnit skyUnit = geom::degrees) const;

    geom::AffineTransform linearizeSkyToPixel(coord::Coord const& coord,
                                              geom::AngleUnit skyUnit = geom::degrees) const;

    geom::AffineTransform linearizeSkyToPixel(geom::Point2D const& pix,
                                              geom::AngleUnit skyUnit = geom::degrees) const;

    // Mutators; the first one is virtual, even though it will never be overridden,
    // to make sure subclasses use the correct version of both

    virtual void shiftReferencePixel(double dx, double dy);

    // Virtual to make sure subclasses use the correct version of both shiftReferencePixel mutators.
    virtual void shiftReferencePixel(geom::Extent2D const& d) { shiftReferencePixel(d.getX(), d.getY()); }

    virtual bool isPersistable() const;

private:
    // Allow the formatter to access private goo
    LSST_PERSIST_FORMATTER(lsst::afw::formatters::WcsFormatter)

    
    void initWcsLib(geom::Point2D const& crval, geom::Point2D const& crpix, Eigen::Matrix2d const& CD,
                    std::string const& ctype1, std::string const& ctype2, double equinox,
                    std::string const& raDecSys, std::string const& cunits1, std::string const& cunits2);

protected:
    friend class WcsFactory;
    bool _mayBePersistable() const;
    // See afw::table::io::Persistable
    virtual std::string getPersistenceName() const;
    virtual std::string getPythonModule() const;
    virtual void write(OutputArchiveHandle& handle) const;

    // Protected virtual implementation for operator== (must be true in both directions for equality).
    virtual bool _isSubset(Wcs const& other) const;

    // Return true if coordinate system is ICRS or FK5 J2000
    bool _isIcrs() const {
        return (getCoordSystem() == afw::coord::ICRS) ||
               ((getCoordSystem() == afw::coord::FK5) && (getEquinox() == 2000));
    }

    // Default constructor, only used by WcsFormatter
    Wcs();

    // If you want to create a Wcs from a FITS header, use makeWcs().
    // This is protected because the derived classes need to be able to see it.
    Wcs(std::shared_ptr<daf::base::PropertySet const> const& fitsMetadata);

    // Construct from a record; used by WcsFactory for afw::table::io persistence.
    explicit Wcs(afw::table::BaseRecord const& record);

    Wcs(Wcs const& rhs);
    Wcs(Wcs&& rhs);
    Wcs& operator=(const Wcs&);
    Wcs& operator=(Wcs&&);

    virtual void pixelToSkyImpl(double pixel1, double pixel2, geom::Angle skyTmp[2]) const;
    virtual geom::Point2D skyToPixelImpl(geom::Angle sky1, geom::Angle sky2) const;

    std::shared_ptr<afw::coord::Coord> makeCorrectCoord(geom::Angle sky0, geom::Angle sky1) const;

    std::shared_ptr<afw::coord::Coord> convertCoordToSky(coord::Coord const& coord) const;

    virtual geom::AffineTransform linearizePixelToSkyInternal(geom::Point2D const& pix,
                                                              coord::Coord const& coord,
                                                              geom::AngleUnit skyUnit) const;

    virtual geom::AffineTransform linearizeSkyToPixelInternal(geom::Point2D const& pix,
                                                              coord::Coord const& coord,
                                                              geom::AngleUnit skyUnit) const;

    void initWcsLibFromFits(std::shared_ptr<lsst::daf::base::PropertySet const> const& fitsMetadata);
    void _initWcs();
    void _setWcslibParams();

    struct wcsprm* _wcsInfo;
    int _nWcsInfo;
    int _relax;       
    int _wcsfixCtrl;  
    int _wcshdrCtrl;  
    int _nReject;
    coord::CoordSystem _coordSystem;
    bool _skyAxesSwapped;  
};

namespace detail {
std::shared_ptr<lsst::daf::base::PropertyList> createTrivialWcsAsPropertySet(std::string const& wcsName,
                                                                             int const x0 = 0,
                                                                             int const y0 = 0);

geom::Point2I getImageXY0FromMetadata(std::string const& wcsName, lsst::daf::base::PropertySet* metadata);
}

std::shared_ptr<Wcs> makeWcs(std::shared_ptr<lsst::daf::base::PropertySet> const& fitsMetadata,
                             bool stripMetadata = false);

std::shared_ptr<Wcs> makeWcs(coord::Coord const& crval, geom::Point2D const& crpix, double CD11, double CD12,
                             double CD21, double CD22);

namespace detail {
int stripWcsKeywords(
        std::shared_ptr<lsst::daf::base::PropertySet> const& metadata,  
        std::shared_ptr<Wcs const> const& wcs                           
        );
}

class XYTransformFromWcsPair : public afw::geom::XYTransform {
public:
    XYTransformFromWcsPair(std::shared_ptr<Wcs const> dst, std::shared_ptr<Wcs const> src);
    virtual ~XYTransformFromWcsPair() {}

    virtual std::shared_ptr<afw::geom::XYTransform> invert() const;

    virtual std::shared_ptr<afw::geom::XYTransform> clone() const;
    virtual Point2D forwardTransform(Point2D const& pixel) const;
    virtual Point2D reverseTransform(Point2D const& pixel) const;

protected:
    std::shared_ptr<Wcs const> _dst;
    std::shared_ptr<Wcs const> _src;
    bool const _isSameSkySystem;
};
}
}
}  // lsst::afw::image

#endif  // LSST_AFW_IMAGE_WCS_H

//  LocalWords:  LSST