Background.h

Go to the documentation of this file.

// -*- LSST-C++ -*-

/*
 * LSST Data Management System
 * Copyright 2008-2015 AURA/LSST.
 *
 * 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 <https://www.lsstcorp.org/LegalNotices/>.
 */

#if !defined(LSST_AFW_MATH_BACKGROUND_H)
#define LSST_AFW_MATH_BACKGROUND_H
/*
 * Estimate image backgrounds
 */
#include <boost/preprocessor/seq.hpp>
#include <memory>
#include "lsst/daf/base/Citizen.h"
#include "lsst/pex/exceptions.h"
#include "lsst/geom/Box.h"
#include "lsst/afw/math/Statistics.h"
#include "lsst/afw/math/Interpolate.h"
#include "lsst/afw/math/Approximate.h"

namespace lsst {
namespace afw {
namespace math {

//
// 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,
                      ApproximateControl const actrl = ApproximateControl(ApproximateControl::UNKNOWN, 1))
            : _style(Interpolate::AKIMA_SPLINE),
              _nxSample(nxSample),
              _nySample(nySample),
              _undersampleStyle(THROW_EXCEPTION),
              _sctrl(new StatisticsControl(sctrl)),
              _prop(prop),
              _actrl(new ApproximateControl(actrl)) {
        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,
                      ApproximateControl const actrl = ApproximateControl(ApproximateControl::UNKNOWN, 1))
            : _style(Interpolate::AKIMA_SPLINE),
              _nxSample(nxSample),
              _nySample(nySample),
              _undersampleStyle(THROW_EXCEPTION),
              _sctrl(new StatisticsControl(sctrl)),
              _prop(stringToStatisticsProperty(prop)),
              _actrl(new ApproximateControl(actrl)) {
        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,
                      ApproximateControl const actrl = ApproximateControl(ApproximateControl::UNKNOWN, 1)

                              )
            : _style(style),
              _nxSample(nxSample),
              _nySample(nySample),
              _undersampleStyle(undersampleStyle),
              _sctrl(new StatisticsControl(sctrl)),
              _prop(prop),
              _actrl(new ApproximateControl(actrl)) {
        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",
                      ApproximateControl const actrl = ApproximateControl(ApproximateControl::UNKNOWN, 1))
            : _style(math::stringToInterpStyle(style)),
              _nxSample(nxSample),
              _nySample(nySample),
              _undersampleStyle(math::stringToUndersampleStyle(undersampleStyle)),
              _sctrl(new StatisticsControl(sctrl)),
              _prop(stringToStatisticsProperty(prop)),
              _actrl(new ApproximateControl(actrl)) {
        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(BackgroundControl const&) = default;
    BackgroundControl(BackgroundControl&&) = default;
    BackgroundControl& operator=(BackgroundControl const&) = default;
    BackgroundControl& operator=(BackgroundControl&&) = default;

    virtual ~BackgroundControl() = default;
    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; }
    std::shared_ptr<StatisticsControl> getStatisticsControl() { return _sctrl; }
    std::shared_ptr<StatisticsControl const> getStatisticsControl() const { return _sctrl; }

    Property getStatisticsProperty() const { return _prop; }
    void setStatisticsProperty(Property prop) { _prop = prop; }
    void setStatisticsProperty(std::string prop) { _prop = stringToStatisticsProperty(prop); }

    void setApproximateControl(std::shared_ptr<ApproximateControl> actrl) { _actrl = actrl; }
    std::shared_ptr<ApproximateControl> getApproximateControl() { return _actrl; }
    std::shared_ptr<ApproximateControl const> getApproximateControl() const { return _actrl; }

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
    std::shared_ptr<StatisticsControl> _sctrl;   // statistics control object
    Property _prop;                              // statistics Property
    std::shared_ptr<ApproximateControl> _actrl;  // approximate control object
};

class Background : public daf::base::Citizen {
protected:
    template <typename ImageT>
    explicit Background(ImageT const& img, BackgroundControl const& bgCtrl);

    explicit Background(lsst::geom::Box2I const imageBBox, int const nx, int const ny);
    virtual ~Background() = default;

public:
    typedef float InternalPixelT;  

    Background(Background const&) = delete;
    Background(Background&&) = delete;
    Background& operator=(Background const&) = delete;
    Background& operator=(Background&&) = delete;

    virtual Background& operator+=(float const delta) = 0;
    virtual Background& operator-=(float const delta) = 0;
    template <typename PixelT>
    std::shared_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>
    std::shared_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>
    std::shared_ptr<lsst::afw::image::Image<PixelT>> getImage(
            lsst::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>
    std::shared_ptr<lsst::afw::image::Image<PixelT>> getImage(
            lsst::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>
    std::shared_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; }
    std::shared_ptr<math::Approximate<InternalPixelT>> getApproximate(
            ApproximateControl const& actrl,
            UndersampleStyle const undersampleStyle = THROW_EXCEPTION) const {
        InternalPixelT disambiguate = 0;
        return _getApproximate(actrl, undersampleStyle, disambiguate);
    }
    lsst::geom::Box2I getImageBBox() const { return _imgBBox; }

    std::shared_ptr<BackgroundControl> getBackgroundControl() { return _bctrl; }
    std::shared_ptr<BackgroundControl const> getBackgroundControl() const { return _bctrl; }

protected:
    lsst::geom::Box2I _imgBBox;                        
    std::shared_ptr<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
                                 */
// 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 (Background::InternalPixelT)
#define LSST_makeBackground_getApproximate_types (Background::InternalPixelT)
#define LSST_makeBackground_getImage(m, v, T)                                      \
    virtual std::shared_ptr<lsst::afw::image::Image<T>> _getImage(                 \
            lsst::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 std::shared_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)
private:
    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(lsst::geom::Box2I const imageDimensions,
                          image::MaskedImage<InternalPixelT> const& statsImage);

    BackgroundMI(BackgroundMI const&) = delete;
    BackgroundMI(BackgroundMI&&) = delete;
    BackgroundMI& operator=(BackgroundMI const&) = delete;
    BackgroundMI& operator=(BackgroundMI&&) = delete;
    ~BackgroundMI() override = default;

    BackgroundMI& operator+=(float const delta) override;
    BackgroundMI& operator-=(float const delta) override;

    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(LSST_makeBackground_getImage)
    BOOST_PP_SEQ_FOR_EACH(LSST_makeBackground_getImage, override, LSST_makeBackground_getImage_types)
    BOOST_PP_SEQ_FOR_EACH(LSST_makeBackground_getApproximate, override,
                          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>
    std::shared_ptr<image::Image<PixelT>> doGetImage(lsst::geom::Box2I const& bbox,
                                                     Interpolate::Style const interpStyle_,
                                                     UndersampleStyle const undersampleStyle) const;
    // and for _getApproximate
    template <typename PixelT>
    std::shared_ptr<Approximate<PixelT>> doGetApproximate(ApproximateControl const& actrl,
                                                          UndersampleStyle const undersampleStyle) const;
};
template <typename ImageT>
std::shared_ptr<Background> makeBackground(ImageT const& img, BackgroundControl const& bgCtrl) {
    return std::shared_ptr<Background>(new BackgroundMI(img, bgCtrl));
}
}  // namespace math
}  // namespace afw
}  // namespace lsst

#endif  //   LSST_AFW_MATH_BACKGROUND_H