Background.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/>.
 */
 
#if !defined(LSST_AFW_MATH_BACKGROUND_H)
#define LSST_AFW_MATH_BACKGROUND_H

#include <boost/preprocessor/seq.hpp>
#include "boost/shared_ptr.hpp"
#include "lsst/daf/base/Citizen.h"
#include "lsst/pex/exceptions.h"
#include "lsst/afw/geom/Box.h"
#include "lsst/afw/math/Statistics.h"
#include "lsst/afw/math/Interpolate.h"

namespace lsst {
namespace afw {
namespace math {
class ApproximateControl;
template<typename T> class Approximate;

//
// Remember to update stringToUndersampleStyle if you change this.
// If this happens often, we can play CPP games to put the definition in exactly one place, although swig
// may not be happy (so we could think m4 thoughts instead)
//
enum UndersampleStyle {
    THROW_EXCEPTION,
    REDUCE_INTERP_ORDER,
    INCREASE_NXNYSAMPLE
};
UndersampleStyle stringToUndersampleStyle(std::string const &style);
    
class BackgroundControl {
public:
    BackgroundControl(
        int const nxSample,                                  
        int const nySample,                                  
        StatisticsControl const sctrl = StatisticsControl(), 
        Property const prop = MEANCLIP 
                     )
        : _style(Interpolate::AKIMA_SPLINE),
          _nxSample(nxSample), _nySample(nySample),
          _undersampleStyle(THROW_EXCEPTION),
          _sctrl(new StatisticsControl(sctrl)),
          _prop(prop) {
        if (nxSample <= 0 || nySample <= 0) {
            throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                              str(boost::format("You must specify at least one point, not %dx%d")
                                  % nxSample % nySample)
                             );
        }
    }
    
    BackgroundControl(
        int const nxSample,             
        int const nySample,             
        StatisticsControl const &sctrl, 
        std::string const &prop         
                     )
        : _style(Interpolate::AKIMA_SPLINE),
          _nxSample(nxSample), _nySample(nySample),
          _undersampleStyle(THROW_EXCEPTION),
          _sctrl(new StatisticsControl(sctrl)),
          _prop(stringToStatisticsProperty(prop)) {
        if (nxSample <= 0 || nySample <= 0) {
            throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                              str(boost::format("You must specify at least one point, not %dx%d")
                                  % nxSample % nySample)
                             );
        }
    }
    // And now the two old APIs (preserved for backward compatibility)
    BackgroundControl(
        Interpolate::Style const style, 
        int const nxSample = 10,        
        int const nySample = 10,        
        UndersampleStyle const undersampleStyle = THROW_EXCEPTION, 
        StatisticsControl const sctrl = StatisticsControl(), 
        Property const prop = MEANCLIP 
                     )
        : _style(style),
          _nxSample(nxSample), _nySample(nySample),
          _undersampleStyle(undersampleStyle),
          _sctrl(new StatisticsControl(sctrl)),
          _prop(prop) {
        if (nxSample <= 0 || nySample <= 0) {
            throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                              str(boost::format("You must specify at least one point, not %dx%d")
                                  % nxSample % nySample)
                             );
        }
    }
    
    BackgroundControl(
        std::string const &style, 
        int const nxSample = 10, 
        int const nySample = 10, 
        std::string const &undersampleStyle = "THROW_EXCEPTION", 
        StatisticsControl const sctrl = StatisticsControl(), 
        std::string const &prop = "MEANCLIP" 
                     )
        : _style(math::stringToInterpStyle(style)),
          _nxSample(nxSample), _nySample(nySample),
          _undersampleStyle(math::stringToUndersampleStyle(undersampleStyle)),
          _sctrl(new StatisticsControl(sctrl)),
          _prop(stringToStatisticsProperty(prop)) {
        if (nxSample <= 0 || nySample <= 0) {
            throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                              str(boost::format("You must specify at least one point, not %dx%d")
                                  % nxSample % nySample)
                             );
        }
    }

    virtual ~BackgroundControl() {}
    void setNxSample (int nxSample) {
        if (nxSample <= 0) {
            throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                              str(boost::format("nxSample must be position, not %d") % nxSample));
        }
        _nxSample = nxSample;
    }
    void setNySample (int nySample) {
        if (nySample <= 0) {
            throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                              str(boost::format("nySample must be position, not %d") % nySample));
        }
        _nySample = nySample;
    }

    void setInterpStyle (Interpolate::Style const style) { _style = style; }
    // overload to take a string
    void setInterpStyle (std::string const &style) { _style = math::stringToInterpStyle(style); }
    
    void setUndersampleStyle (UndersampleStyle const undersampleStyle) {
        _undersampleStyle = undersampleStyle;
    }
    // overload to take a string
    void setUndersampleStyle (std::string const &undersampleStyle) {
        _undersampleStyle = math::stringToUndersampleStyle(undersampleStyle);
    }
    
    int getNxSample() const { return _nxSample; }
    int getNySample() const { return _nySample; }
    Interpolate::Style getInterpStyle() const {
        if (_style < 0 || _style >= Interpolate::NUM_STYLES) {
            throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                              str(boost::format("Style %d is invalid") % _style));

        }
        return _style;
    }
    UndersampleStyle getUndersampleStyle() const { return _undersampleStyle; }
    PTR(StatisticsControl) getStatisticsControl() { return _sctrl; }
    CONST_PTR(StatisticsControl) getStatisticsControl() const { return _sctrl; }

    Property getStatisticsProperty() const { return _prop; }
    void setStatisticsProperty(Property prop) { _prop = prop; }
    void setStatisticsProperty(std::string prop) { _prop = stringToStatisticsProperty(prop); }
    
