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 <memory>
#include <vector>
 
#include "lsst/base.h"
#include "Eigen/Core"
 
#include "lsst/afw/table/Match.h"
#include "lsst/geom/Angle.h"
#include "lsst/geom/Box.h"
 
namespace lsst {
namespace meas {
namespace astrom {
namespace sip {

template <class MatchT>
class CreateWcsWithSip {
public:
    typedef std::shared_ptr<CreateWcsWithSip> Ptr;
    typedef std::shared_ptr<CreateWcsWithSip const> ConstPtr;

    CreateWcsWithSip(std::vector<MatchT> const& matches, afw::geom::SkyWcs const& linearWcs, int const order,
                     geom::Box2I const& bbox = geom::Box2I(), int const ngrid = 0);
 
    std::shared_ptr<afw::geom::SkyWcs> getNewWcs() { return _newWcs; }

    double getScatterInPixels() const;

    geom::Angle getScatterOnSky() const;

    double getLinearScatterInPixels() const;

    geom::Angle getLinearScatterOnSky() const;

    int getOrder() const { return _sipA.rows(); }
    int getNPoints() const { return _matches.size(); }
    int getNGrid() const { return _ngrid; }
 
    // Return the SIP A matrix
    Eigen::MatrixXd const getSipA() { return _sipA; }
    // Return the SIP B matrix
    Eigen::MatrixXd const getSipB() { return _sipB; }
    // Return the SIP Ap matrix
    Eigen::MatrixXd const getSipAp() { return _sipAp; }
    // Return the SIP Bp matrix
    Eigen::MatrixXd const getSipBp() { return _sipBp; }
 
private:
    std::vector<MatchT> const _matches;
    geom::Box2I mutable _bbox;
    int _ngrid;  // grid size to calculate inverse SIP coefficients (1-D)
    std::shared_ptr<const afw::geom::SkyWcs> _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;
 
    std::shared_ptr<afw::geom::SkyWcs> _newWcs;
 
    void _calculateForwardMatrices();
    void _calculateReverseMatrices();
};

template <class MatchT>
CreateWcsWithSip<MatchT> makeCreateWcsWithSip(std::vector<MatchT> const& matches,
                                              afw::geom::SkyWcs const& linearWcs, int const order,
                                              geom::Box2I const& bbox = geom::Box2I(), int const ngrid = 0) {
    return CreateWcsWithSip<MatchT>(matches, linearWcs, order, bbox, ngrid);
}
 
}  // namespace sip
}  // namespace astrom
}  // namespace meas
}  // namespace lsst
 
#endif