raster.hpp

Go to the documentation of this file.

#ifndef POPS_RASTER_HPP
#define POPS_RASTER_HPP
 
/*
 * PoPS model - native raster manipulation
 *
 * Copyright (C) 2015-2020 by the authors.
 *
 * Authors: Vaclav Petras <wenzeslaus gmail com>
 *          Completely rewritten by Vaclav Petras based on
 *          version by Zexi Chen <zchen22 ncsu edu>.
 *
 * The code contained herein is licensed under the GNU General Public
 * License. You may obtain a copy of the GNU General Public License
 * Version 2 or later at the following locations:
 *
 * http://www.opensource.org/licenses/gpl-license.html
 * http://www.gnu.org/copyleft/gpl.html
 */
 
#include <iostream>
#include <cmath>
#include <algorithm>
#include <stdexcept>
#include <initializer_list>
#include <type_traits>
 
namespace pops {
 
template<class InputIt1, class InputIt2, class BinaryOperation>
BinaryOperation
for_each_zip(InputIt1 first1, InputIt1 last1, InputIt2 first2, BinaryOperation f)
{
    for (; first1 != last1; ++first1, ++first2) {
        f(*first1, *first2);
    }
    return f;
}
 
template<typename Number, typename Index = int>
class Raster
{
protected:
    Index rows_;
    Index cols_;
    Number* data_;
    // owning is true for any state which is not using someone's data
    bool owns_;
 
public:
    typedef Number NumberType;
    typedef Index IndexType;
 
    Raster() : owns_(true)
    {
        cols_ = 0;
        rows_ = 0;
        data_ = NULL;
    }
 
    Raster(const Raster& other) : owns_(true)
    {
        cols_ = other.cols_;
        rows_ = other.rows_;
        data_ = new Number[cols_ * rows_];
        std::copy(other.data_, other.data_ + (cols_ * rows_), data_);
    }
 
    Raster(const Raster& other, Number value) : owns_(true)
    {
        cols_ = other.cols_;
        rows_ = other.rows_;
        data_ = new Number[cols_ * rows_]{value};
    }
 
    Raster(Raster&& other) : owns_(other.owns_)
    {
        cols_ = other.cols_;
        rows_ = other.rows_;
        data_ = other.data_;
        other.data_ = nullptr;
    }
 
    Raster(Index rows, Index cols) : owns_(true)
    {
        this->cols_ = cols;
        this->rows_ = rows;
        this->data_ = new Number[cols_ * rows_];
    }
 
    Raster(Index rows, Index cols, Number value) : owns_(true)
    {
        this->cols_ = cols;
        this->rows_ = rows;
        this->data_ = new Number[cols_ * rows_]{value};
    }
 
    Raster(Number* data, Index rows, Index cols)
        : rows_(rows), cols_(cols), data_(data), owns_(false)
    {}
 
    // maybe remove from the class, or make it optional together with
    // a reference
    Raster(std::initializer_list<std::initializer_list<Number>> l)
        : Raster(l.size(), l.begin()->size())
    {
        Index i = 0;
        Index j = 0;
        for (const auto& subl : l) {
            for (const auto& value : subl) {
                data_[cols_ * i + j] = value;
                ++j;
            }
            j = 0;
            ++i;
        }
    }
 
    ~Raster()
    {
        if (data_ && owns_) {
            delete[] data_;
        }
    }
 
    Index cols() const
    {
        return cols_;
    }
 
    Index rows() const
    {
        return rows_;
    }
 
    Number* data() noexcept
    {
        return data_;
    }
 
    const Number* data() const noexcept
    {
        return data_;
    }
 
    void fill(Number value)
    {
        std::fill(data_, data_ + (cols_ * rows_), value);
    }
 
    void zero()
    {
        std::fill(data_, data_ + (cols_ * rows_), 0);
    }
 
    template<class UnaryOperation>
    void for_each(UnaryOperation op) const
    {
        std::for_each(data_, data_ + (cols_ * rows_), op);
    }
 
    const Number& operator()(Index row, Index col) const
    {
        return data_[row * cols_ + col];
    }
 
    Number& operator()(Index row, Index col)
    {
        return data_[row * cols_ + col];
    }
 
    Raster& operator=(const Raster& other)
    {
        if (this != &other) {
            if (data_ && owns_)
                delete[] data_;
            cols_ = other.cols_;
            rows_ = other.rows_;
            data_ = new Number[cols_ * rows_];
            std::copy(other.data_, other.data_ + (cols_ * rows_), data_);
        }
        return *this;
    }
 
    Raster& operator=(Raster&& other)
    {
        if (this != &other) {
            if (data_ && owns_)
                delete[] data_;
            cols_ = other.cols_;
            rows_ = other.rows_;
            data_ = other.data_;
            owns_ = other.owns_;
            other.data_ = nullptr;
        }
        return *this;
    }
 
    template<typename OtherNumber>
    Raster& operator+=(OtherNumber value)
    {
        std::for_each(
            data_, data_ + (cols_ * rows_), [&value](Number& a) { a += value; });
        return *this;
    }
 