private:
    Interpolate::Style _style;          // style of interpolation to use
    int _nxSample;                      // number of grid squares to divide image into to sample in x
    int _nySample;                      // number of grid squares to divide image into to sample in y
    UndersampleStyle _undersampleStyle; // what to do when nx,ny are too small for the requested interp style
    PTR(StatisticsControl) _sctrl;           // statistics control object
    Property _prop;                          // statistics Property
};
    
class Background : public daf::base::Citizen {
protected:
    template<typename ImageT>
    explicit Background(ImageT const& img, BackgroundControl const& bgCtrl);

    explicit Background(geom::Box2I const imageBBox, int const nx, int const ny);
    virtual ~Background() { }
public:
    typedef float InternalPixelT;    

    virtual void operator+=(float const delta) = 0;
    virtual void operator-=(float const delta) = 0;
    template<typename PixelT>
    PTR(lsst::afw::image::Image<PixelT>) getImage(
        Interpolate::Style const interpStyle,                   
        UndersampleStyle const undersampleStyle=THROW_EXCEPTION 
                                                 ) const {
        return getImage<PixelT>(_imgBBox, interpStyle, undersampleStyle);
    }
    template<typename PixelT>
    PTR(lsst::afw::image::Image<PixelT>) getImage(
        std::string const &interpStyle, 
        std::string const &undersampleStyle="THROW_EXCEPTION"   
                                                 ) const {
        return getImage<PixelT>(math::stringToInterpStyle(interpStyle),
                                stringToUndersampleStyle(undersampleStyle));
    }
    template<typename PixelT>
    PTR(lsst::afw::image::Image<PixelT>) getImage(
        lsst::afw::geom::Box2I const& bbox,                     
        Interpolate::Style const interpStyle,                   
        UndersampleStyle const undersampleStyle=THROW_EXCEPTION 
        ) const {
        return _getImage(bbox, interpStyle, undersampleStyle, static_cast<PixelT>(0));
    }
    template<typename PixelT>
    PTR(lsst::afw::image::Image<PixelT>) getImage(
        lsst::afw::geom::Box2I const& bbox,               
        std::string const& interpStyle,                   
        std::string const& undersampleStyle="THROW_EXCEPTION" 
        ) const {
        return _getImage(bbox, math::stringToInterpStyle(interpStyle),
                         stringToUndersampleStyle(undersampleStyle), static_cast<PixelT>(0));
    }

    template<typename PixelT>
    PTR(lsst::afw::image::Image<PixelT>) getImage() const {
        return getImage<PixelT>(_bctrl.getInterpStyle(), _bctrl.getUndersampleStyle());
    }
    Interpolate::Style getAsUsedInterpStyle() const {
        return _asUsedInterpStyle;
    }
    UndersampleStyle getAsUsedUndersampleStyle() const {
        return _asUsedUndersampleStyle;
    }
    PTR(math::Approximate<InternalPixelT>) getApproximate(
        ApproximateControl const& actrl,                        
        UndersampleStyle const undersampleStyle=THROW_EXCEPTION 
                                                 ) const {
        InternalPixelT disambiguate = 0;
        return _getApproximate(actrl, undersampleStyle, disambiguate);
    }
    geom::Box2I getImageBBox() const { return _imgBBox; }

protected:
    geom::Box2I _imgBBox;                             
    BackgroundControl _bctrl;                         
    mutable Interpolate::Style _asUsedInterpStyle;    
    mutable UndersampleStyle _asUsedUndersampleStyle; 

