Mat3.h

Go to the documentation of this file.

// Copyright Contributors to the OpenVDB Project
// SPDX-License-Identifier: MPL-2.0

#ifndef OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED
#define OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED

#include <openvdb/Exceptions.h>
#include "Vec3.h"
#include "Mat.h"
#include <algorithm> // for std::copy()
#include <cassert>
#include <cmath>
#include <iomanip>


namespace openvdb {
OPENVDB_USE_VERSION_NAMESPACE
namespace OPENVDB_VERSION_NAME {
namespace math {

template<typename T> class Vec3;
template<typename T> class Mat4;
template<typename T> class Quat;

template<typename T>
class Mat3: public Mat<3, T>
{
public:
    using value_type = T;
    using ValueType = T;
    using MyBase = Mat<3, T>;
    Mat3() {}

    Mat3(const Quat<T> &q)
    { setToRotation(q); }



    template<typename Source>
    Mat3(Source a, Source b, Source c,
         Source d, Source e, Source f,
         Source g, Source h, Source i)
    {
        MyBase::mm[0] = static_cast<T>(a);
        MyBase::mm[1] = static_cast<T>(b);
        MyBase::mm[2] = static_cast<T>(c);
        MyBase::mm[3] = static_cast<T>(d);
        MyBase::mm[4] = static_cast<T>(e);
        MyBase::mm[5] = static_cast<T>(f);
        MyBase::mm[6] = static_cast<T>(g);
        MyBase::mm[7] = static_cast<T>(h);
        MyBase::mm[8] = static_cast<T>(i);
    } // constructor1Test

    template<typename Source>
    Mat3(const Vec3<Source> &v1, const Vec3<Source> &v2, const Vec3<Source> &v3, bool rows = true)
    {
        if (rows) {
            this->setRows(v1, v2, v3);
        } else {
            this->setColumns(v1, v2, v3);
        }
    }

    template<typename Source>
    Mat3(Source *a)
    {
        MyBase::mm[0] = static_cast<T>(a[0]);
        MyBase::mm[1] = static_cast<T>(a[1]);
        MyBase::mm[2] = static_cast<T>(a[2]);
        MyBase::mm[3] = static_cast<T>(a[3]);
        MyBase::mm[4] = static_cast<T>(a[4]);
        MyBase::mm[5] = static_cast<T>(a[5]);
        MyBase::mm[6] = static_cast<T>(a[6]);
        MyBase::mm[7] = static_cast<T>(a[7]);
        MyBase::mm[8] = static_cast<T>(a[8]);
    } // constructor1Test

    Mat3(const Mat<3, T> &m)
    {
        for (int i=0; i<3; ++i) {
            for (int j=0; j<3; ++j) {
                MyBase::mm[i*3 + j] = m[i][j];
            }
        }
    }

    template<typename Source>
    explicit Mat3(const Mat3<Source> &m)
    {
        for (int i=0; i<3; ++i) {
            for (int j=0; j<3; ++j) {
                MyBase::mm[i*3 + j] = static_cast<T>(m[i][j]);
            }
        }
    }

    explicit Mat3(const Mat4<T> &m)
    {
        for (int i=0; i<3; ++i) {
            for (int j=0; j<3; ++j) {
                MyBase::mm[i*3 + j] = m[i][j];
            }
        }
    }

    static const Mat3<T>& identity() {
        static const Mat3<T> sIdentity = Mat3<T>(
            1, 0, 0,
            0, 1, 0,
            0, 0, 1
        );
        return sIdentity;
    }

    static const Mat3<T>& zero() {
        static const Mat3<T> sZero = Mat3<T>(
            0, 0, 0,
            0, 0, 0,
            0, 0, 0
        );
        return sZero;
    }

