CreateWcsWithSip.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 CREATE_WCS_WITH_SIP
#define CREATE_WCS_WITH_SIP

#include <vector>

#include "boost/shared_ptr.hpp"
#include "lsst/base.h"
#include "Eigen/Core"

#include "lsst/afw/table/Match.h"
#include "lsst/afw/geom/Box.h"
#include "lsst/afw/geom/Point.h"
#include "lsst/afw/geom/Angle.h"
#include "lsst/pex/logging/Log.h"

namespace lsst { 
    namespace afw {
        namespace image {
            class Wcs;
            class TanWcs;
        }
    }

namespace meas { 
namespace astrom { 
namespace sip {


template<class MatchT>
class CreateWcsWithSip {
public:

    typedef boost::shared_ptr<CreateWcsWithSip> Ptr;
    typedef boost::shared_ptr<CreateWcsWithSip const> ConstPtr;

    CreateWcsWithSip(
        std::vector<MatchT> const & matches,
        afw::image::Wcs const & linearWcs,
        int const order,
        afw::geom::Box2I const& bbox = afw::geom::Box2I(),
        int const ngrid=0
    );

    PTR(afw::image::TanWcs) getNewWcs() { return _newWcs; }

    /*
     Returns the median separation between points in this object's match list,
     projecting reference sources from RA,Dec to pixels using the SIP WCS, and
     comparing with the matched source pixel positions.
     */
    double getScatterInPixels();

    /*
     Returns the median separation between points in this object's match list,
     projecting sources from pixel space to RA,Dec using the SIP WCS, and
     comparing with the reference source RA,Dec positions.
     */
    afw::geom::Angle getScatterOnSky();

    /*
     Returns the median separation between points in this object's match list,
     projecting reference sources from RA,Dec to pixels using the input TAN
     (linear) WCS, and comparing with the matched source pixel positions.
     */
    double getLinearScatterInPixels();

    /*
     Returns the median separation between points in this object's match list,
     projecting sources from pixel space to RA,Dec using the input TAN (linear)
     WCS, and comparing with the reference source RA,Dec positions.
     */
    afw::geom::Angle getLinearScatterOnSky();

    int getOrder() const { return  _sipA.rows(); }
    int getNPoints() const { return _matches.size(); }
    int getNGrid() const { return _ngrid; }

private:

    lsst::pex::logging::Log _log;
    
    std::vector<MatchT> const _matches;
    afw::geom::Box2I mutable _bbox;
    int _ngrid;                         // grid size to calculate inverse SIP coefficients (1-D)
    CONST_PTR(afw::image::Wcs) _linearWcs;
    // _sipOrder is polynomial order for forward transform.
    // _reverseSipOrder is order for reverse transform, not necessarily the same.
    int const _sipOrder, _reverseSipOrder;      

    Eigen::MatrixXd _sipA, _sipB;
    Eigen::MatrixXd _sipAp, _sipBp;

    PTR(afw::image::TanWcs) _newWcs;

    double _getScatterPixels(afw::image::Wcs const& wcs,
                             std::vector<MatchT> const & matches);
    afw::geom::Angle _getScatterSky(afw::image::Wcs const& wcs,
                                    std::vector<MatchT> const & matches);
    
    void _calculateForwardMatrices();
    void _calculateReverseMatrices();
    
    afw::geom::Point2D _getCrvalAsGeomPoint();
};    

template<class MatchT>
CreateWcsWithSip<MatchT> makeCreateWcsWithSip(
    std::vector<MatchT> const & matches,
    afw::image::Wcs const& linearWcs,
    int const order,
    afw::geom::Box2I const& bbox = afw::geom::Box2I(),
    int const ngrid=0
) {
    return CreateWcsWithSip<MatchT>(matches, linearWcs, order, bbox, ngrid);
}

}}}}

#endif