vec.h

/********************************************************************
 * Copyright (C) 2015-2021 by Liangliang Nan <liangliang.nan@gmail.com>
 * Copyright (C) 2000-2005 INRIA - Project ALICE
 *
 * The code in this file is partly from OGF/Graphite (2.0 alpha-4) with
 * modifications and enhancement:
 *      https://gforge.inria.fr/forum/forum.php?forum_id=11459
 * The original code was distributed under the GNU GPL license.
 ********************************************************************/
 
#ifndef EASY3D_CORE_VEC_H
#define EASY3D_CORE_VEC_H
 
#include <cassert>
#include <iostream>
#include <cmath>
#include <cstring>
#include <limits>
 
 
namespace easy3d {

    template <size_t N, class T>
    class Vec {
    public:
        typedef Vec<N, T>   thisclass;

        Vec() { for (size_t i = 0; i < N; i++) { data_[i] = T(0); } }

        explicit Vec(const T& s) { for (size_t i = 0; i < N; i++) { data_[i] = s; } }

        // This one should never be called : a template constructor cannot be a copy constructor
        template<class T2> explicit Vec(const Vec<N, T2>& rhs) {
            for (size_t i = 0; i < N; i++) {
                data_[i] = T(rhs[i]);
            }
        }

        // to avoid compilation problems
        template<class T2, size_t M> explicit Vec(const Vec<M, T2>& rhs) {
            assert(M == N);
            for (size_t i = 0; i < N; i++) {
                data_[i] = T(rhs[i]);
            }
        }

        template<class T2> explicit Vec(const T2* rhs) {
            for (size_t i = 0; i < N; i++) {
                data_[i] = T(rhs[i]);
            }
        }

        thisclass& operator=(const thisclass& rhs) {
            // Check for self-assignment
            if (this == &rhs)
                return *this; // Return *this if rhs is the same as *this
            memcpy(data_, rhs.data(), N*sizeof(T));
            return *this;
        }

        static size_t dimension() { return (size_t)N; }
        static size_t size() { return (size_t)N; }

        T* data()             { return data_; }
        const T* data() const { return data_; }

        operator const T*() const { return data_; }
        operator T*() { return data_; }
 
        // Liangliang: The compiler can't decide whether to use your overloaded 
        //             operator[] or the built-in operator[] on the const T*.
        // See https://stackoverflow.com/questions/1726740/c-error-operator-2-overloads-have-similar-conversions 
        //T& operator[](size_t idx) {
        //    assert(idx < N);
        //    return data()[idx];
        //}
        //const T& operator[](size_t idx) const {
        //    assert(idx < N);
        //    return data()[idx];
        //}

        T length2() const {
            T result = T(0);
            for (size_t i = 0; i < N; i++) {
                result += data_[i] * data_[i];
            }
            return result;
        }

        T length() const {
            return std::sqrt(length2());
        }

        T norm() const {
            return length();
        }

        T distance2(const thisclass &rhs) const {
            T result = T(0);
            for (size_t i = 0; i < N; i++) {
                T val = rhs.data_[i] - data_[i];
                result += val*val;
            }
            return result;
        }

        thisclass& normalize() {
            T s = length();
            s = (s > std::numeric_limits<T>::min()) ? T(1.0) / s : T(0.0);
            *this *= s;
            return *this;
        }

        thisclass& operator+=(const thisclass& v) {
            for (size_t i = 0; i < N; i++) {
                data_[i] += v.data_[i];
            }
            return *this;
        }

        thisclass& operator-=(const thisclass& v) {
            for (size_t i = 0; i < N; i++) {
                data_[i] -= v.data_[i];
            }
            return *this;
        }

        thisclass& operator*=(const thisclass& v) {
            for (size_t i = 0; i < N; i++) {
                data_[i] *= v.data_[i];
            }
            return *this;
        }

        thisclass& operator/=(const thisclass& v) {
            for (size_t i = 0; i < N; i++) {
                data_[i] /= v.data_[i];
            }
            return *this;
        }

        template <class T2> thisclass& operator*=(T2 s) {
            for (size_t i = 0; i < N; i++) {
                data_[i] *= T(s);
            }
            return *this;
        }

        template <class T2> thisclass& operator/=(T2 s) {
            for (size_t i = 0; i < N; i++) {
                data_[i] /= T(s);
            }
            return *this;
        }