    void setRow(int i, const Vec3<T> &v)
    {
        // assert(i>=0 && i<3);
        int i3 = i * 3;

        MyBase::mm[i3+0] = v[0];
        MyBase::mm[i3+1] = v[1];
        MyBase::mm[i3+2] = v[2];
    } // rowColumnTest

    Vec3<T> row(int i) const
    {
        // assert(i>=0 && i<3);
        return Vec3<T>((*this)(i,0), (*this)(i,1), (*this)(i,2));
    } // rowColumnTest

    void setCol(int j, const Vec3<T>& v)
    {
        // assert(j>=0 && j<3);
        MyBase::mm[0+j] = v[0];
        MyBase::mm[3+j] = v[1];
        MyBase::mm[6+j] = v[2];
    } // rowColumnTest

    Vec3<T> col(int j) const
    {
        // assert(j>=0 && j<3);
        return Vec3<T>((*this)(0,j), (*this)(1,j), (*this)(2,j));
    } // rowColumnTest

    T& operator()(int i, int j)
    {
        // assert(i>=0 && i<3);
        // assert(j>=0 && j<3);
        return MyBase::mm[3*i+j];
    } // trivial

    T operator()(int i, int j) const
    {
        // assert(i>=0 && i<3);
        // assert(j>=0 && j<3);
        return MyBase::mm[3*i+j];
    } // trivial

    void setRows(const Vec3<T> &v1, const Vec3<T> &v2, const Vec3<T> &v3)
    {
        MyBase::mm[0] = v1[0];
        MyBase::mm[1] = v1[1];
        MyBase::mm[2] = v1[2];
        MyBase::mm[3] = v2[0];
        MyBase::mm[4] = v2[1];
        MyBase::mm[5] = v2[2];
        MyBase::mm[6] = v3[0];
        MyBase::mm[7] = v3[1];
        MyBase::mm[8] = v3[2];
    } // setRows

    void setColumns(const Vec3<T> &v1, const Vec3<T> &v2, const Vec3<T> &v3)
    {
        MyBase::mm[0] = v1[0];
        MyBase::mm[1] = v2[0];
        MyBase::mm[2] = v3[0];
        MyBase::mm[3] = v1[1];
        MyBase::mm[4] = v2[1];
        MyBase::mm[5] = v3[1];
        MyBase::mm[6] = v1[2];
        MyBase::mm[7] = v2[2];
        MyBase::mm[8] = v3[2];
    } // setColumns

    void setSymmetric(const Vec3<T> &vdiag, const Vec3<T> &vtri)
    {
        MyBase::mm[0] = vdiag[0];
        MyBase::mm[1] = vtri[0];
        MyBase::mm[2] = vtri[1];
        MyBase::mm[3] = vtri[0];
        MyBase::mm[4] = vdiag[1];
        MyBase::mm[5] = vtri[2];
        MyBase::mm[6] = vtri[1];
        MyBase::mm[7] = vtri[2];
        MyBase::mm[8] = vdiag[2];
    } // setSymmetricTest

    static Mat3 symmetric(const Vec3<T> &vdiag, const Vec3<T> &vtri)
    {
        return Mat3(
                    vdiag[0], vtri[0], vtri[1],
                    vtri[0], vdiag[1], vtri[2],
                    vtri[1], vtri[2], vdiag[2]
                    );
    }

    void setSkew(const Vec3<T> &v)
    {*this = skew(v);}

    void setToRotation(const Quat<T> &q)
    {*this = rotation<Mat3<T> >(q);}

    void setToRotation(const Vec3<T> &axis, T angle)
    {*this = rotation<Mat3<T> >(axis, angle);}

    void setZero()
    {
        MyBase::mm[0] = 0;
        MyBase::mm[1] = 0;
        MyBase::mm[2] = 0;
        MyBase::mm[3] = 0;
        MyBase::mm[4] = 0;
        MyBase::mm[5] = 0;
        MyBase::mm[6] = 0;
        MyBase::mm[7] = 0;
        MyBase::mm[8] = 0;
    } // trivial

