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:
    typedef std::shared_ptr<Wcs> Ptr;
    typedef std::shared_ptr<Wcs const> ConstPtr;

    friend PTR(Wcs) makeWcs(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 Ptr clone(void) const;

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

    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 PTR(lsst::daf::base::PropertyList) getFitsMetadata() const;

    bool isFlipped() const;

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

    geom::Angle pixelScale() const;

    PTR(coord::Coord) pixelToSky(double pix1, double pix2) const;

    PTR(coord::Coord) pixelToSky(lsst::afw::geom::Point2D const & pixel) const;

    /*
     *  @brief Convert from pixel position to sky coordinates (e.g. RA/dec)
     *
     *  @note This routine is designed for the knowledgeable user in need of performance;
     *  it's safer to call the version that returns a PTR(Coord).
     */
    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(CONST_PTR(lsst::daf::base::PropertySet) 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& operator= (const Wcs &);

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

    PTR(afw::coord::Coord) makeCorrectCoord(geom::Angle sky0, geom::Angle sky1) const;

    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(CONST_PTR(lsst::daf::base::PropertySet) 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 {
    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);
}

PTR(Wcs) makeWcs(PTR(lsst::daf::base::PropertySet) const& fitsMetadata, bool stripMetadata=false);

/*
 Note, CD matrix elements must be in degrees/pixel.
 */
PTR(Wcs) makeWcs(coord::Coord const & crval, geom::Point2D const & crpix,
                 double CD11, double CD12, double CD21, double CD22);

namespace detail {
    int stripWcsKeywords(PTR(lsst::daf::base::PropertySet) const& metadata, 
                         CONST_PTR(Wcs) const& wcs 
                        );
}


class XYTransformFromWcsPair : public afw::geom::XYTransform
{
public:
    XYTransformFromWcsPair(CONST_PTR(Wcs) dst, CONST_PTR(Wcs) src);
    virtual ~XYTransformFromWcsPair() { }

    virtual PTR(afw::geom::XYTransform) invert() const;

    virtual PTR(afw::geom::XYTransform) clone() const;
    virtual Point2D forwardTransform(Point2D const &pixel) const;
    virtual Point2D reverseTransform(Point2D const &pixel) const;

protected:
    CONST_PTR(Wcs) _dst;
    CONST_PTR(Wcs) _src;
    bool const _isSameSkySystem;
};


}}} // lsst::afw::image

#endif // LSST_AFW_IMAGE_WCS_H

//  LocalWords:  LSST