Vector11T.hh

/* ========================================================================= *
 *                                                                           *
 *                            OPENVOLUMEMESHMesh                                 *
 *           Copyright (c) 2001-2016, RWTH-Aachen University                 *
 *           Department of Computer Graphics and Multimedia                  *
 *                          All rights reserved.                             *
 *                         www.OPENVOLUMEMESHmesh.org                            *
 *                                                                           *
 *---------------------------------------------------------------------------*
 * This file is part of OPENVOLUMEMESHMesh.                                      *
 * This file was originally taken from OpenMesh                              *
 *---------------------------------------------------------------------------*
 *                                                                           *
 * Redistribution and use in source and binary forms, with or without        *
 * modification, are permitted provided that the following conditions        *
 * are met:                                                                  *
 *                                                                           *
 * 1. Redistributions of source code must retain the above copyright notice, *
 *    this list of conditions and the following disclaimer.                  *
 *                                                                           *
 * 2. Redistributions in binary form must reproduce the above copyright      *
 *    notice, this list of conditions and the following disclaimer in the    *
 *    documentation and/or other materials provided with the distribution.   *
 *                                                                           *
 * 3. Neither the name of the copyright holder nor the names of its          *
 *    contributors may be used to endorse or promote products derived from   *
 *    this software without specific prior written permission.               *
 *                                                                           *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS       *
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A           *
 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,  *
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,       *
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR        *
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF    *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING      *
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS        *
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.              *
 *                                                                           *
 * ========================================================================= */

#ifndef OPENVOLUMEMESH_SRC_OPENVOLUMEMESH_GEOMETRY_VECTOR11T_HH_
#define OPENVOLUMEMESH_SRC_OPENVOLUMEMESH_GEOMETRY_VECTOR11T_HH_

#include <array>
#include <utility>
#include <algorithm>
#include <numeric>
#include <type_traits>
#include <cmath>
#include <ostream>
#include <istream>
#include <cassert>
#include <cstdlib>

namespace OpenVolumeMesh {

namespace Geometry {

/*
 * Helpers for VectorT
 */
template<typename ... Ts>
struct are_convertible_to;

template<typename To, typename From, typename ... Froms>
struct are_convertible_to<To, From, Froms...> {
    static constexpr bool value = std::is_convertible<From, To>::value
            && are_convertible_to<To, Froms...>::value;
};

template<typename To, typename From>
struct are_convertible_to<To, From> : public std::is_convertible<From, To> {
};

template<typename Scalar, int DIM>
class VectorT {

        static_assert(DIM >= 1, "VectorT requires positive dimensionality.");

    private:
        using container = std::array<Scalar, DIM>;
        container values_;

    public:

        //---------------------------------------------------------------- class info

        typedef Scalar value_type;

        typedef VectorT<Scalar, DIM> vector_type;

        static constexpr int dim() {
            return DIM;
        }

        static constexpr size_t size() {
            return DIM;
        }

        static constexpr const size_t size_ = DIM;

        //-------------------------------------------------------------- constructors

        constexpr VectorT() {}

        explicit VectorT(const Scalar &v) {
            vectorize(v);
        }

        template<typename ... T,
            typename = typename std::enable_if<sizeof...(T) == DIM>::type,
            typename = typename std::enable_if<
                are_convertible_to<Scalar, T...>::value>::type>
        constexpr VectorT(T... vs) : values_ { {static_cast<Scalar>(vs)...} } {
            static_assert(sizeof...(T) == DIM,
                    "Invalid number of components specified in constructor.");
            static_assert(are_convertible_to<Scalar, T...>::value,
                    "Not all components are convertible to Scalar.");
        }

        VectorT(const VectorT &rhs) = default;
        VectorT(VectorT &&rhs) = default;
        VectorT &operator=(const VectorT &rhs) = default;
        VectorT &operator=(VectorT &&rhs) = default;