    void setIdentity()
    {
        MyBase::mm[0] = 1;
        MyBase::mm[1] = 0;
        MyBase::mm[2] = 0;
        MyBase::mm[3] = 0;
        MyBase::mm[4] = 1;
        MyBase::mm[5] = 0;
        MyBase::mm[6] = 0;
        MyBase::mm[7] = 0;
        MyBase::mm[8] = 1;
    } // trivial

    template<typename Source>
    const Mat3& operator=(const Mat3<Source> &m)
    {
        const Source *src = m.asPointer();

        // don't suppress type conversion warnings
        std::copy(src, (src + this->numElements()), MyBase::mm);
        return *this;
    } // opEqualToTest

    bool eq(const Mat3 &m, T eps=1.0e-8) const
    {
        return (isApproxEqual(MyBase::mm[0],m.mm[0],eps) &&
                isApproxEqual(MyBase::mm[1],m.mm[1],eps) &&
                isApproxEqual(MyBase::mm[2],m.mm[2],eps) &&
                isApproxEqual(MyBase::mm[3],m.mm[3],eps) &&
                isApproxEqual(MyBase::mm[4],m.mm[4],eps) &&
                isApproxEqual(MyBase::mm[5],m.mm[5],eps) &&
                isApproxEqual(MyBase::mm[6],m.mm[6],eps) &&
                isApproxEqual(MyBase::mm[7],m.mm[7],eps) &&
                isApproxEqual(MyBase::mm[8],m.mm[8],eps));
    } // trivial

    Mat3<T> operator-() const
    {
        return Mat3<T>(
                       -MyBase::mm[0], -MyBase::mm[1], -MyBase::mm[2],
                       -MyBase::mm[3], -MyBase::mm[4], -MyBase::mm[5],
                       -MyBase::mm[6], -MyBase::mm[7], -MyBase::mm[8]
                       );
    } // trivial

    // friend Mat3 operator*(T scalar, const Mat3& m) {
    //     return m*scalar;
    // }

    template <typename S>
    const Mat3<T>& operator*=(S scalar)
    {
        MyBase::mm[0] *= scalar;
        MyBase::mm[1] *= scalar;
        MyBase::mm[2] *= scalar;
        MyBase::mm[3] *= scalar;
        MyBase::mm[4] *= scalar;
        MyBase::mm[5] *= scalar;
        MyBase::mm[6] *= scalar;
        MyBase::mm[7] *= scalar;
        MyBase::mm[8] *= scalar;
        return *this;
    }

    template <typename S>
    const Mat3<T> &operator+=(const Mat3<S> &m1)
    {
        const S *s = m1.asPointer();

        MyBase::mm[0] += s[0];
        MyBase::mm[1] += s[1];
        MyBase::mm[2] += s[2];
        MyBase::mm[3] += s[3];
        MyBase::mm[4] += s[4];
        MyBase::mm[5] += s[5];
        MyBase::mm[6] += s[6];
        MyBase::mm[7] += s[7];
        MyBase::mm[8] += s[8];
        return *this;
    }

    template <typename S>
    const Mat3<T> &operator-=(const Mat3<S> &m1)
    {
        const S *s = m1.asPointer();

        MyBase::mm[0] -= s[0];
        MyBase::mm[1] -= s[1];
        MyBase::mm[2] -= s[2];
        MyBase::mm[3] -= s[3];
        MyBase::mm[4] -= s[4];
        MyBase::mm[5] -= s[5];
        MyBase::mm[6] -= s[6];
        MyBase::mm[7] -= s[7];
        MyBase::mm[8] -= s[8];
        return *this;
    }

