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
/*
  Implements a compact representation of a collection of pixels
 */

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

namespace lsst { namespace afw { namespace geom { namespace details {
    /* Functor object to be used with maskToSpanSet function
     */
    template <typename T>
    class AnyBitSetFunctor {
     public:
        bool operator()(T const & pixelValue) {
            return pixelValue !=0;
        }
    };

}}}} // end lsst::afw::geom::details



namespace lsst { namespace afw { namespace geom {

enum class Stencil { CIRCLE, BOX, MANHATTAN };

// Forward declaration of the SpanSet class
class SpanSet;

template <typename T, typename UnaryPredicate = details::AnyBitSetFunctor<T>>
std::shared_ptr<geom::SpanSet> maskToSpanSet(image::Mask<T> const & mask,
                               UnaryPredicate p = details::AnyBitSetFunctor<T>() ) {
    std::vector<Span> tempVec;
    std::size_t startValue{0};
    bool started{false};
    auto const & maskArray = mask.getArray();
    auto dimensions = maskArray.getShape();
    for (size_t y = 0; y < dimensions[0]; ++y) {
        startValue = 0;
        started = false;
        for (size_t x = 0; x < dimensions[1]; ++x) {
            // If a new span has not been started, and a given x matches the functor condition
            // start a new span
            if (p(maskArray[y][x]) && !started) {
                started = true;
                startValue = x;
            }
            // If a span has been started, and the functor condition is false, that means the
            // Span being created should be stopped, and appended to the Span vector
            else if (started && !p(maskArray[y][x])) {
                tempVec.push_back(Span(y, startValue, x-1));
                started = false;
            }
            // If this is the last value in the Span's x range, and started is still true
            // that means the last value does not evaluate false in the functor and should be
            // included in the Span under construction. The Span should be completed and added
            // to the Span Vector before the next span is concidered.
            if (started && x == dimensions[1]) {
                tempVec.push_back(Span(y, startValue, x));
            }
        }
    }

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

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;
    const_iterator begin() const;
    const_iterator end() const;
    const_iterator cbegin() const;
    const_iterator cend() const;
    const_reference front() const;
    const_reference back() const;
    size_type size() const;
    bool empty() const;

    SpanSet();

    explicit SpanSet(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 = geom::Box2I();
            _area = 0;
        } else {
            if (normalize) {
                _runNormalize();
            }
            _initialize();
        }
    }

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


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

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

    // Define class methods
    size_type getArea() const;

    Box2I getBBox() const;

    bool isContiguous() const;

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

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

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

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

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

    bool overlaps(SpanSet const & other) const;

    bool contains(SpanSet const & other) const;

    bool contains(Point2I const & point) const;

    Point2D computeCentroid() const;

    ellipses::Quadrupole computeShape() const;

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

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

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

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

    template <typename Pixel, int inC>
    ndarray::Array<Pixel, 1, 1> flatten(ndarray::Array<Pixel, 2, inC>  & input,
                                        Point2I const & xy0 = Point2I()) const {
        // Populate a one dimensional array with the values from input taken at the points of SpanSet
        ndarray::Array<Pixel, 1, 1> outputArray = ndarray::allocate(ndarray::makeVector(getArea()));
        outputArray.deep() = 0;
        flatten(outputArray, input, xy0);
        return outputArray;
    }