    template<typename OtherNumber>
    Raster& operator-=(OtherNumber value)
    {
        std::for_each(
            data_, data_ + (cols_ * rows_), [&value](Number& a) { a -= value; });
        return *this;
    }
 
    template<typename OtherNumber>
    Raster& operator*=(OtherNumber value)
    {
        std::for_each(
            data_, data_ + (cols_ * rows_), [&value](Number& a) { a *= value; });
        return *this;
    }
 
    template<typename OtherNumber>
    Raster& operator/=(OtherNumber value)
    {
        std::for_each(
            data_, data_ + (cols_ * rows_), [&value](Number& a) { a /= value; });
        return *this;
    }
 
    template<typename OtherNumber>
    typename std::enable_if<
        std::is_floating_point<Number>::value
            || std::is_same<Number, OtherNumber>::value,
        Raster&>::type
    operator+=(const Raster<OtherNumber>& image)
    {
        for_each_zip(
            data_,
            data_ + (cols_ * rows_),
            image.data(),
            [](Number& a, const OtherNumber& b) { a += b; });
        return *this;
    }
 
    template<typename OtherNumber>
    typename std::enable_if<
        std::is_floating_point<Number>::value
            || std::is_same<Number, OtherNumber>::value,
        Raster&>::type
    operator-=(const Raster<OtherNumber>& image)
    {
        for_each_zip(
            data_,
            data_ + (cols_ * rows_),
            image.data(),
            [](Number& a, const OtherNumber& b) { a -= b; });
        return *this;
    }
 
    template<typename OtherNumber>
    typename std::enable_if<
        std::is_floating_point<Number>::value
            || std::is_same<Number, OtherNumber>::value,
        Raster&>::type
    operator*=(const Raster<OtherNumber>& image)
    {
        for_each_zip(
            data_,
            data_ + (cols_ * rows_),
            image.data(),
            [](Number& a, const OtherNumber& b) { a *= b; });
        return *this;
    }
 
    template<typename OtherNumber>
    typename std::enable_if<
        std::is_floating_point<Number>::value
            || std::is_same<Number, OtherNumber>::value,
        Raster&>::type
    operator/=(const Raster<OtherNumber>& image)
    {
        for_each_zip(
            data_,
            data_ + (cols_ * rows_),
            image.data(),
            [](Number& a, const OtherNumber& b) { a /= b; });
        return *this;
    }
 
    bool operator==(const Raster& other) const
    {
        // TODO: assumes same sizes
        for (Index i = 0; i < cols_; i++) {
            for (Index j = 0; j < cols_; j++) {
                if (this->data_[i * cols_ + j] != other.data_[i * cols_ + j])
                    return false;
            }
        }
        return true;
    }
 
    bool operator!=(const Raster& other) const
    {
        // TODO: assumes same sizes
        for (Index i = 0; i < cols_; i++) {
            for (Index j = 0; j < cols_; j++) {
                if (this->data_[i * cols_ + j] != other.data_[i * cols_ + j])
                    return true;
            }
        }
        return false;
    }
 
    template<typename OtherNumber>
    friend inline
        typename std::enable_if<std::is_arithmetic<OtherNumber>::value, Raster>::type
        operator+(const Raster& raster, OtherNumber value)
    {
        auto out = Raster(raster.rows(), raster.cols());
 
        std::transform(
            raster.data(),
            raster.data() + (raster.cols() * raster.rows()),
            out.data(),
            [&value](const Number& a) { return a + value; });
        return out;
    }
 
    template<typename OtherNumber>
    friend inline
        typename std::enable_if<std::is_arithmetic<OtherNumber>::value, Raster>::type
        operator-(const Raster& raster, OtherNumber value)
    {
        auto out = Raster(raster.rows(), raster.cols());
 
        std::transform(
            raster.data(),
            raster.data() + (raster.cols() * raster.rows()),
            out.data(),
            [&value](const Number& a) { return a - value; });
        return out;
    }
 
    template<typename OtherNumber>
    friend inline
        typename std::enable_if<std::is_arithmetic<OtherNumber>::value, Raster>::type
        operator*(const Raster& raster, OtherNumber value)
    {
        auto out = Raster(raster.rows(), raster.cols());
 
        std::transform(
            raster.data(),
            raster.data() + (raster.cols() * raster.rows()),
            out.data(),
            [&value](const Number& a) { return a * value; });
        return out;
    }
 
    template<typename OtherNumber>
    friend inline
        typename std::enable_if<std::is_arithmetic<OtherNumber>::value, Raster>::type
        operator/(const Raster& raster, OtherNumber value)
    {
        auto out = Raster(raster.rows(), raster.cols());
 
        std::transform(
            raster.data(),
            raster.data() + (raster.cols() * raster.rows()),
            out.data(),
            [&value](const Number& a) { return a / value; });
        return out;
    }
 
    template<typename OtherNumber>
    friend inline
        typename std::enable_if<std::is_arithmetic<OtherNumber>::value, Raster>::type
        operator+(OtherNumber value, const Raster& raster)
    {
        return raster + value;
    }
 