    template <typename S>
    const Mat3<T> &operator*=(const Mat3<S> &m1)
    {
        Mat3<T> m0(*this);
        const T* s0 = m0.asPointer();
        const S* s1 = m1.asPointer();

        MyBase::mm[0] = static_cast<T>(s0[0] * s1[0] +
                                       s0[1] * s1[3] +
                                       s0[2] * s1[6]);
        MyBase::mm[1] = static_cast<T>(s0[0] * s1[1] +
                                       s0[1] * s1[4] +
                                       s0[2] * s1[7]);
        MyBase::mm[2] = static_cast<T>(s0[0] * s1[2] +
                                       s0[1] * s1[5] +
                                       s0[2] * s1[8]);

        MyBase::mm[3] = static_cast<T>(s0[3] * s1[0] +
                                       s0[4] * s1[3] +
                                       s0[5] * s1[6]);
        MyBase::mm[4] = static_cast<T>(s0[3] * s1[1] +
                                       s0[4] * s1[4] +
                                       s0[5] * s1[7]);
        MyBase::mm[5] = static_cast<T>(s0[3] * s1[2] +
                                       s0[4] * s1[5] +
                                       s0[5] * s1[8]);

        MyBase::mm[6] = static_cast<T>(s0[6] * s1[0] +
                                       s0[7] * s1[3] +
                                       s0[8] * s1[6]);
        MyBase::mm[7] = static_cast<T>(s0[6] * s1[1] +
                                       s0[7] * s1[4] +
                                       s0[8] * s1[7]);
        MyBase::mm[8] = static_cast<T>(s0[6] * s1[2] +
                                       s0[7] * s1[5] +
                                       s0[8] * s1[8]);

        return *this;
    }

    Mat3 cofactor() const
    {
        return Mat3<T>(
          MyBase::mm[4] * MyBase::mm[8] - MyBase::mm[5] * MyBase::mm[7],
          MyBase::mm[5] * MyBase::mm[6] - MyBase::mm[3] * MyBase::mm[8],
          MyBase::mm[3] * MyBase::mm[7] - MyBase::mm[4] * MyBase::mm[6],
          MyBase::mm[2] * MyBase::mm[7] - MyBase::mm[1] * MyBase::mm[8],
          MyBase::mm[0] * MyBase::mm[8] - MyBase::mm[2] * MyBase::mm[6],
          MyBase::mm[1] * MyBase::mm[6] - MyBase::mm[0] * MyBase::mm[7],
          MyBase::mm[1] * MyBase::mm[5] - MyBase::mm[2] * MyBase::mm[4],
          MyBase::mm[2] * MyBase::mm[3] - MyBase::mm[0] * MyBase::mm[5],
          MyBase::mm[0] * MyBase::mm[4] - MyBase::mm[1] * MyBase::mm[3]);
    }

    Mat3 adjoint() const
    {
        return Mat3<T>(
          MyBase::mm[4] * MyBase::mm[8] - MyBase::mm[5] * MyBase::mm[7],
          MyBase::mm[2] * MyBase::mm[7] - MyBase::mm[1] * MyBase::mm[8],
          MyBase::mm[1] * MyBase::mm[5] - MyBase::mm[2] * MyBase::mm[4],
          MyBase::mm[5] * MyBase::mm[6] - MyBase::mm[3] * MyBase::mm[8],
          MyBase::mm[0] * MyBase::mm[8] - MyBase::mm[2] * MyBase::mm[6],
          MyBase::mm[2] * MyBase::mm[3] - MyBase::mm[0] * MyBase::mm[5],
          MyBase::mm[3] * MyBase::mm[7] - MyBase::mm[4] * MyBase::mm[6],
          MyBase::mm[1] * MyBase::mm[6] - MyBase::mm[0] * MyBase::mm[7],
          MyBase::mm[0] * MyBase::mm[4] - MyBase::mm[1] * MyBase::mm[3]);

    } // adjointTest

    Mat3 transpose() const
    {
        return Mat3<T>(
          MyBase::mm[0], MyBase::mm[3], MyBase::mm[6],
          MyBase::mm[1], MyBase::mm[4], MyBase::mm[7],
          MyBase::mm[2], MyBase::mm[5], MyBase::mm[8]);

    } // transposeTest

