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

    friend Wcs::Ptr 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); }

    lsst::afw::coord::Coord::Ptr 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;};

    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;

    // 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;
    
    afw::coord::Coord::Ptr makeCorrectCoord(geom::Angle sky0, geom::Angle sky1) const;

    afw::coord::Coord::Ptr 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;
};

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);
}

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

/*
 Note, CD matrix elements must be in degrees/pixel.
 */
Wcs::Ptr 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;
};  


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

#endif // LSST_AFW_IMAGE_WCS_H

//  LocalWords:  LSST