        template<typename S = Scalar, int D = DIM>
        auto homogenized() const ->
            typename std::enable_if<D == 4,
                VectorT<decltype(std::declval<S>()/std::declval<S>()), DIM>>::type {
            static_assert(D == DIM, "D and DIM need to be identical. (Never "
                    "override the default template arguments.)");
            static_assert(std::is_same<S, Scalar>::value, "S and Scalar need "
                    "to be the same type. (Never override the default template "
                    "arguments.)");
            return VectorT(
                    values_[0]/values_[3],
                    values_[1]/values_[3],
                    values_[2]/values_[3],
                    1);
        }

        template<typename Iterator,
            typename = decltype(
                    *std::declval<Iterator&>(), void(),
                    ++std::declval<Iterator&>(), void())>
        explicit VectorT(Iterator it) {
            std::copy_n(it, DIM, values_.begin());
        }

        template<typename otherScalarType,
            typename = typename std::enable_if<
                std::is_convertible<otherScalarType, Scalar>::value>>
        explicit VectorT(const VectorT<otherScalarType, DIM>& _rhs) {
            operator=(_rhs);
        }

        //--------------------------------------------------------------------- casts

        template<typename OtherScalar,
            typename = typename std::enable_if<
                std::is_convertible<OtherScalar, Scalar>::value>>
        vector_type& operator=(const VectorT<OtherScalar, DIM>& _rhs) {
            std::transform(_rhs.data(), _rhs.data() + DIM,
                    data(), [](OtherScalar rhs) {
                        return static_cast<Scalar>(std::move(rhs));
                    });
            return *this;
        }

        Scalar* data() { return values_.data(); }

        const Scalar *data() const { return values_.data(); }

        //----------------------------------------------------------- element access

        Scalar& operator[](size_t _i) {
            assert(_i < DIM);
            return values_[_i];
        }

        const Scalar& operator[](size_t _i) const {
            assert(_i < DIM);
            return values_[_i];
        }

        //---------------------------------------------------------------- comparsion

        bool operator==(const vector_type& _rhs) const {
            return std::equal(_rhs.values_.cbegin(), _rhs.values_.cend(), values_.cbegin());
        }

        bool operator!=(const vector_type& _rhs) const {
            return !std::equal(_rhs.values_.cbegin(), _rhs.values_.cend(), values_.cbegin());
        }

        //---------------------------------------------------------- scalar operators

        template<typename OtherScalar>
        auto operator*=(const OtherScalar& _s) ->
            typename std::enable_if<std::is_convertible<
                decltype(this->values_[0] * _s), Scalar>::value,
                VectorT<Scalar, DIM>&>::type {
            for (auto& e : *this) {
                e *= _s;
            }
            return *this;
        }

        template<typename OtherScalar>
        auto operator/=(const OtherScalar& _s) ->
            typename std::enable_if<std::is_convertible<
                decltype(this->values_[0] / _s), Scalar>::value,
                VectorT<Scalar, DIM>&>::type {
            for (auto& e : *this) {
                e /= _s;
            }
            return *this;
        }

        template<typename OtherScalar>
        typename std::enable_if<std::is_convertible<
                decltype(std::declval<Scalar>() * std::declval<OtherScalar>()),
                Scalar>::value,
            VectorT<Scalar, DIM>>::type
        operator*(const OtherScalar& _s) const {
            return vector_type(*this) *= _s;
        }

        template<typename OtherScalar>
        typename std::enable_if<std::is_convertible<
                decltype(std::declval<Scalar>() / std::declval<OtherScalar>()),
                Scalar>::value,
            VectorT<Scalar, DIM>>::type
        operator/(const OtherScalar& _s) const {
            return vector_type(*this) /= _s;
        }

        //---------------------------------------------------------- vector operators

        template<typename OtherScalar>
        auto operator*=(const VectorT<OtherScalar, DIM>& _rhs) ->
            typename std::enable_if<
                sizeof(decltype(this->values_[0] * *_rhs.data())) >= 0,
                vector_type&>::type {
            for (int i = 0; i < DIM; ++i) {
                data()[i] *= _rhs.data()[i];
            }
            return *this;
        }