    template<typename OtherNumber>
    friend inline
        typename std::enable_if<std::is_arithmetic<OtherNumber>::value, Raster>::type
        operator-(OtherNumber value, const Raster& raster)
    {
        auto out = Raster(raster.rows(), raster.cols());
 
        std::transform(
            raster.data(),
            raster.data() + (raster.cols() * raster.rows()),
            out.data(),
            [&value](const Number& a) { return value - a; });
        return out;
    }
 
    template<typename OtherNumber>
    friend inline
        typename std::enable_if<std::is_arithmetic<OtherNumber>::value, Raster>::type
        operator*(OtherNumber value, const Raster& raster)
    {
        return raster * value;
    }
 
    template<typename OtherNumber>
    friend inline
        typename std::enable_if<std::is_arithmetic<OtherNumber>::value, Raster>::type
        operator/(OtherNumber value, const Raster& raster)
    {
        auto out = Raster(raster.rows(), raster.cols());
 
        std::transform(
            raster.data(),
            raster.data() + (raster.cols() * raster.rows()),
            out.data(),
            [&value](const Number& a) { return value / a; });
        return out;
    }
 
    friend inline Raster pow(Raster image, double value)
    {
        image.for_each([value](Number& a) { a = std::pow(a, value); });
        return image;
    }
    friend inline Raster sqrt(Raster image)
    {
        image.for_each([](Number& a) { a = std::sqrt(a); });
        return image;
    }
 
    friend inline std::ostream& operator<<(std::ostream& stream, const Raster& image)
    {
        stream << "[[";
        for (Index i = 0; i < image.rows_; i++) {
            if (i != 0)
                stream << "],\n [";
            for (Index j = 0; j < image.cols_; j++) {
                if (j != 0)
                    stream << ", ";
                stream << image.data_[i * image.cols_ + j];
            }
        }
        stream << "]]\n";
        return stream;
    }
};
 
template<
    typename LeftNumber,
    typename RightNumber,
    typename ResultNumber = typename std::common_type<LeftNumber, RightNumber>::type>
Raster<ResultNumber>
operator+(const Raster<LeftNumber>& lhs, const Raster<RightNumber>& rhs)
{
    if (lhs.cols() != rhs.cols() || lhs.rows() != rhs.rows()) {
        throw std::invalid_argument(
            "Raster::operator+: The number of rows or columns does not match");
    }
    auto out = Raster<ResultNumber>(lhs.rows(), lhs.cols());
 
    std::transform(
        lhs.data(),
        lhs.data() + (lhs.cols() * lhs.rows()),
        rhs.data(),
        out.data(),
        [](const LeftNumber& a, const RightNumber& b) { return a + b; });
    return out;
}
 
template<
    typename LeftNumber,
    typename RightNumber,
    typename ResultNumber = typename std::common_type<LeftNumber, RightNumber>::type>
Raster<ResultNumber>
operator-(const Raster<LeftNumber>& lhs, const Raster<RightNumber>& rhs)
{
    if (lhs.cols() != rhs.cols() || lhs.rows() != rhs.rows()) {
        throw std::invalid_argument(
            "Raster::operator-: The number of rows or columns does not match");
    }
    auto out = Raster<ResultNumber>(lhs.rows(), lhs.cols());
 
    std::transform(
        lhs.data(),
        lhs.data() + (lhs.cols() * lhs.rows()),
        rhs.data(),
        out.data(),
        [](const LeftNumber& a, const RightNumber& b) { return a - b; });
    return out;
}
 
template<
    typename LeftNumber,
    typename RightNumber,
    typename ResultNumber = typename std::common_type<LeftNumber, RightNumber>::type>
Raster<ResultNumber>
operator*(const Raster<LeftNumber>& lhs, const Raster<RightNumber>& rhs)
{
    if (lhs.cols() != rhs.cols() || lhs.rows() != rhs.rows()) {
        throw std::invalid_argument(
            "Raster::operator*: The number of rows or columns does not match");
    }
    auto out = Raster<ResultNumber>(lhs.rows(), lhs.cols());
 
    std::transform(
        lhs.data(),
        lhs.data() + (lhs.cols() * lhs.rows()),
        rhs.data(),
        out.data(),
        [](const LeftNumber& a, const RightNumber& b) { return a * b; });
    return out;
}
 
template<
    typename LeftNumber,
    typename RightNumber,
    typename ResultNumber = typename std::common_type<LeftNumber, RightNumber>::type>
Raster<ResultNumber>
operator/(const Raster<LeftNumber>& lhs, const Raster<RightNumber>& rhs)
{
    if (lhs.cols() != rhs.cols() || lhs.rows() != rhs.rows()) {
        throw std::invalid_argument(
            "Raster::operator/: The number of rows or columns does not match");
    }
    auto out = Raster<ResultNumber>(lhs.rows(), lhs.cols());
 
    std::transform(
        lhs.data(),
        lhs.data() + (lhs.cols() * lhs.rows()),
        rhs.data(),
        out.data(),
        [](const LeftNumber& a, const RightNumber& b) { return a / b; });
    return out;
}
 
}  // namespace pops
 
#endif  // POPS_RASTER_HPP