    template <typename Pixel, int outC, int inC>
    void flatten(
            ndarray::Array<Pixel, 1, outC> const & output,
            ndarray::Array<Pixel, 2, inC> const & input,
            Point2I const & xy0 =  Point2I()) const {
        auto ndAssigner = []
                          (Point2I const & point,
                           typename details::FlatNdGetter<Pixel>::Reference out,
                           typename details::ImageNdGetter<Pixel, 2, 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 inC>
    ndarray::Array<Pixel, 2, 2> unflatten(ndarray::Array<Pixel, 1, inC> & input) const {
        // Create a two dimensional array the size of the bounding box. Populate values from input, placed at
        // locations corresponding to SpanSet, offset by the lower corner of the bounding box
        ndarray::Array<Pixel, 2, 2> outputArray = ndarray::allocate(_bBox.getHeight(), _bBox.getWidth());
        outputArray.deep() = 0;
        unflatten(outputArray, input, Point2I(_bBox.getMinX(), _bBox.getMinY()));
        return outputArray;
    }

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

    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 {
        // Use a lambda to set bits in a mask at the locations given by SpanSet
        auto targetArray = target.getArray();
        auto maskFunctor = []
                           (Point2I const & point,
                            typename details::ImageNdGetter<T, 2, 1>::Reference maskVal,
                            T bitmask)
                           {maskVal |= bitmask;};
        applyFunctor(maskFunctor, ndarray::ndImage(targetArray), bitmask);
    }

    template <typename T>
    void clearMask(lsst::afw::image::Mask<T> & target, T bitmask) const {
        // Use a lambda to clear bits in a mask at the locations given by SpanSet
        auto targetArray = target.getArray();
        auto clearMaskFunctor = []
                                (Point2I const & point,
                                 typename details::ImageNdGetter<T, 2, 1>::Reference maskVal,
                                 T bitmask)
                                {maskVal &= ~bitmask;};
        applyFunctor(clearMaskFunctor, ndarray::ndImage(targetArray), bitmask);
    }

    // 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 const & bitmask) const {
        auto comparator = [bitmask]
                          (T pixelValue)
                          {return (pixelValue & bitmask) == bitmask;};
        auto spanSetFromMask = geom::maskToSpanSet(other, comparator);
        return intersect(*spanSetFromMask);
    }

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

    template <typename T>
    std::shared_ptr<SpanSet> intersectNot(image::Mask<T> const & other, T const & bitmask) const {
        auto comparator = [bitmask]
                          (T pixelValue)
                          {return (pixelValue & bitmask) == bitmask;};
        auto spanSetFromMask = geom::maskToSpanSet(other, comparator);
        return intersectNot(*spanSetFromMask);
    }

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


    template <typename T>
    std::shared_ptr<SpanSet> union_(image::Mask<T> const & other, T const & bitmask) const {
        auto comparator = [bitmask]
                          (T pixelValue)
                          {return (pixelValue & bitmask) == bitmask;};
        auto spanSetFromMask = geom::maskToSpanSet(other, comparator);
        return union_(*spanSetFromMask);
    }

    // Comparison Operators
    bool operator==(SpanSet const & other) const;

    /* @brief 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> spanSetFromShape(int r, Stencil s = Stencil::CIRCLE);

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

    bool isPersistable() const { return true; }

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) const;
    std::pair<std::vector<std::size_t>, std::size_t> _makeLabels() const;

    /* Determine if two spans overlap
     *
     * a - First Span in comparison
     * b - Second Span in comparison
     * compareY - a boolean to control if the comparison takes into account the y position of the spans
     */
    inline bool spansOverlap(Span const & a, Span const & b, bool compareY = true) const {
        bool yTruth(true);
        if (compareY) {
            yTruth = a.getY() == b.getY();
        }
        return (yTruth && ((a.getMaxX() >= b.getMinX() && a.getMinX() <= b.getMinX()) ||
                (b.getMaxX() >= a.getMinX() && b.getMinX() <= a.getMinX())))
                ? true : false;
    }


    /* Determine if two spans are contiguous, that is they can overlap or the end of one span is
     * one pixel before the beginning of the next
     *
     * a - First Span in comparison
     * b - Second Span in comparison
     * compareY - a boolean to control if the comparison takes into account the y position of the spans
     */
    inline bool spansContiguous(Span const & a, Span const & b, bool compareY = true) const {
        bool yTruth(true);
        if (compareY) {
            yTruth = a.getY() == b.getY();
        }
        return (yTruth && ((a.getMaxX()+1 >= b.getMinX() && a.getMinX() <= b.getMinX()) ||
                (b.getMaxX()+1 >= a.getMinX() && b.getMinX() <= a. getMinX())))
                ? true: false;
    }
    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 (auto point : spn) {
                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
    Box2I _bBox;

    // Number of pixels in the SpanSet
    std::size_t _area;
};

 }}} // Close namespace lsst::afw::geom

#endif // LSST_AFW_GEOM_SPANSET_H