        template<typename OtherScalar>
        auto operator/=(const VectorT<OtherScalar, DIM>& _rhs) ->
            typename std::enable_if<
                sizeof(decltype(this->values_[0] / *_rhs.data())) >= 0,
                vector_type&>::type {
            for (int i = 0; i < DIM; ++i) {
                data()[i] /= _rhs.data()[i];
            }
            return *this;
        }

        template<typename OtherScalar>
        auto operator-=(const VectorT<OtherScalar, DIM>& _rhs) ->
            typename std::enable_if<
                sizeof(decltype(this->values_[0] - *_rhs.data())) >= 0,
                vector_type&>::type {
            for (int i = 0; i < DIM; ++i) {
                data()[i] -= _rhs.data()[i];
            }
            return *this;
        }

        template<typename OtherScalar>
        auto operator+=(const VectorT<OtherScalar, DIM>& _rhs) ->
            typename std::enable_if<
                sizeof(decltype(this->values_[0] + *_rhs.data())) >= 0,
                vector_type&>::type {
            for (int i = 0; i < DIM; ++i) {
                data()[i] += _rhs.data()[i];
            }
            return *this;
        }

        template<typename OtherScalar>
        auto operator*(const VectorT<OtherScalar, DIM>& _rhs) const ->
            typename std::enable_if<
                sizeof(decltype(this->values_[0] * *_rhs.data())) >= 0,
                vector_type>::type {
            return vector_type(*this) *= _rhs;
        }

        template<typename OtherScalar>
        auto operator/(const VectorT<OtherScalar, DIM>& _rhs) const ->
            typename std::enable_if<
                sizeof(decltype(this->values_[0] / *_rhs.data())) >= 0,
                vector_type>::type {
            return vector_type(*this) /= _rhs;
        }

        template<typename OtherScalar>
        auto operator+(const VectorT<OtherScalar, DIM>& _rhs) const ->
            typename std::enable_if<
                sizeof(decltype(this->values_[0] + *_rhs.data())) >= 0,
                vector_type>::type {
            return vector_type(*this) += _rhs;
        }

        template<typename OtherScalar>
        auto operator-(const VectorT<OtherScalar, DIM>& _rhs) const ->
            typename std::enable_if<
                sizeof(decltype(this->values_[0] - *_rhs.data())) >= 0,
                vector_type>::type {
            return vector_type(*this) -= _rhs;
        }

        vector_type operator-(void) const {
            vector_type v;
            std::transform(values_.begin(), values_.end(), v.values_.begin(),
                    [](const Scalar &s) { return -s; });
            return v;
        }

        template<typename OtherScalar>
        auto operator% (const VectorT<OtherScalar, DIM> &_rhs) const ->
            typename std::enable_if<DIM == 3,
                VectorT<decltype(this->values_[0] * _rhs[0] -
                                 this->values_[0] * _rhs[0]),
                        DIM>>::type {
            return {
                values_[1] * _rhs[2] - values_[2] * _rhs[1],
                values_[2] * _rhs[0] - values_[0] * _rhs[2],
                values_[0] * _rhs[1] - values_[1] * _rhs[0]
            };
        }

        template<typename OtherScalar>
        auto operator|(const VectorT<OtherScalar, DIM>& _rhs) const ->
            decltype(*this->data() * *_rhs.data()) {

            return std::inner_product(data() + 1, data() + DIM, _rhs.data() + 1,
                    *data() * *_rhs.data());
        }

        //------------------------------------------------------------ euclidean norm



        template<typename S = Scalar>
        decltype(std::declval<S>() * std::declval<S>()) sqrnorm() const {
            static_assert(std::is_same<S, Scalar>::value, "S and Scalar need "
                    "to be the same type. (Never override the default template "
                    "arguments.)");
            typedef decltype(values_[0] * values_[0]) RESULT;
            return std::accumulate(values_.cbegin() + 1, values_.cend(),
                    values_[0] * values_[0],
                    [](const RESULT &l, const Scalar &r) { return l + r * r; });
        }

        template<typename S = Scalar>
        auto norm() const ->
            decltype(std::sqrt(std::declval<VectorT<S, DIM>>().sqrnorm())) {
            static_assert(std::is_same<S, Scalar>::value, "S and Scalar need "
                    "to be the same type. (Never override the default template "
                    "arguments.)");
            return std::sqrt(sqrnorm());
        }

        template<typename S = Scalar>
        auto length() const ->
            decltype(std::declval<VectorT<S, DIM>>().norm()) {
            static_assert(std::is_same<S, Scalar>::value, "S and Scalar need "
                    "to be the same type. (Never override the default template "
                    "arguments.)");
            return norm();
        }

        template<typename S = Scalar>
        auto normalize() ->
            decltype(*this /= std::declval<VectorT<S, DIM>>().norm()) {
            static_assert(std::is_same<S, Scalar>::value, "S and Scalar need "
                    "to be the same type. (Never override the default template "
                    "arguments.)");
            return *this /= norm();
        }

        template<typename S = Scalar>
        auto normalized() const ->
            decltype(*this / std::declval<VectorT<S, DIM>>().norm()) {
            static_assert(std::is_same<S, Scalar>::value, "S and Scalar need "
                    "to be the same type. (Never override the default template "
                    "arguments.)");
            return *this / norm();
        }

        template<typename S = Scalar>
        typename std::enable_if<
            sizeof(decltype(
                static_cast<S>(0),
                std::declval<VectorT<S, DIM>>().norm())) >= 0,
            vector_type&>::type
        normalize_cond() {
            static_assert(std::is_same<S, Scalar>::value, "S and Scalar need "
                    "to be the same type. (Never override the default template "
                    "arguments.)");
            auto n = norm();
            if (n != static_cast<decltype(norm())>(0)) {
                *this /= n;
            }
            return *this;
        }


        //------------------------------------------------------------ euclidean norm



        Scalar l1_norm() const {
            return std::accumulate(
                    values_.cbegin() + 1, values_.cend(), values_[0]);
        }

        Scalar l8_norm() const {
            return max_abs();
        }


        //------------------------------------------------------------ max, min, mean



        Scalar max() const {
            return *std::max_element(values_.cbegin(), values_.cend());
        }

        Scalar max_abs() const {
            return std::abs(
                *std::max_element(values_.cbegin(), values_.cend(),
                    [](const Scalar &a, const Scalar &b) {
                return std::abs(a) < std::abs(b);
            }));
        }

        Scalar min() const {
            return *std::min_element(values_.cbegin(), values_.cend());
        }

        Scalar min_abs() const {
            return std::abs(
                *std::min_element(values_.cbegin(), values_.cend(),
                    [](const Scalar &a, const Scalar &b) {
                return std::abs(a) < std::abs(b);
            }));
        }

        Scalar mean() const {
            return l1_norm()/DIM;
        }

        Scalar mean_abs() const {
            return std::accumulate(values_.cbegin() + 1, values_.cend(),
                    std::abs(values_[0]),
                    [](const Scalar &l, const Scalar &r) {
                        return l + std::abs(r);
                    }) / DIM;
        }

        vector_type& minimize(const vector_type& _rhs) {
            std::transform(values_.cbegin(), values_.cend(),
                    _rhs.values_.cbegin(),
                    values_.begin(),
                    [](const Scalar &l, const Scalar &r) {
                        return std::min(l, r);
                    });
            return *this;
        }

        bool minimized(const vector_type& _rhs) {
            bool result = false;
            std::transform(values_.cbegin(), values_.cend(),
                    _rhs.values_.cbegin(),
                    values_.begin(),
                    [&result](const Scalar &l, const Scalar &r) {
                        if (l < r) {
                            return l;
                        } else {
                            result = true;
                            return r;
                        }
                    });
            return result;
        }

        vector_type& maximize(const vector_type& _rhs) {
            std::transform(values_.cbegin(), values_.cend(),
                    _rhs.values_.cbegin(),
                    values_.begin(),
                    [](const Scalar &l, const Scalar &r) {
                        return std::max(l, r);
                    });
            return *this;
        }

