SpanSet.h

Go to the documentation of this file.

// -*- lsst-c++ -*-

/*
 * LSST Data Management System
 * Copyright 2008-2016  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/>.
 */

#ifndef LSST_AFW_GEOM_SPANSET_H
#define LSST_AFW_GEOM_SPANSET_H

#include <vector>
#include <algorithm>
#include <functional>
#include <memory>
#include <utility>
#include "lsst/pex/exceptions.h"
#include "lsst/afw/geom/Span.h"
#include "lsst/geom/Box.h"
#include "lsst/afw/image/Mask.h"
#include "lsst/afw/table/io/Persistable.h"
#include "lsst/afw/geom/ellipses/Ellipse.h"
#include "lsst/afw/geom/ellipses/Quadrupole.h"
#include "lsst/afw/geom/SpanSetFunctorGetters.h"
#include "lsst/afw/geom/Transform.h"
#include "lsst/afw/image/Image.h"
#include "lsst/afw/image/MaskedImage.h"

namespace lsst {
namespace afw {
namespace geom {
namespace details {

/* Functor object to be used with fromMask function
 */
template <typename T>
class AnyBitSetFunctor {
public:
    bool operator()(T pixelValue) { return pixelValue != 0; }
};

}  // namespace details

enum class Stencil { CIRCLE, BOX, MANHATTAN };

class SpanSet : public afw::table::io::PersistableFacade<lsst::afw::geom::SpanSet>,
                public afw::table::io::Persistable {
public:
    typedef std::vector<Span>::const_iterator const_iterator;
    typedef std::vector<Span>::size_type size_type;
    typedef Span value_type;
    typedef value_type const &const_reference;

    // Expose properties of the underlying vector containing spans such that the
    // SpanSet can be considered a container.
    // Return the constant versions as SpanSets should be immutable
    const_iterator begin() const { return _spanVector.cbegin(); }
    const_iterator end() const { return _spanVector.cend(); }
    const_iterator cbegin() const { return _spanVector.cbegin(); }
    const_iterator cend() const { return _spanVector.cend(); }
    const_reference front() const { return const_cast<geom::Span &>(_spanVector.front()); }
    const_reference back() const { return const_cast<geom::Span &>(_spanVector.back()); }
    size_type size() const { return _spanVector.size(); }
    bool empty() const { return _spanVector.empty(); }

    SpanSet();

    explicit SpanSet(lsst::geom::Box2I const &box);

    template <typename iter>
    SpanSet(iter begin, iter end, bool normalize = true) : _spanVector(begin, end) {
        // Return a null SpanSet if spanVector is 0
        if (_spanVector.size() == 0) {
            _bbox = lsst::geom::Box2I();
            _area = 0;
        } else {
            if (normalize) {
                _runNormalize();
            }
            _initialize();
        }
    }

    explicit SpanSet(std::vector<Span> const &vec, bool normalize = true);

    explicit SpanSet(std::vector<Span> &&vec, bool normalize = true);

    // Explicitly delete copy and move constructors
    SpanSet(SpanSet const &other) = delete;
    SpanSet(SpanSet &&other) = delete;
    ~SpanSet() override = default;

    SpanSet &operator=(SpanSet const &) = default;
    // Delegate to copy-assignment for backwards compatibility
    SpanSet &operator=(SpanSet &&other) { return *this = other; }

    // Define class methods
    size_type getArea() const;

    lsst::geom::Box2I getBBox() const;

    bool isContiguous() const;

    std::shared_ptr<SpanSet> shiftedBy(int x, int y) const;

    std::shared_ptr<SpanSet> shiftedBy(lsst::geom::Extent2I const &offset) const;

    std::shared_ptr<SpanSet> clippedTo(lsst::geom::Box2I const &box) const;

    std::shared_ptr<SpanSet> transformedBy(lsst::geom::LinearTransform const &t) const;

    std::shared_ptr<SpanSet> transformedBy(lsst::geom::AffineTransform const &t) const;

    std::shared_ptr<SpanSet> transformedBy(TransformPoint2ToPoint2 const &t) const;

    bool overlaps(SpanSet const &other) const;

    bool contains(SpanSet const &other) const;

    bool contains(lsst::geom::Point2I const &point) const;

    lsst::geom::Point2D computeCentroid() const;

    ellipses::Quadrupole computeShape() const;

    std::shared_ptr<SpanSet> dilated(int r, Stencil s = Stencil::CIRCLE) const;

    std::shared_ptr<SpanSet> dilated(SpanSet const &other) const;

    std::shared_ptr<SpanSet> eroded(int r, Stencil s = Stencil::CIRCLE) const;

    std::shared_ptr<SpanSet> eroded(SpanSet const &other) const;

    template <typename Pixel, int inN, int inC>
    ndarray::Array<typename std::remove_const<Pixel>::type, inN - 1, inN - 1> flatten(
            ndarray::Array<Pixel, inN, inC> const &input,
            lsst::geom::Point2I const &xy0 = lsst::geom::Point2I()) const {
        // Populate a lower dimensional array with the values from input taken at the points of SpanSet
        auto outputShape = ndarray::concatenate(ndarray::makeVector(getArea()),
                                                input.getShape().template last<inN - 2>());
        ndarray::Array<typename std::remove_const<Pixel>::type, inN - 1, inN - 1> outputArray =
                ndarray::allocate(outputShape);
        outputArray.deep() = 0;
        flatten(outputArray, input, xy0);
        return outputArray;
    }

    template <typename PixelIn, typename PixelOut, int inA, int outC, int inC>
    void flatten(ndarray::Array<PixelOut, inA - 1, outC> const &output,
                 ndarray::Array<PixelIn, inA, inC> const &input,
                 lsst::geom::Point2I const &xy0 = lsst::geom::Point2I()) const {
        auto ndAssigner = [](lsst::geom::Point2I const &point,
                             typename details::FlatNdGetter<PixelOut, inA - 1, outC>::Reference out,
                             typename details::ImageNdGetter<PixelIn, inA, inC>::Reference in) { out = in; };
        // Populate array output with values from input at positions given by SpanSet
        applyFunctor(ndAssigner, ndarray::ndFlat(output), ndarray::ndImage(input, xy0));
    }

    template <typename Pixel, int inA, int inC>
    ndarray::Array<typename std::remove_const<Pixel>::type, inA + 1, inA + 1> unflatten(
            ndarray::Array<Pixel, inA, inC> const &input) const {
        // Create a higher dimensional array the size of the bounding box and extra dimensions of input.
        // Populate values from input, placed at locations corresponding to SpanSet, offset by the
        // lower corner of the bounding box
        auto existingShape = input.getShape();
        typename decltype(existingShape)::Element height = _bbox.getHeight();
        typename decltype(existingShape)::Element width = _bbox.getWidth();
        auto outputShape = ndarray::concatenate(ndarray::makeVector(height, width),
                                                input.getShape().template last<inA - 1>());
        ndarray::Array<typename std::remove_const<Pixel>::type, inA + 1, inA + 1> outputArray =
                ndarray::allocate(outputShape);
        outputArray.deep() = 0;
        unflatten(outputArray, input, lsst::geom::Point2I(_bbox.getMinX(), _bbox.getMinY()));
        return outputArray;
    }

    template <typename PixelIn, typename PixelOut, int inA, int outC, int inC>
    void unflatten(ndarray::Array<PixelOut, inA + 1, outC> const &output,
                   ndarray::Array<PixelIn, inA, inC> const &input,
                   lsst::geom::Point2I const &xy0 = lsst::geom::Point2I()) const {
        // Populate 2D ndarray output with values from input, at locations defined by SpanSet, optionally
        // offset by xy0
        auto ndAssigner = [](lsst::geom::Point2I const &point,
                             typename details::ImageNdGetter<PixelOut, inA + 1, outC>::Reference out,
                             typename details::FlatNdGetter<PixelIn, inA, inC>::Reference in) { out = in; };
        applyFunctor(ndAssigner, ndarray::ndImage(output, xy0), ndarray::ndFlat(input));
    }

    template <typename ImageT>
    void copyImage(image::Image<ImageT> const &src, image::Image<ImageT> &dest) {
        auto copyFunc = [](lsst::geom::Point2I const &point, ImageT const &srcPix, ImageT &destPix) {
            destPix = srcPix;
        };
        applyFunctor(copyFunc, src, dest);
    }

    template <typename ImageT, typename MaskT, typename VarT>
    void copyMaskedImage(image::MaskedImage<ImageT, MaskT, VarT> const &src,
                         image::MaskedImage<ImageT, MaskT, VarT> &dest) {
        auto copyFunc = [](lsst::geom::Point2I const &point, ImageT const &srcPix, MaskT const &srcMask,
                           VarT const &srcVar, ImageT &destPix, MaskT &destMask, VarT &destVar) {
            destPix = srcPix;
            destMask = srcMask;
            destVar = srcVar;
        };
        applyFunctor(copyFunc, *(src.getImage()), *(src.getMask()), *(src.getVariance()), *(dest.getImage()),
                     *(dest.getMask()), *(dest.getVariance()));
    }

    template <typename ImageT>
    void setImage(image::Image<ImageT> &image, ImageT val,
                  lsst::geom::Box2I const &region = lsst::geom::Box2I(), bool doClip = false) const;

    template <typename Functor, typename... Args>
    // Normally std::forward would be used with a universal reference, however
    // this function does not use one because without std::forward the
    // compiler is forced to keep any r-value references alive for the
    // duration of the function call
    void applyFunctor(Functor &&func, Args &&... args) const {
        /* Use a variadic template to take a functor object, and an arbitrary number
           of parameters. For each of the arguments, construct a Getter class using
           a function (makeGetter) which is overloaded to all the types applyFunctorImpl
           supports: Images, MaskedImages, Exposures, ndarrays, numeric values, and
           iterators. The functor and the getters are then passed to the implementation of
           applyFunctor where the values of the input arguments are intelligently
           generated at each point in SpanSet, and passed to the functor object for evaluation.
         */
        applyFunctorImpl(func, details::makeGetter(args)...);
    }

    template <typename T>
    void setMask(lsst::afw::image::Mask<T> &target, T bitmask) const;

    template <typename T>
    void clearMask(lsst::afw::image::Mask<T> &target, T bitmask) const;

    // SpanSet functions
    std::shared_ptr<SpanSet> intersect(SpanSet const &other) const;

    template <typename T>
    std::shared_ptr<SpanSet> intersect(image::Mask<T> const &other, T bitmask) const;

    std::shared_ptr<SpanSet> intersectNot(SpanSet const &other) const;

    template <typename T>
    std::shared_ptr<SpanSet> intersectNot(image::Mask<T> const &other, T bitmask) const;

    std::shared_ptr<SpanSet> union_(SpanSet const &other) const;

    template <typename T>
    std::shared_ptr<SpanSet> union_(image::Mask<T> const &other, T bitmask) const;

    // Comparison Operators

    bool operator==(SpanSet const &other) const;

    /* Compute inequality between two SpanSets
     *
     * @param other The SpanSet for which inequality will be computed
     */
    bool operator!=(SpanSet const &other) const;

    static std::shared_ptr<geom::SpanSet> fromShape(int r, Stencil s = Stencil::CIRCLE,
                                                    lsst::geom::Point2I offset = lsst::geom::Point2I());

    static std::shared_ptr<geom::SpanSet> fromShape(geom::ellipses::Ellipse const &ellipse);

    template <typename T, typename UnaryPredicate = details::AnyBitSetFunctor<T>>
    static std::shared_ptr<geom::SpanSet> fromMask(
            image::Mask<T> const &mask, UnaryPredicate comparator = details::AnyBitSetFunctor<T>()) {
        // Create a vector which will hold all the spans created from the mask
        std::vector<Span> tempVec;
        // Grab some variables that will be used in the loop, so that they do not need to be fetched
        // every iteration.
        auto const maskArray = mask.getArray();
        auto const minPoint = mask.getBBox().getMin();
        auto const dimensions = maskArray.getShape();
        auto const minY = minPoint.getY();
        auto const minX = minPoint.getX();
        auto const dimMinusOne = dimensions[1] - 1;
        auto const yDim = dimensions[0];
        auto const xDim = dimensions[1];
        auto arrIter = maskArray.begin();
        for (size_t y = 0; y < yDim; ++y) {
            auto yWithOffset = y + minY;
            bool inSpan = false;  // are we currently in a span of interest?
            std::size_t start;    // starting x value of span we're currently in
            for (size_t x = 0; x < xDim; ++x) {
                bool compareValue = comparator((*arrIter)[x]);
                if (inSpan) {
                    // If we were in a span but now the condition isn't satisfied, it means the Span stops
                    // at the previous pixel.
                    if (!compareValue) {
                        tempVec.push_back(Span(yWithOffset, start + minX, x - 1 + minX));
                        inSpan = false;
                    }
                } else if (compareValue) {
                    // If we weren't in a span but now the condition is satisfied, we start a new span.
                    inSpan = true;
                    start = x;
                }
            }
            // Since the x loop is over, if we're still in a span, this means the Span was not
            // closed out and added to the vector. The last pixel should be included in the Span
            // and the Span should be closed and added to the vector of spans.
            if (inSpan) {
                tempVec.push_back(Span(yWithOffset, start + minX, dimMinusOne + minX));
            }
            ++arrIter;
        }

        // construct a SpanSet from the spans determined above
        return std::make_shared<SpanSet>(std::move(tempVec), false);
    }

    template <typename T>
    static std::shared_ptr<geom::SpanSet> fromMask(image::Mask<T> const &mask, T bitmask) {
        return fromMask(mask, [bitmask](T const &bitPattern) { return bitPattern & bitmask; });
    }

    std::vector<std::shared_ptr<geom::SpanSet>> split() const;

    bool isPersistable() const noexcept override { return true; }

    std::shared_ptr<geom::SpanSet> findEdgePixels() const;

private:
    /* Returns the name used by the persistence layer to identify the SpanSet class
     */
    std::string getPersistenceName() const override;

    /* Return a string corresponding to the python module that SpanSets lives in
     */
    inline std::string getPythonModule() const override { return "lsst.afw.geom"; }

    /* Writes the representation of the class out to an output archive
     */
    void write(OutputArchiveHandle &handle) const override;

    /* A class which is used by the persistence layer to restore SpanSets from an archive
     */
    friend class SpansSetFactory;

    /* A function to combine overlapping Spans in a SpanSet into a single Span
     */
    void _runNormalize();

    /* Initializes the SpanSet class. Contains code that is common to multiple constructors
     */
    void _initialize();

    /* Label Spans according to contiguous group. If the SpanSet is contiguous, all Spans will be labeled 1.
     * If there is more than one group each group will receive a label one higher than the previous.
     */
    void _label(geom::Span const &spn, std::vector<std::size_t> &labelVector, std::size_t currentLabel,
                std::unordered_map<int, std::vector<std::pair<std::size_t, const Span *>>> &sortVector) const;
    std::pair<std::vector<std::size_t>, std::size_t> _makeLabels() const;

    std::shared_ptr<SpanSet> makeShift(int x, int y) const;

    template <typename F, typename... T>
    void applyFunctorImpl(F &&f, T... args) const {
        /* Implementation for applying functors, loop over each of the spans, and then
         * each point. Use the get method in the getters to fetch the value and pass
         * the point, and the values to the functor
         */
        // make sure that the SpanSet is within the bounds of functor arguments
        details::variadicBoundChecker(_bbox, _area, args...);
        for (auto const &spn : _spanVector) {
            // Set the current span in the getter, useful for optimizing value lookups
            details::variadicSpanSetter(spn, args...);
            for (int x = spn.getX0(); x <= spn.getX1(); ++x) {
                lsst::geom::Point2I point(x, spn.getY());
                f(point, args.get()...);
                details::variadicIncrementPosition(args...);
            }
        }
    }

    // Vector to hold the Spans contained in the SpanSet
    std::vector<Span> _spanVector;

    // Box that is large enough to bound all pixels in the SpanSet
    lsst::geom::Box2I _bbox;

    // Number of pixels in the SpanSet
    std::size_t _area;
};
}  // namespace geom
}  // namespace afw
}  // namespace lsst

#endif  // LSST_AFW_GEOM_SPANSET_H