    Mat3 inverse(T tolerance = 0) const
    {
        Mat3<T> inv(this->adjoint());

        const T det = inv.mm[0]*MyBase::mm[0] + inv.mm[1]*MyBase::mm[3] + inv.mm[2]*MyBase::mm[6];

        // If the determinant is 0, m was singular and the result will be invalid.
        if (isApproxEqual(det,T(0.0),tolerance)) {
            OPENVDB_THROW(ArithmeticError, "Inversion of singular 3x3 matrix");
        }
        return inv * (T(1)/det);
    } // invertTest

    T det() const
    {
        const T co00 = MyBase::mm[4]*MyBase::mm[8] - MyBase::mm[5]*MyBase::mm[7];
        const T co10 = MyBase::mm[5]*MyBase::mm[6] - MyBase::mm[3]*MyBase::mm[8];
        const T co20 = MyBase::mm[3]*MyBase::mm[7] - MyBase::mm[4]*MyBase::mm[6];
        return MyBase::mm[0]*co00  + MyBase::mm[1]*co10 + MyBase::mm[2]*co20;
    } // determinantTest

    T trace() const
    {
        return MyBase::mm[0]+MyBase::mm[4]+MyBase::mm[8];
    }

    Mat3 snapBasis(Axis axis, const Vec3<T> &direction)
    {
        return snapMatBasis(*this, axis, direction);
    }

    template<typename T0>
    Vec3<T0> transform(const Vec3<T0> &v) const
    {
        return static_cast< Vec3<T0> >(v * *this);
    } // xformVectorTest

    template<typename T0>
    Vec3<T0> pretransform(const Vec3<T0> &v) const
    {
        return static_cast< Vec3<T0> >(*this * v);
    } // xformTVectorTest