        bool maximized(const vector_type& _rhs) {
            bool result = false;
            std::transform(values_.cbegin(), values_.cend(),
                    _rhs.values_.cbegin(),
                    values_.begin(),
                    [&result](const Scalar &l, const Scalar &r) {
                        if (l > r) {
                            return l;
                        } else {
                            result = true;
                            return r;
                        }
                    });
            return result;
        }

        inline vector_type min(const vector_type& _rhs) const {
            return vector_type(*this).minimize(_rhs);
        }

        inline vector_type max(const vector_type& _rhs) const {
            return vector_type(*this).maximize(_rhs);
        }


        //------------------------------------------------------------ misc functions

        template<typename Functor>
        inline vector_type apply(const Functor& _func) const {
            vector_type result;
            std::transform(result.values_.begin(), result.values_.end(),
                    result.values_.begin(), _func);
            return result;
        }

        vector_type& vectorize(const Scalar& _s) {
            std::fill(values_.begin(), values_.end(), _s);
            return *this;
        }

        static vector_type vectorized(const Scalar& _s) {
            return vector_type().vectorize(_s);
        }

        bool operator<(const vector_type& _rhs) const {
            return std::lexicographical_compare(
                    values_.begin(), values_.end(),
                    _rhs.values_.begin(), _rhs.values_.end());
        }

        void swap(VectorT& _other)
        noexcept(noexcept(std::swap(values_, _other.values_))) {
            std::swap(values_, _other.values_);
        }

        //------------------------------------------------------------ component iterators



        using iterator               = typename container::iterator;
        using const_iterator         = typename container::const_iterator;
        using reverse_iterator       = typename container::reverse_iterator;
        using const_reverse_iterator = typename container::const_reverse_iterator;

        iterator               begin()         noexcept { return values_.begin();   }
        const_iterator         begin()   const noexcept { return values_.cbegin();  }
        const_iterator         cbegin()  const noexcept { return values_.cbegin();  }

        iterator               end()           noexcept { return values_.end();     }
        const_iterator         end()     const noexcept { return values_.cend();    }
        const_iterator         cend()    const noexcept { return values_.cend();    }

        reverse_iterator       rbegin()        noexcept { return values_.rbegin();  }
        const_reverse_iterator rbegin()  const noexcept { return values_.crbegin(); }
        const_reverse_iterator crbegin() const noexcept { return values_.crbegin(); }

        reverse_iterator       rend()          noexcept { return values_.rend();    }
        const_reverse_iterator rend()    const noexcept { return values_.crend();   }
        const_reverse_iterator crend()   const noexcept { return values_.crend();   }

};

template<typename Scalar, int DIM, typename OtherScalar>
auto operator*(const OtherScalar& _s, const VectorT<Scalar, DIM> &rhs) ->
    decltype(rhs.operator*(_s)) {

    return rhs * _s;
}

template<typename Scalar, int DIM>
auto operator<<(std::ostream& os, const VectorT<Scalar, DIM> &_vec) ->
    typename std::enable_if<
        sizeof(decltype(os << _vec[0])) >= 0, std::ostream&>::type {

    os << _vec[0];
    for (int i = 1; i < DIM; ++i) {
        os << " " << _vec[i];
    }
    return os;
}

template<typename Scalar, int DIM>
auto operator>> (std::istream& is, VectorT<Scalar, DIM> &_vec) ->
    typename std::enable_if<
        sizeof(decltype(is >> _vec[0])) >= 0, std::istream &>::type {
    for (int i = 0; i < DIM; ++i)
        is >> _vec[i];
    return is;
}

template<typename Scalar, int DIM>
Scalar dot(const VectorT<Scalar, DIM>& _v1, const VectorT<Scalar, DIM>& _v2) {
  return (_v1 | _v2);
}

template<typename LScalar, typename RScalar, int DIM>
auto
cross(const VectorT<LScalar, DIM>& _v1, const VectorT<RScalar, DIM>& _v2) ->
    decltype(_v1 % _v2) {
  return (_v1 % _v2);
}

template<typename Scalar, int DIM>
void swap(VectorT<Scalar, DIM>& _v1, VectorT<Scalar, DIM>& _v2)
noexcept(noexcept(_v1.swap(_v2))) {
    _v1.swap(_v2);
}

//== TYPEDEFS =================================================================

typedef VectorT<signed char,1> Vec1c;
typedef VectorT<unsigned char,1> Vec1uc;
typedef VectorT<signed short int,1> Vec1s;
typedef VectorT<unsigned short int,1> Vec1us;
typedef VectorT<signed int,1> Vec1i;
typedef VectorT<unsigned int,1> Vec1ui;
typedef VectorT<float,1> Vec1f;
typedef VectorT<double,1> Vec1d;

typedef VectorT<signed char,2> Vec2c;
typedef VectorT<unsigned char,2> Vec2uc;
typedef VectorT<signed short int,2> Vec2s;
typedef VectorT<unsigned short int,2> Vec2us;
typedef VectorT<signed int,2> Vec2i;
typedef VectorT<unsigned int,2> Vec2ui;
typedef VectorT<float,2> Vec2f;
typedef VectorT<double,2> Vec2d;

typedef VectorT<signed char,3> Vec3c;
typedef VectorT<unsigned char,3> Vec3uc;
typedef VectorT<signed short int,3> Vec3s;
typedef VectorT<unsigned short int,3> Vec3us;
typedef VectorT<signed int,3> Vec3i;
typedef VectorT<unsigned int,3> Vec3ui;
typedef VectorT<float,3> Vec3f;
typedef VectorT<double,3> Vec3d;
typedef VectorT<bool,3> Vec3b;

typedef VectorT<signed char,4> Vec4c;
typedef VectorT<unsigned char,4> Vec4uc;
typedef VectorT<signed short int,4> Vec4s;
typedef VectorT<unsigned short int,4> Vec4us;
typedef VectorT<signed int,4> Vec4i;
typedef VectorT<unsigned int,4> Vec4ui;
typedef VectorT<float,4> Vec4f;
typedef VectorT<double,4> Vec4d;

typedef VectorT<signed char, 5> Vec5c;
typedef VectorT<unsigned char, 5> Vec5uc;
typedef VectorT<signed short int, 5> Vec5s;
typedef VectorT<unsigned short int, 5> Vec5us;
typedef VectorT<signed int, 5> Vec5i;
typedef VectorT<unsigned int, 5> Vec5ui;
typedef VectorT<float, 5> Vec5f;
typedef VectorT<double, 5> Vec5d;

typedef VectorT<signed char,6> Vec6c;
typedef VectorT<unsigned char,6> Vec6uc;
typedef VectorT<signed short int,6> Vec6s;
typedef VectorT<unsigned short int,6> Vec6us;
typedef VectorT<signed int,6> Vec6i;
typedef VectorT<unsigned int,6> Vec6ui;
typedef VectorT<float,6> Vec6f;
typedef VectorT<double,6> Vec6d;


} // namespace Geometry

using namespace Geometry;

template <class T>
const std::string typeName();

template <> const std::string typeName<Vec2f>();
template <> const std::string typeName<Vec2d>();
template <> const std::string typeName<Vec2i>();
template <> const std::string typeName<Vec2ui>();

template <> const std::string typeName<Vec3f>();
template <> const std::string typeName<Vec3d>();
template <> const std::string typeName<Vec3i>();
template <> const std::string typeName<Vec3ui>();

template <> const std::string typeName<Vec4f>();
template <> const std::string typeName<Vec4d>();
template <> const std::string typeName<Vec4i>();
template <> const std::string typeName<Vec4ui>();

} // namespace OpenVolumeMesh


#endif /* OPENVOLUMEMESH_SRC_OPENVOLUMEMESH_GEOMETRY_VECTOR11T_HH_ */