        thisclass operator+ (const thisclass& v) const {
            thisclass result(*this);
            for (size_t i = 0; i < N; i++) {
                result.data_[i] += v.data_[i];
            }
            return result;
        }

        thisclass operator- (const thisclass& v) const {
            thisclass result(*this);
            for (size_t i = 0; i < N; i++) {
                result.data_[i] -= v.data_[i];
            }
            return result;
        }

        template <class T2> thisclass operator* (T2 s) const {
            thisclass result(*this);
            for (size_t i = 0; i < N; i++) {
                result.data_[i] *= T(s);
            }
            return result;
        }

        template <class T2> thisclass operator/ (T2 s) const {
            thisclass result(*this);
            for (size_t i = 0; i < N; i++) {
                result.data_[i] /= T(s);
            }
            return result;
        }

        thisclass operator- () const {
            thisclass result;
            for (size_t i = 0; i < N; i++) {
                result.data_[i] = -data_[i];
            }
            return result;
        }
 
 
    private:
        T data_[N];
    };
 

    template <size_t N, class T> T dot(const Vec<N, T>& v1, const Vec<N, T>& v2) {
        T result = 0;
        for (size_t i = 0; i < N; i++) {
            result += v1[i] * v2[i];
        }
        return result;
    }

    template <size_t N, class T> Vec<N, T> operator-(const Vec<N, T>& v1) {
        Vec<N, T> result;
        for (size_t i = 0; i < N; i++) {
            result[i] = -v1[i];
        }
        return result;
    }

    template <class T2, size_t N, class T> Vec<N, T> operator*(T2 s, const Vec<N, T>& v) {
        Vec<N, T> result;
        for (size_t i = 0; i < N; i++) {
            result[i] = T(s)*v[i];
        }
        return result;
    }

    template <size_t N, class T> Vec<N, T> operator+(const Vec<N, T>& v1, const Vec<N, T>& v2) {
        Vec<N, T> result;
        for (size_t i = 0; i < N; i++) {
            result[i] = v1[i] + v2[i];
        }
        return result;
    }

    template <size_t N, class T> Vec<N, T> operator-(const Vec<N, T>& v1, const Vec<N, T>& v2) {
        Vec<N, T> result;
        for (size_t i = 0; i < N; i++) {
            result[i] = v1[i] - v2[i];
        }
        return result;
    }
 
    // Global functions for vectors (compatible with GLSL)

    template <size_t N, class T> T length(const Vec<N, T>& v) { return v.length(); }
    template <size_t N, class T> T norm(const Vec<N, T>& v) { return v.length(); }
    template <size_t N, class T> T length2(const Vec<N, T>& v) { return v.length2(); }
    template <size_t N, class T> T distance(const Vec<N, T>& v1, const Vec<N, T>& v2) { return length(v2 - v1); }
    template <size_t N, class T> T distance2(const Vec<N, T>& v1, const Vec<N, T>& v2)   { return v2.distance2(v1); }
    template <size_t N, class T> Vec<N, T> normalize(const Vec<N, T>& v) {
        T s = v.length();
        s = (s > std::numeric_limits<T>::min()) ? T(1.0) / s : T(0.0);
        return v * s;
    }

    template <size_t N, class T> Vec<N, T> mix(const Vec<N, T>& v1, const Vec<N, T>& v2, T w) {
        return (T(1) - w) * v1 + w * v2;
    }
 
    //-------------------- vec2 -------------------------------------------------------------------

    template <class T>
    class Vec<2, T> {
    public:
        typedef Vec<2, T>   thisclass;

        Vec() : x(0), y(0) { }

        Vec(T x_in, T y_in) : x(x_in), y(y_in) { }

        explicit Vec(const Vec<3, T>& v) : x(v.x), y(v.y) { }

        explicit Vec(const T& s) : x(s), y(s) { }

        template<class T2> explicit Vec(const Vec<2, T2> & v)
            : x(v.x), y(v.y) {}

        template<class T2> explicit Vec(const T2* v)
            : x(T(v[0])), y(T(v[1])) {}

        T length2() const { return x*x + y*y; }

        T length() const { return std::sqrt(x*x + y*y); }

        T norm() const {    return length(); }

        T distance2(const thisclass& rhs) const {
            T dx = rhs.x - x;
            T dy = rhs.y - y;
            return dx*dx + dy*dy;
        }

        thisclass& normalize() {
            T s = length();
            s = (s > std::numeric_limits<T>::min()) ? T(1.0) / s : T(0.0);
            *this *= s;
            return *this;
        }

        thisclass& operator+=(const thisclass& v) { x += v.x; y += v.y; return *this; }
        thisclass& operator-=(const thisclass& v) { x -= v.x; y -= v.y; return *this; }
        thisclass& operator*=(const thisclass& v) { x *= v.x; y *= v.y; return *this; }
        thisclass& operator/=(const thisclass& v) { x /= v.x; y /= v.y; return *this; }
        template <class T2> thisclass& operator*=(T2 s) { x *= T(s); y *= T(s); return *this; }
        template <class T2> thisclass& operator/=(T2 s) { x /= T(s); y /= T(s); return *this; }

        thisclass operator+ (const thisclass& v) const { return thisclass(x + v.x, y + v.y); }
        thisclass operator- (const thisclass& v) const { return thisclass(x - v.x, y - v.y); }

        template <class T2> thisclass operator* (T2 s) const { return thisclass(x*T(s), y*T(s)); }
        template <class T2> thisclass operator/ (T2 s) const { return thisclass(x / T(s), y / T(s)); }

        thisclass operator- () const { return thisclass(-x, -y); }

        static size_t dimension() { return (size_t)2; }
        static size_t size() { return (size_t)2; }

        T* data() { return _array; }
        const T* data() const { return _array; }

        operator const T*() const { return _array; }
        operator T*() { return _array; }
 
        // Liangliang: The compiler can't decide whether to use your overloaded 
        //             operator[] or the built-in operator[] on the const T*.
        // See https://stackoverflow.com/questions/1726740/c-error-operator-2-overloads-have-similar-conversions 
        //T& operator[](size_t idx) {
        //    assert(idx < 2);
        //    return _array[idx];
        //}
        //const T& operator[](size_t idx) const {
        //    assert(idx < 2);
        //    return _array[idx];
        //}

        union {
            T _array[2];        
            struct { T x, y; }; 
            struct { T u, v; }; 
        };
    };

    template <class T> T dot(const Vec<2, T>& v1, const Vec<2, T>& v2) { return v1.x*v2.x + v1.y*v2.y; }

    template <class T> T det(const Vec<2, T>& v1, const Vec<2, T>& v2) {
        return v1.x*v2.y - v1.y*v2.x;
    }

    template <class T> Vec<2, T> operator-(const Vec<2, T>& v1) {
        return Vec<2, T>(-v1.x, -v1.y);
    }

    template <class T2, class T> Vec<2, T> operator*(T2 s, const Vec<2, T>& v) {
        return Vec<2, T>(T(s)*v.x, T(s)*v.y);
    }

    template <class T> Vec<2, T> operator+(const Vec<2, T>& v1, const Vec<2, T>& v2) {
        return Vec<2, T>(v1.x + v2.x, v1.y + v2.y);
    }

    template <class T> Vec<2, T> operator-(const Vec<2, T>& v1, const Vec<2, T>& v2) {
        return Vec<2, T>(v1.x - v2.x, v1.y - v2.y);
    }
 
 
    //---------------- vec3 ------------------------------------------------------------------------

    template <class T>
    class Vec<3, T> {
    public:
        typedef Vec<3, T>   thisclass;

        Vec() : x(0), y(0), z(0) {}

        explicit Vec(const Vec<2, T>& v, const T& s = 1) : x(v.x), y(v.y), z(s) {}

        explicit Vec(const Vec<4, T>& v) : x(v.x), y(v.y), z(v.z) {}

        Vec(T x_in, T y_in, T z_in) : x(x_in), y(y_in), z(z_in) {}

        explicit Vec(const T& s) : x(s), y(s), z(s) {  }

        template<class T2> explicit Vec(const Vec<3, T2> & v) : x(v.x), y(v.y), z(v.z) {}

        template<class T2> explicit Vec(const T2* v)
            : x(T(v[0])), y(T(v[1])), z(T(v[2])) {}

        T length2() const { return x*x + y*y + z*z; }
        T length() const { return std::sqrt(x*x + y*y + z*z); }
        T norm() const { return length(); }
        T distance2(const thisclass& rhs) const {
            T dx = rhs.x - x;
            T dy = rhs.y - y;
            T dz = rhs.z - z;
            return dx*dx + dy*dy + dz*dz;
        }

        thisclass& normalize() {
            T s = length();
            s = (s > std::numeric_limits<T>::min()) ? T(1.0) / s : T(0.0);
            *this *= s;
            return *this;
        }

        thisclass& operator+=(const thisclass& v) { x += v.x; y += v.y; z += v.z; return *this; }
        thisclass& operator-=(const thisclass& v) { x -= v.x; y -= v.y; z -= v.z; return *this; }
        thisclass& operator*=(const thisclass& v) { x *= v.x; y *= v.y; z *= v.z; return *this; }
        thisclass& operator/=(const thisclass& v) { x /= v.x; y /= v.y; z /= v.z; return *this; }
        template <class T2> thisclass& operator*=(T2 s) { x *= T(s); y *= T(s); z *= T(s); return *this; }
        template <class T2> thisclass& operator/=(T2 s) { x /= T(s); y /= T(s); z /= T(s); return *this; }

        thisclass operator+ (const thisclass& v) const { return thisclass(x + v.x, y + v.y, z + v.z); }
        thisclass operator- (const thisclass& v) const { return thisclass(x - v.x, y - v.y, z - v.z); }

        template <class T2> thisclass operator* (T2 s) const { return thisclass(x*T(s), y*T(s), z*T(s)); }
        template <class T2> thisclass operator/ (T2 s) const { return thisclass(x / T(s), y / T(s), z / T(s)); }

        thisclass operator- () const { return thisclass(-x, -y, -z); }

        static size_t dimension() { return (size_t)3; }
        static size_t size() { return (size_t)3; }

        T* data() { return _array; }
        const T* data() const { return _array; }

        Vec<2, T> xy() const { return  Vec<2, T>(x, y); }

        operator const T*() const { return _array; }
        operator T*() { return _array; }
 
        // Liangliang: The compiler can't decide whether to use your overloaded 
        //             operator[] or the built-in operator[] on the const T*.
        // See https://stackoverflow.com/questions/1726740/c-error-operator-2-overloads-have-similar-conversions 
        //T& operator[](size_t idx) {
        //    assert(idx < 3);
        //    return _array[idx];
        //}
        //const T& operator[](size_t idx) const {
        //    assert(idx < 3);
        //    return _array[idx];
        //}

        union {
            T _array[3];            
            struct { T x, y, z; };  
            struct { T r, g, b; };  
        };
    };

    template <class T> T dot(const Vec<3, T>& v1, const Vec<3, T>& v2) {
        return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;
    }

    template <class T> Vec<3, T> cross(const Vec<3, T>& v1, const Vec<3, T>& v2) {
        return Vec<3, T>(
            v1.y*v2.z - v1.z*v2.y,
            v1.z*v2.x - v1.x*v2.z,
            v1.x*v2.y - v1.y*v2.x
            );
    }

    template <class T> Vec<3, T> operator-(const Vec<3, T>& v1) { return Vec<3, T>(-v1.x, -v1.y, -v1.z); }

    template <class T2, class T> Vec<3, T> operator*(T2 s, const Vec<3, T>& v) {
        return Vec<3, T>(T(s)*v.x, T(s)*v.y, T(s)*v.z);
    }

    template <class T> Vec<3, T> operator+(const Vec<3, T>& v1, const Vec<3, T>& v2) {
        return Vec<3, T>(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
    }

    template <class T> Vec<3, T> operator-(const Vec<3, T>& v1, const Vec<3, T>& v2) {
        return Vec<3, T>(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
    }

    template <class T> Vec<3, T> orthogonal(const Vec<3, T>& v) {
        T absx = std::fabs(v.x);
        T absy = std::fabs(v.y);
        T absz = std::fabs(v.z);
        // Find the smallest component. Keep equal case for null values.
        if ((absy >= absx) && (absz >= absx))
            return Vec<3, T>(0.0f, -v.z, v.y);
        else
            if ((absx >= absy) && (absz >= absy))
                return Vec<3, T>(-v.z, 0.0f, v.x);
            else
                return Vec<3, T>(-v.y, v.x, 0.0f);
    }
 
    // ----------------- vec4 ----------------------------------------------------------------------------------

    template <class T>
    class Vec<4, T> {
    public:
        typedef Vec<4, T>   thisclass;

        Vec() : x(0), y(0), z(0), w(0) {}

        explicit Vec(const Vec<3, T>& v, const T& s = 1) : x(v.x), y(v.y), z(v.z), w(s) {}

        Vec(T x_in, T y_in, T z_in, T w_in) : x(x_in), y(y_in), z(z_in), w(w_in) {}

        explicit Vec(const T& s) : x(s), y(s), z(s), w(s) {  }

        template<class T2> explicit Vec(const Vec<4, T2> & v)
            : x(v.x), y(v.y), z(v.z), w(v.w) {}

        template<class T2> explicit Vec(const T2* v)
            : x(v[0]), y(v[1]), z(v[2]), w(v[3]) {}

        T length2() const { return x*x + y*y + z*z + w*w; }
        T length() const { return std::sqrt(x*x + y*y + z*z + w*w); }
        T norm() const { return length(); }
        T distance2(const thisclass& rhs) const {
            T dx = rhs.x - x;
            T dy = rhs.y - y;
            T dz = rhs.z - z;
            T dw = rhs.w - w;
            return dx*dx + dy*dy + dz*dz + dw*dw;
        }

        thisclass& normalize() {
            T s = length();
            s = (s > std::numeric_limits<T>::min()) ? T(1.0) / s : T(0.0);
            *this *= s;
            return *this;
        }

        static size_t dimension() { return (size_t)4; }
        static size_t size() { return (size_t)4; }

        thisclass& operator+=(const thisclass& v) { x += v.x; y += v.y; z += v.z; w += v.w; return *this; }
        thisclass& operator-=(const thisclass& v) { x -= v.x; y -= v.y; z -= v.z; w -= v.w; return *this; }
        thisclass& operator*=(const thisclass& v) { x *= v.x; y *= v.y; z *= v.z; w *= v.w; return *this; }
        thisclass& operator/=(const thisclass& v) { x /= v.x; y /= v.y; z /= v.z; w /= v.w; return *this; }

        template <class T2> thisclass& operator*=(T2 s) {
            x *= T(s); y *= T(s); z *= T(s); w *= T(s); return *this;
        }
        template <class T2> thisclass& operator/=(T2 s) {
            x /= T(s); y /= T(s); z /= T(s); w /= T(s); return *this;
        }

        thisclass operator+ (const thisclass& v) const { return thisclass(x + v.x, y + v.y, z + v.z, w + v.w); }
        thisclass operator- (const thisclass& v) const { return thisclass(x - v.x, y - v.y, z - v.z, w - v.w); }

        template <class T2> thisclass operator* (T2 s) const { return thisclass(x*T(s), y*T(s), z*T(s), w*T(s)); }
        template <class T2> thisclass operator/ (T2 s) const { return thisclass(x / T(s), y / T(s), z / T(s), w / T(s)); }

        thisclass operator- () const { return thisclass(-x, -y, -z, -w); }

        T* data() { return _array; }
        const T* data() const { return _array; }

        Vec<3, T> xyz() const { return  Vec<3, T>(x, y, z); }

        operator const T*() const { return _array; }
        operator T*() { return _array; }
 
        // Liangliang: The compiler can't decide whether to use your overloaded 
        //             operator[] or the built-in operator[] on the const T*.
        // See https://stackoverflow.com/questions/1726740/c-error-operator-2-overloads-have-similar-conversions 
        //T& operator[](size_t idx) {
        //    assert(idx < 4);
        //    return _array[idx];
        //}
        //const T& operator[](size_t idx) const {
        //    assert(idx < 4);
        //    return _array[idx];
        //}

        union {
            T _array[4];                
            struct { T x, y, z, w; };   
            struct { T r, g, b, a; };   
        };
    };

    template <class T> T dot(const Vec<4, T>& v1, const Vec<4, T>& v2) {
        return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z + v1.w*v2.w;
    }

    template <class T> Vec<4, T> operator-(const Vec<4, T>& v1) { return Vec<4, T>(-v1.x, -v1.y, -v1.z, -v1.w); }

    template <class T2, class T> Vec<4, T> operator*(T2 s, const Vec<4, T>& v) {
        return Vec<4, T>(T(s)*v.x, T(s)*v.y, T(s)*v.z, T(s)*v.w);
    }

    template <class T> Vec<4, T> operator+(const Vec<4, T>& v1, const Vec<4, T>& v2) {
        return Vec<4, T>(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z, v1.w + v2.w);
    }

    template <class T> Vec<4, T> operator-(const Vec<4, T>& v1, const Vec<4, T>& v2) {
        return Vec<4, T>(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z, v1.w - v2.w);
    }
 
 
 
    //------------------------------- IO (input/output) ----------------------------

    template <size_t N, class T> std::ostream& operator<<(std::ostream& out, const Vec<N, T>& v) {
        for (size_t i = 0; i < N; i++) {
            out << v[i] << " ";
        }
        return out;
    }

    template <size_t N, class T> std::istream& operator>>(std::istream& in, Vec<N, T>& v) {
        for (size_t i = 0; i < N; i++) {
            in >> v[i];
        }
        return in;
    }

    template <class T> std::ostream& operator<<(std::ostream& out, const Vec<2, T>& v) {
        return out << v.x << " " << v.y;
    }

    template <class T> std::istream& operator>>(std::istream& in, Vec<2, T>& v) {
        return in >> v.x >> v.y;
    }

    template <class T> std::ostream& operator<<(std::ostream& out, const Vec<3, T>& v) {
        return out << v.x << " " << v.y << " " << v.z;
    }

    template <class T> std::istream& operator>>(std::istream& in, Vec<3, T>& v) {
        return in >> v.x >> v.y >> v.z;
    }

    template <class T> std::ostream& operator<<(std::ostream& out, const Vec<4, T>& v) {
        return out << v.x << " " << v.y << " " << v.z << " " << v.w;
    }

    template <class T> std::istream& operator>>(std::istream& in, Vec<4, T>& v) {
        return in >> v.x >> v.y >> v.z >> v.w;
    }
 
    //----------------------------------------------------------------------

    template <size_t N, class T>
    bool has_nan(const Vec<N, T> &v) {
        for (std::size_t i = 0; i < N; ++i) {
            if (std::isnan(v[i]) || std::isinf(v[i]))
                return true;
        }
        return false;
    }

    template <size_t N, class T>
    bool operator==(const Vec<N,T> &a, const Vec<N,T> &b) {
        for (size_t i = 0; i < N; i++) {
            if (a[i] != b[i]) return false;
        }
        return true;
    }

    template <size_t N, class T>
    bool operator!=(const Vec<N,T> &a, const Vec<N,T> &b) {
        for (size_t i = 0; i < N; i++) {
            if (a[i] != b[i]) return true;
        }
        return false;
    }

    template <size_t N, class T>
    bool operator<(const Vec<N,T> &a, const Vec<N,T> &b) {
        for(unsigned int i=0; i<N; ++i){
            if(a[i]<b[i]) return true;
            if(a[i]>b[i]) return false;
        }
        return false;
    }

    template <size_t N, class T>
    Vec<N, T> comp_product(const Vec<N, T> &v1, const Vec<N, T> &v2) {
        Vec<N, T> result;
        for (int i = 0; i < N; ++i)
            result[i] = v1[i] * v2[i];
        return result;
    }

    template <size_t N, class T>
    Vec<N, T> comp_min(const Vec<N, T> &v1, const Vec<N, T> &v2) {
        Vec<N, T> result;
        for (int i = 0; i < N; ++i)
            result[i] = std::min(v1[i], v2[i]);
        return result;
    }

    template <size_t N, class T>
    Vec<N, T> comp_max(const Vec<N, T> &v1, const Vec<N, T> &v2) {
        Vec<N, T> result;
        for (int i = 0; i < N; ++i)
            result[i] = std::max(v1[i], v2[i]);
        return result;
    }

    template<size_t N, class T>
    T min_coord(const Vec<N, T> &a) {
        T result = a[0];
        for (unsigned int i = 1; i < N; ++i) if (a[i] < result) result = a[i];
        return result;
    }

    template<size_t N, class T>
    T max_coord(const Vec<N, T> &a) {
        T result = a[0];
        for (unsigned int i = 1; i < N; ++i) if (a[i] > result) result = a[i];
        return result;
    }

    template<size_t N, class T>
    Vec<N, T> clamp(const Vec<N, T> &a, const Vec<N, T> &lower, const Vec<N, T> &upper) {
        Vec<N, T> result;
        for (unsigned int i = 0; i < N; ++i) result[i] = std::min(upper[i], std::max(a[i], lower[i]));
        return result;
    }
 
}
 
#endif  // EASY3D_CORE_VEC_H