    Mat3 timesDiagonal(const Vec3<T>& diag) const
    {
        Mat3 ret(*this);

        ret.mm[0] *= diag(0);
        ret.mm[1] *= diag(1);
        ret.mm[2] *= diag(2);
        ret.mm[3] *= diag(0);
        ret.mm[4] *= diag(1);
        ret.mm[5] *= diag(2);
        ret.mm[6] *= diag(0);
        ret.mm[7] *= diag(1);
        ret.mm[8] *= diag(2);
        return ret;
    }
}; // class Mat3


template <typename T0, typename T1>
bool operator==(const Mat3<T0> &m0, const Mat3<T1> &m1)
{
    const T0 *t0 = m0.asPointer();
    const T1 *t1 = m1.asPointer();

    for (int i=0; i<9; ++i) {
        if (!isExactlyEqual(t0[i], t1[i])) return false;
    }
    return true;
}

template <typename T0, typename T1>
bool operator!=(const Mat3<T0> &m0, const Mat3<T1> &m1) { return !(m0 == m1); }

template <typename S, typename T>
Mat3<typename promote<S, T>::type> operator*(S scalar, const Mat3<T> &m)
{ return m*scalar; }

template <typename S, typename T>
Mat3<typename promote<S, T>::type> operator*(const Mat3<T> &m, S scalar)
{
    Mat3<typename promote<S, T>::type> result(m);
    result *= scalar;
    return result;
}

template <typename T0, typename T1>
Mat3<typename promote<T0, T1>::type> operator+(const Mat3<T0> &m0, const Mat3<T1> &m1)
{
    Mat3<typename promote<T0, T1>::type> result(m0);
    result += m1;
    return result;
}

template <typename T0, typename T1>
Mat3<typename promote<T0, T1>::type> operator-(const Mat3<T0> &m0, const Mat3<T1> &m1)
{
    Mat3<typename promote<T0, T1>::type> result(m0);
    result -= m1;
    return result;
}


template <typename T0, typename T1>
Mat3<typename promote<T0, T1>::type>operator*(const Mat3<T0> &m0, const Mat3<T1> &m1)
{
    Mat3<typename promote<T0, T1>::type> result(m0);
    result *= m1;
    return result;
}

template<typename T, typename MT>
inline Vec3<typename promote<T, MT>::type>
operator*(const Mat3<MT> &_m, const Vec3<T> &_v)
{
    MT const *m = _m.asPointer();
    return Vec3<typename promote<T, MT>::type>(
        _v[0]*m[0] + _v[1]*m[1] + _v[2]*m[2],
        _v[0]*m[3] + _v[1]*m[4] + _v[2]*m[5],
        _v[0]*m[6] + _v[1]*m[7] + _v[2]*m[8]);
}

template<typename T, typename MT>
inline Vec3<typename promote<T, MT>::type>
operator*(const Vec3<T> &_v, const Mat3<MT> &_m)
{
    MT const *m = _m.asPointer();
    return Vec3<typename promote<T, MT>::type>(
        _v[0]*m[0] + _v[1]*m[3] + _v[2]*m[6],
        _v[0]*m[1] + _v[1]*m[4] + _v[2]*m[7],
        _v[0]*m[2] + _v[1]*m[5] + _v[2]*m[8]);
}

template<typename T, typename MT>
inline Vec3<T> &operator *= (Vec3<T> &_v, const Mat3<MT> &_m)
{
    Vec3<T> mult = _v * _m;
    _v = mult;
    return _v;
}

template <typename T>
Mat3<T> outerProduct(const Vec3<T>& v1, const Vec3<T>& v2)
{
    return Mat3<T>(v1[0]*v2[0], v1[0]*v2[1], v1[0]*v2[2],
                   v1[1]*v2[0], v1[1]*v2[1], v1[1]*v2[2],
                   v1[2]*v2[0], v1[2]*v2[1], v1[2]*v2[2]);
}// outerProduct


template<typename T, typename T0>
Mat3<T> powLerp(const Mat3<T0> &m1, const Mat3<T0> &m2, T t)
{
    Mat3<T> x = m1.inverse() * m2;
    powSolve(x, x, t);
    Mat3<T> m = m1 * x;
    return m;
}


namespace mat3_internal {

template<typename T>
inline void
pivot(int i, int j, Mat3<T>& S, Vec3<T>& D, Mat3<T>& Q)
{
    const int& n = Mat3<T>::size;  // should be 3
    T temp;
    double cotan_of_2_theta;
    double tan_of_theta;
    double cosin_of_theta;
    double sin_of_theta;
    double z;

    double Sij = S(i,j);

    double Sjj_minus_Sii = D[j] - D[i];

    if (fabs(Sjj_minus_Sii) * (10*math::Tolerance<T>::value()) > fabs(Sij)) {
        tan_of_theta = Sij / Sjj_minus_Sii;
    } else {
        cotan_of_2_theta = 0.5*Sjj_minus_Sii / Sij ;

        if (cotan_of_2_theta < 0.) {
            tan_of_theta =
                -1./(sqrt(1. + cotan_of_2_theta*cotan_of_2_theta) - cotan_of_2_theta);
        } else {
            tan_of_theta =
                1./(sqrt(1. + cotan_of_2_theta*cotan_of_2_theta) + cotan_of_2_theta);
        }
    }

    cosin_of_theta = 1./sqrt( 1. + tan_of_theta * tan_of_theta);
    sin_of_theta = cosin_of_theta * tan_of_theta;
    z = tan_of_theta * Sij;
    S(i,j) = 0;
    D[i] -= z;
    D[j] += z;
    for (int k = 0; k < i; ++k) {
        temp = S(k,i);
        S(k,i) = cosin_of_theta * temp - sin_of_theta * S(k,j);
        S(k,j)= sin_of_theta * temp + cosin_of_theta * S(k,j);
    }
    for (int k = i+1; k < j; ++k) {
        temp = S(i,k);
        S(i,k) = cosin_of_theta * temp - sin_of_theta * S(k,j);
        S(k,j) = sin_of_theta * temp + cosin_of_theta * S(k,j);
    }
    for (int k = j+1; k < n; ++k) {
        temp = S(i,k);
        S(i,k) = cosin_of_theta * temp - sin_of_theta * S(j,k);
        S(j,k) = sin_of_theta * temp + cosin_of_theta * S(j,k);
    }
    for (int k = 0; k < n; ++k)
        {
            temp = Q(k,i);
            Q(k,i) = cosin_of_theta * temp - sin_of_theta*Q(k,j);
            Q(k,j) = sin_of_theta * temp + cosin_of_theta*Q(k,j);
        }
}

} // namespace mat3_internal


template<typename T>
inline bool
diagonalizeSymmetricMatrix(const Mat3<T>& input, Mat3<T>& Q, Vec3<T>& D,
    unsigned int MAX_ITERATIONS=250)
{
    Q  = Mat3<T>::identity();
    int n = Mat3<T>::size;  // should be 3

    Mat3<T> S(input);

    for (int i = 0; i < n; ++i) {
        D[i] = S(i,i);
    }

    unsigned int iterations(0);
    do {
        double er = 0;
        for (int i = 0; i < n; ++i) {
            for (int j = i+1; j < n; ++j) {
                er += fabs(S(i,j));
            }
        }
        if (std::abs(er) < math::Tolerance<T>::value()) {
            return true;
        }
        iterations++;

        T max_element = 0;
        int ip = 0;
        int jp = 0;
        for (int i = 0; i < n; ++i) {
            for (int j = i+1; j < n; ++j){

                if ( fabs(D[i]) * (10*math::Tolerance<T>::value()) > fabs(S(i,j))) {
                    S(i,j) = 0;
                }
                if (fabs(S(i,j)) > max_element) {
                    max_element = fabs(S(i,j));
                    ip = i;
                    jp = j;
                }
            }
        }
        mat3_internal::pivot(ip, jp, S, D, Q);
    } while (iterations < MAX_ITERATIONS);

    return false;
}

template<typename T>
inline Mat3<T>
Abs(const Mat3<T>& m)
{
    Mat3<T> out;
    const T* ip = m.asPointer();
    T* op = out.asPointer();
    for (unsigned i = 0; i < 9; ++i, ++op, ++ip) *op = math::Abs(*ip);
    return out;
}

template<typename Type1, typename Type2>
inline Mat3<Type1>
cwiseAdd(const Mat3<Type1>& m, const Type2 s)
{
    Mat3<Type1> out;
    const Type1* ip = m.asPointer();
    Type1* op = out.asPointer();
    for (unsigned i = 0; i < 9; ++i, ++op, ++ip) {
        OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
        *op = *ip + s;
        OPENVDB_NO_TYPE_CONVERSION_WARNING_END
    }
    return out;
}

template<typename T>
inline bool
cwiseLessThan(const Mat3<T>& m0, const Mat3<T>& m1)
{
    return cwiseLessThan<3, T>(m0, m1);
}

template<typename T>
inline bool
cwiseGreaterThan(const Mat3<T>& m0, const Mat3<T>& m1)
{
    return cwiseGreaterThan<3, T>(m0, m1);
}

using Mat3s = Mat3<float>;
using Mat3d = Mat3<double>;
using Mat3f = Mat3d;

} // namespace math


template<> inline math::Mat3s zeroVal<math::Mat3s>() { return math::Mat3s::zero(); }
template<> inline math::Mat3d zeroVal<math::Mat3d>() { return math::Mat3d::zero(); }

} // namespace OPENVDB_VERSION_NAME
} // namespace openvdb

#endif // OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED