Mat22.h

Go to the documentation of this file.

/*
 * Mat22.h
 *
 *  Created on: Jul 25, 2019
 *      Author: mschefer
 */

#ifndef _SEUTILS_MAT22_H_
#define _SEUTILS_MAT22_H_

#include <tuple>
#include <assert.h>

namespace SExtractor {

class Mat22;
Mat22 operator*(const Mat22& a, const Mat22& b);

class Mat22 {
public:

  Mat22() : m_mat { 1.0, 0.0, 0.0, 1.0 } {
  }

  Mat22(double m1, double m2, double m3, double m4) : m_mat { m1, m2, m3, m4 } {
  }

  explicit Mat22(std::tuple<double, double, double, double> tuple) : m_mat {
    std::get<0>(tuple),
    std::get<1>(tuple),
    std::get<2>(tuple),
    std::get<3>(tuple)} {}

  Mat22(const Mat22&) = default;

  Mat22& operator=(const Mat22& other) {
    m_mat[0] = other[0];
    m_mat[1] = other[1];
    m_mat[2] = other[2];
    m_mat[3] = other[3];

    return *this;
  }

  double operator[](size_t i) const {
    return m_mat[i];
  }

  double& operator[](size_t i) {
    return m_mat[i];
  }

  Mat22& operator*=(const Mat22& other) {
    *this = *this * other;
    return *this;
  }

  Mat22 GetInverse() const {
    Mat22 out;

    double inv_det = 1.  / (m_mat[0] * m_mat[3] - m_mat[2] * m_mat[1]);
    out[0] = m_mat[3] * inv_det; out[1] = -m_mat[1] * inv_det;
    out[2] = -m_mat[2] * inv_det; out[3] = m_mat[0] * inv_det;

    assert(m_mat[0] * out[0] + m_mat[1] * out[2] >= 1 - 1e-8);
    assert(m_mat[0] * out[1] + m_mat[1] * out[3] <= 1e-8);
    assert(m_mat[2] * out[0] + m_mat[3] * out[2] <= 1e-8);
    assert(m_mat[2] * out[1] + m_mat[3] * out[3] >= 1 - 1e-8);

    return out;
  }

  Mat22 GetTranspose() const {
    Mat22 out { m_mat[0], m_mat[2],
                m_mat[1], m_mat[3] };

    return out;
  }

  double m_mat[4];
};

inline Mat22 operator*(const Mat22& a, const Mat22& b) {
  Mat22 out;

  out[0] = a[0] * b[0] + a[1] * b[2];
  out[1] = a[0] * b[1] + a[1] * b[3];
  out[2] = a[2] * b[0] + a[3] * b[2];
  out[3] = a[2] * b[1] + a[3] * b[3];

  return out;
}

inline Mat22 operator*(const Mat22& a, double b) {
  Mat22 out;

  out[0] = a[0] * b;
  out[1] = a[1] * b;
  out[2] = a[2] * b;
  out[3] = a[3] * b;

  return out;
}

}


#endif /* _SEUTILS_MAT22_H_ */