    std::vector<double> _xcen;          
    std::vector<double> _ycen;          
    std::vector<int> _xorig;            
    std::vector<int> _yorig;            
    std::vector<int> _xsize;            
    std::vector<int> _ysize;            
    /*
     * We want getImage to be present in the base class, but a templated virtual function
     * is impossible.  So we'll solve the dilemma with a hack: explicitly defined
     * virtual functions for the image types we need
     */
#if !defined(SWIG)
// We'll evaluate LSST_makeBackground_get{Approximation,Image} for each type in
// LSST_makeBackground_get{Approximation,Image}_types,
// setting v to the second arg (i.e. "= 0" for the first invocation).  The first agument, m, is ignores

// Desired types
#define LSST_makeBackground_getImage_types            (double)(float)(int)
#define LSST_makeBackground_getApproximate_types      (Background::InternalPixelT)
#define LSST_makeBackground_getImage(m, v, T)                \
    virtual PTR(lsst::afw::image::Image<T>) _getImage( \
        lsst::afw::geom::Box2I const& bbox, \
        Interpolate::Style const interpStyle,                     /* Style of the interpolation */ \
        UndersampleStyle const undersampleStyle=THROW_EXCEPTION,  /* Behaviour if there are too few points */\
        T = 0                                                     /* disambiguate */ \
                                                     ) const v;

#define LSST_makeBackground_getApproximate(m, v, T)           \
    virtual PTR(Approximate<T>) _getApproximate( \
        ApproximateControl const& actrl,                          /* Approximation style */ \
        UndersampleStyle const undersampleStyle=THROW_EXCEPTION,  /* Behaviour if there are too few points */\
        T = 0                                                     /* disambiguate */ \
                                                      ) const v;

    BOOST_PP_SEQ_FOR_EACH(LSST_makeBackground_getImage, = 0, LSST_makeBackground_getImage_types)
    BOOST_PP_SEQ_FOR_EACH(LSST_makeBackground_getApproximate, = 0, LSST_makeBackground_getApproximate_types)
#endif
private:
    Background(Background const&);
    Background& operator=(Background const&);    
    void _setCenOrigSize(int const width, int const height, int const nxSample, int const nySample);
};
    
class BackgroundMI : public Background {
public:
    template<typename ImageT>
    explicit BackgroundMI(ImageT const& img,
                        BackgroundControl const& bgCtrl);
    explicit BackgroundMI(geom::Box2I const imageDimensions,
                          image::MaskedImage<InternalPixelT> const& statsImage);
    
    virtual void operator+=(float const delta);
    virtual void operator-=(float const delta);


    double getPixel(Interpolate::Style const style, int const x, int const y) const;
    double getPixel(int const x, int const y) const {
        return getPixel(_bctrl.getInterpStyle(), x, y);
    }
    lsst::afw::image::MaskedImage<InternalPixelT> getStatsImage() const {
        return _statsImage;
    }

private:
    lsst::afw::image::MaskedImage<InternalPixelT> _statsImage; // statistical properties for the grid of subimages
    mutable std::vector<std::vector<double> > _gridColumns; // interpolated columns for the bicubic spline

    void _setGridColumns(Interpolate::Style const interpStyle,
                         UndersampleStyle const undersampleStyle,
                         int const iX, std::vector<int> const& ypix) const;

#if !defined(SWIG) && defined(LSST_makeBackground_getImage)
    BOOST_PP_SEQ_FOR_EACH(LSST_makeBackground_getImage, , LSST_makeBackground_getImage_types)
    BOOST_PP_SEQ_FOR_EACH(LSST_makeBackground_getApproximate, , LSST_makeBackground_getApproximate_types)
#if 0                                   // keep for use in Background instantiations
#undef LSST_makeBackground_getImage_types
#undef LSST_makeBackground_getApproximate_types
#endif
#undef LSST_makeBackground_getImage
#undef LSST_makeBackground_getApproximate
#endif
    // Here's the worker function for _getImage (non-virtual; it's templated in BackgroundMI, not Background)
    template<typename PixelT>
    PTR(image::Image<PixelT>) doGetImage(geom::Box2I const& bbox,
                                         Interpolate::Style const interpStyle_,
                                         UndersampleStyle const undersampleStyle) const;
    // and for _getApproximate
    template<typename PixelT>
    PTR(Approximate<PixelT>) doGetApproximate(ApproximateControl const& actrl,
                                              UndersampleStyle const undersampleStyle) const;
};
template<typename ImageT>
PTR(Background) makeBackground(ImageT const& img, BackgroundControl const& bgCtrl) {
    return PTR(Background)(new BackgroundMI(img, bgCtrl));
}
    
}}}

#endif  //   LSST_AFW_MATH_BACKGROUND_H