ShapeletKernel.h

Go to the documentation of this file.

// -*- LSST-C++ -*-
#ifndef LSST_MEAS_ALGORITHMS_SHAPELETKERNEL_H
#define LSST_MEAS_ALGORITHMS_SHAPELETKERNEL_H

/* 
 * 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/>.
 */

#include "boost/shared_ptr.hpp"

#include "lsst/afw/math/Kernel.h"
#include "lsst/afw/image/Image.h"
#include "lsst/afw/image/Wcs.h"
#include "lsst/afw/geom/Point.h"
#include "lsst/meas/algorithms/Shapelet.h"
#include "lsst/meas/algorithms/ShapeletInterpolation.h"

namespace lsst {
namespace meas {
namespace algorithms {

    class LocalShapeletKernel : public lsst::afw::math::AnalyticKernel
    {
        /* @brief LocalShapeletKernel has no spatial variation
         *
         * LocalShapeletKernel is appropriate for a small patch of an image 
         * where the variation is expected to be minimal (e.g. over the size 
         * of a galaxy).
         *
         * It uses a shapelet description of the PSF for its underlying
         * implementation.  As such, the functional form is an analytic
         * function f(x,y) = Sum_pq b_pq psi_pq(x,y,sigma).
         * (See the Shapelet class for more details.)
         *
         * The method computeImage is the best way to use this class.
         *
         * The Kernel classes allows for individual pixels to be 
         * calculated one at a time, so we are forced to implement
         * this as well, but it is much less efficient than calculating
         * the whole image at once.
         *
         * There is a suggestion that we also include a method to 
         * calculate the Fourier transform directly.  This would likely
         * be more efficient than going through the normal image
         * via an FFT, since shapelets are their own Fourier transforms
         * (modulo factors of i).  
         * One concern about this, which needs to be considered more
         * carefully, is the effect of the telescope distortion (the Wcs).
         * The shapelets are defined in sky coordinates.  However, the 
         * appropriate image is in chip coordinates.  So there is a
         * local distortion that is applied.  This means that the 
         * direct Fourier transform is not quite as simple as I had
         * initially thought.  However, I do suspect that this detail
         * can be correctly calculated analytically and the method 
         * could be written that is significantly more efficient than
         * the FFT.
         */

    public :
        typedef boost::shared_ptr<LocalShapeletKernel> Ptr;
        typedef boost::shared_ptr<const LocalShapeletKernel> ConstPtr;

        typedef lsst::afw::math::AnalyticKernel base;
        typedef lsst::afw::geom::Point2D Point;
        typedef lsst::afw::geom::Extent2I Extent;
        typedef lsst::afw::image::Image<double> Image;
        typedef lsst::afw::image::Wcs Wcs;

        LocalShapeletKernel(
            Shapelet::ConstPtr shapelet,    
            const Wcs::ConstPtr& wcsPtr,    
            const Extent& size  
        );

        LocalShapeletKernel(
            Shapelet::ConstPtr shapelet, 
            const Wcs::ConstPtr& wcsPtr     
        );

        double computeImage(
            Image& image,       
            bool doNormalize,   
            double x = 0.0,     
            double y = 0.0      
        ) const;

    private :

        Shapelet::ConstPtr _shapelet;
        const Wcs::Ptr _wcsPtr;
    };


    class ShapeletKernel : public lsst::afw::math::AnalyticKernel
    {
        /* @brief ShapeletKernel includes spatial variation
         *
         * A ShapeletKernel is basically a function that can return
         * a LocalShapeletKernel for any location on an Image.
         * This is the most efficient way to use this class.
         * Namely to get the LocalShapeletKernel appropriate for a 
         * particular galaxy, and then convolve with that over a 
         * small patch around the galaxy.  
         *
         * Convolving a larger image is probably not efficient the
         * way that AnalyticKernel seems to be implemented.  AnalyticKernel
         * defines a number of functions that seem to be used for 
         * convolution, like determining the parameters of the 
         * analytic function separately.  The interpolation is faster
         * to calculate the full local function as a complete vector
         * rather than one component at a time.
         */
    public :
        typedef boost::shared_ptr<ShapeletKernel> Ptr;
        typedef boost::shared_ptr<const ShapeletKernel> ConstPtr;

        typedef lsst::afw::math::AnalyticKernel base;
        typedef lsst::afw::geom::Point2D Point;
        typedef lsst::afw::geom::Extent2I Extent;
        typedef lsst::afw::image::Image<double> Image;
        typedef lsst::afw::image::Wcs Wcs;

        ShapeletKernel(
            ShapeletInterpolation::ConstPtr interp,  
            const Wcs::ConstPtr& wcsPtr,    
            const Extent& size  
        );
        ShapeletKernel(
            ShapeletInterpolation::ConstPtr interp,  
            const Wcs::ConstPtr& wcsPtr     
        );

        LocalShapeletKernel::ConstPtr getLocalKernel(
            const Point& pos   
        ) const;

        double computeImage(
            Image& image,       
            bool doNormalize,   
            double x = 0.0,     
            double y = 0.0      
        ) const;

    private :

        ShapeletInterpolation::ConstPtr _interp;
        const Wcs::Ptr _wcsPtr;
    };

}}}


#endif