$search

mat3x3.cc

Go to the documentation of this file.

/*
  James R. Diebel
  Stanford University

  Started: 3 December 2004

  mat3x3.cc - implementation for basic 3x3 matrix class

  Depends on vec3d.hh
*/

#include "bmtk/mat3x3.hh"

namespace bmtk {

// Constructors //

// identity matrix
Mat3x3::Mat3x3() {
  x[0] = x[4] = x[8] = 1.0f;
  x[1] = x[2] = x[3] = x[5] = x[6] = x[7] = 0;
  r[0].x = &x[0]; r[1].x = &x[3]; r[2].x = &x[6];
}

// infinite matrix in case of inverse of singular matrix
Mat3x3::Mat3x3(bool inf) {
  if (inf) x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = 1e16;
  else x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = -1e16;
  r[0].x = &x[0]; r[1].x = &x[3]; r[2].x = &x[6];
}

// each value passed
Mat3x3::Mat3x3(float x11, float x12, float x13,
         float x21, float x22, float x23,
         float x31, float x32, float x33) {
  x[0] = x11; x[1] = x12; x[2] = x13;
  x[3] = x21; x[4] = x22; x[5] = x23;
  x[6] = x31; x[7] = x32; x[8] = x33;
  r[0].x = &x[0]; r[1].x = &x[3]; r[2].x = &x[6];
}

// scalar multiply of identity matrix
Mat3x3::Mat3x3(const float a) {
  x[0] = x[4] = x[8] = a;
  x[1] = x[2] = x[3] = x[5] = x[6] = x[7] = 0;
  r[0].x = &x[0]; r[1].x = &x[3]; r[2].x = &x[6];
}

Mat3x3::Mat3x3(const float* x_) {
  x[0] = x_[0]; x[1] = x_[1]; x[2] = x_[2];
  x[3] = x_[3]; x[4] = x_[4]; x[5] = x_[5];
  x[6] = x_[6]; x[7] = x_[7]; x[8] = x_[8];
  r[0].x = &x[0]; r[1].x = &x[3]; r[2].x = &x[6];
}

// diagonal matrix defined by vector
Mat3x3::Mat3x3(const Vec3d& v) {
  x[0] = v.x[0]; x[4] = v.x[1]; x[8] = v.x[2];
  x[1] = x[2] = x[3] = x[5] = x[6] = x[7] = 0;
  r[0].x = &x[0]; r[1].x = &x[3]; r[2].x = &x[6];
}

// matrix formed as outer product of two vectors
Mat3x3::Mat3x3(const Vec3d& v0, const Vec3d& v1) {
  x[0] = v0.x[0]*v1.x[0]; x[1] = v0.x[0]*v1.x[1]; x[2] = v0.x[0]*v1.x[2];
  x[3] = v0.x[1]*v1.x[0]; x[4] = v0.x[1]*v1.x[1]; x[5] = v0.x[1]*v1.x[2];
  x[6] = v0.x[2]*v1.x[0]; x[7] = v0.x[2]*v1.x[1]; x[8] = v0.x[2]*v1.x[2];
  r[0].x = &x[0]; r[1].x = &x[3]; r[2].x = &x[6];
}

// columns of matrix passed as vectors
Mat3x3::Mat3x3(const Vec3d& v0, const Vec3d& v1, const Vec3d& v2) {
  x[0] = v0.x[0]; x[3] = v0.x[1]; x[6] = v0.x[2];
  x[1] = v1.x[0]; x[4] = v1.x[1]; x[7] = v1.x[2];
  x[2] = v2.x[0]; x[5] = v2.x[1]; x[8] = v2.x[2];
  r[0].x = &x[0]; r[1].x = &x[3]; r[2].x = &x[6];
}

Mat3x3::Mat3x3(const Vec3d& v0, const Vec3d& v1, const Vec3d& v2,
         const Vec3d& l) {
  x[0] = l.x[0]*v0.x[0]*v0.x[0]+l.x[1]*v1.x[0]*v1.x[0]+l.x[2]*v2.x[0]*v2.x[0];
  x[3] = l.x[0]*v0.x[1]*v0.x[0]+l.x[1]*v1.x[1]*v1.x[0]+l.x[2]*v2.x[1]*v2.x[0];
  x[4] = l.x[0]*v0.x[1]*v0.x[1]+l.x[1]*v1.x[1]*v1.x[1]+l.x[2]*v2.x[1]*v2.x[1];
  x[6] = l.x[0]*v0.x[2]*v0.x[0]+l.x[1]*v1.x[2]*v1.x[0]+l.x[2]*v2.x[2]*v2.x[0];
  x[7] = l.x[0]*v0.x[2]*v0.x[1]+l.x[1]*v1.x[2]*v1.x[1]+l.x[2]*v2.x[2]*v2.x[1];
  x[8] = l.x[0]*v0.x[2]*v0.x[2]+l.x[1]*v1.x[2]*v1.x[2]+l.x[2]*v2.x[2]*v2.x[2];
  x[1] = x[3]; x[2] = x[6]; x[5] = x[7];
  r[0].x = &x[0]; r[1].x = &x[3]; r[2].x = &x[6];
}
// Assignments //

// Assign an entire matrix
Mat3x3 Mat3x3::operator = (const Mat3x3& m) {
  x[0] = m.x[0]; x[1] = m.x[1]; x[2] = m.x[2];
  x[3] = m.x[3]; x[4] = m.x[4]; x[5] = m.x[5];
  x[6] = m.x[6]; x[7] = m.x[7]; x[8] = m.x[8];
  return *this;
}

// Reference operator: zero-based square bracket referencing
Row1x3 Mat3x3::operator [] (const int index) const {
  if (index < 0 || index > 2) cerr << "Index our of bounds" << endl << flush;
  return r[index];
}
Row1x3& Mat3x3::operator [] (const int index) {
  if (index < 0 || index > 2) cerr << "Index our of bounds" << endl << flush;
  return r[index];
}
float Row1x3::operator [] (const int index) const {
  if (index < 0 || index > 2) cerr << "Index our of bounds" << endl << flush;
  return x[index];
}
float& Row1x3::operator [] (const int index) {
  if (index < 0 || index > 2) cerr << "Index our of bounds" << endl << flush;
  return x[index];
}

// Addition, subtraction, multiplication, division //

// <matrix> +/- <matrix>
Mat3x3 Mat3x3::operator + (const Mat3x3& m) const {
  return Mat3x3(x[0]+m.x[0],x[1]+m.x[1],x[2]+m.x[2],
    x[3]+m.x[3],x[4]+m.x[4],x[5]+m.x[5],
    x[6]+m.x[6],x[7]+m.x[7],x[8]+m.x[8]);
}
void Mat3x3::operator += (const Mat3x3& m) {
  x[0] += m.x[0]; x[1] += m.x[1]; x[2] += m.x[2];
  x[3] += m.x[3]; x[4] += m.x[4]; x[5] += m.x[5];
  x[6] += m.x[6]; x[7] += m.x[7]; x[8] += m.x[8];
}
Mat3x3 Mat3x3::operator - (const Mat3x3& m) const {
  return Mat3x3(x[0]-m.x[0],x[1]-m.x[1],x[2]-m.x[2],
    x[3]-m.x[3],x[4]-m.x[4],x[5]-m.x[5],
    x[6]-m.x[6],x[7]-m.x[7],x[8]-m.x[8]);
}
void Mat3x3::operator -= (const Mat3x3& m) {
  x[0] -= m.x[0]; x[1] -= m.x[1]; x[2] -= m.x[2];
  x[3] -= m.x[3]; x[4] -= m.x[4]; x[5] -= m.x[5];
  x[6] -= m.x[6]; x[7] -= m.x[7]; x[8] -= m.x[8];
}

// <matrix> *// <scalar>
Mat3x3 Mat3x3::operator * (const float a) const {
  return Mat3x3(x[0]*a,x[1]*a,x[2]*a,
    x[3]*a,x[4]*a,x[5]*a,
    x[6]*a,x[7]*a,x[8]*a);
}
void Mat3x3::operator *= (const float a) {
  x[0] *= a; x[1] *= a; x[2] *= a;
  x[3] *= a; x[4] *= a; x[5] *= a;
  x[6] *= a; x[7] *= a; x[8] *= a;
}
Mat3x3 Mat3x3::operator / (const float a) const {
  return Mat3x3(x[0]/a,x[1]/a,x[2]/a,
    x[3]/a,x[4]/a,x[5]/a,
    x[6]/a,x[7]/a,x[8]/a);
}
void Mat3x3::operator /= (const float a) {
  x[0] /= a; x[1] /= a; x[2] /= a;
  x[3] /= a; x[4] /= a; x[5] /= a;
  x[6] /= a; x[7] /= a; x[8] /= a;
}
Mat3x3 operator * (const float a, const Mat3x3 &m) {
  return Mat3x3(a*m.x[0],a*m.x[1],a*m.x[2],
    a*m.x[3],a*m.x[4],a*m.x[5],
    a*m.x[6],a*m.x[7],a*m.x[8]);
}

// <matrix> * <vector>
Vec3d Mat3x3::operator * (const Vec3d& v) const {
  return Vec3d(x[0]*v.x[0]+x[1]*v.x[1]+x[2]*v.x[2],
         x[3]*v.x[0]+x[4]*v.x[1]+x[5]*v.x[2],
         x[6]*v.x[0]+x[7]*v.x[1]+x[8]*v.x[2]);
}
Vec3d operator * (const Vec3d& v, const Mat3x3& m) {
  return Vec3d(m.x[0]*v.x[0]+m.x[3]*v.x[1]+m.x[6]*v.x[2],
         m.x[1]*v.x[0]+m.x[4]*v.x[1]+m.x[7]*v.x[2],
         m.x[2]*v.x[0]+m.x[5]*v.x[1]+m.x[8]*v.x[2]);
}
void operator *= (Vec3d& v, const Mat3x3& m) {
  float v0 = m.x[0]*v.x[0]+m.x[3]*v.x[1]+m.x[6]*v.x[2];
  float v1 = m.x[1]*v.x[0]+m.x[4]*v.x[1]+m.x[7]*v.x[2];
  float v2 = m.x[2]*v.x[0]+m.x[5]*v.x[1]+m.x[8]*v.x[2];
  v.x[0] = v0; v.x[1] = v1; v.x[2] = v2;
}

// <matrix> * <matrix>
Mat3x3 Mat3x3::operator * (const Mat3x3& m) const {
  return Mat3x3(x[0]*m.x[0]+x[1]*m.x[3]+x[2]*m.x[6],
    x[0]*m.x[1]+x[1]*m.x[4]+x[2]*m.x[7],
    x[0]*m.x[2]+x[1]*m.x[5]+x[2]*m.x[8],
    x[3]*m.x[0]+x[4]*m.x[3]+x[5]*m.x[6],
    x[3]*m.x[1]+x[4]*m.x[4]+x[5]*m.x[7],
    x[3]*m.x[2]+x[4]*m.x[5]+x[5]*m.x[8],
    x[6]*m.x[0]+x[7]*m.x[3]+x[8]*m.x[6],
    x[6]*m.x[1]+x[7]*m.x[4]+x[8]*m.x[7],
    x[6]*m.x[2]+x[7]*m.x[5]+x[8]*m.x[8]);
}
void Mat3x3::operator *= (const Mat3x3& m) {
  float y[9];
  y[0] = x[0]; y[1] = x[1]; y[2] = x[2];
  y[3] = x[3]; y[4] = x[4]; y[5] = x[5];
  y[6] = x[6]; y[7] = x[7]; y[8] = x[8];
  x[0] = y[0]*m.x[0]+y[1]*m.x[3]+y[2]*m.x[6];
  x[1] = y[0]*m.x[1]+y[1]*m.x[4]+y[2]*m.x[7];
  x[2] = y[0]*m.x[2]+y[1]*m.x[5]+y[2]*m.x[8];
  x[3] = y[3]*m.x[0]+y[4]*m.x[3]+y[5]*m.x[6];
  x[4] = y[3]*m.x[1]+y[4]*m.x[4]+y[5]*m.x[7];
  x[5] = y[3]*m.x[2]+y[4]*m.x[5]+y[5]*m.x[8];
  x[6] = y[6]*m.x[0]+y[7]*m.x[3]+y[8]*m.x[6];
  x[7] = y[6]*m.x[1]+y[7]*m.x[4]+y[8]*m.x[7];
  x[8] = y[6]*m.x[2]+y[7]*m.x[5]+y[8]*m.x[8];
}

// Determinant, inverse, transpose, etc. //

// Determinant
float Mat3x3::det() const {
  return (x[0]*(x[4]*x[8] - x[7]*x[5]) -
    x[1]*(x[3]*x[8] - x[6]*x[5]) +
    x[2]*(x[3]*x[7] - x[6]*x[4]));
}

// Inverse
Mat3x3 Mat3x3::inv() const {
  float d = det();
  if (d != 0.0f) {
    Mat3x3 minv((x[4] * x[8] - x[5] * x[7]),
    (x[2] * x[7] - x[1] * x[8]),
    (x[1] * x[5] - x[2] * x[4]),
    (x[5] * x[6] - x[3] * x[8]),
    (x[0] * x[8] - x[2] * x[6]),
    (x[2] * x[3] - x[0] * x[5]),
    (x[3] * x[7] - x[4] * x[6]),
    (x[1] * x[6] - x[0] * x[7]),
    (x[0] * x[4] - x[1] * x[3]));
    minv /= d;
    return minv;
  }
  else {
    cout << "Warning: matrix singular, inverse set to 1e16" << endl << flush;
    return Mat3x3(true);
  }
}

// Transpose
void Mat3x3::transIt() {
  float temp;
  temp = x[1]; x[1] = x[3]; x[3] = temp;
  temp = x[2]; x[2] = x[6]; x[6] = temp;
  temp = x[5]; x[5] = x[7]; x[7] = temp;
}
Mat3x3 Mat3x3::trans() const {
  return Mat3x3(x[0],x[3],x[6],x[1],x[4],x[7],x[2],x[5],x[8]);
}

// Outer product
Mat3x3 Mat3x3::outer(const Vec3d& v0, const Vec3d& v1) {
  return Mat3x3(v0,v1);
}

// Eigenvector, eigenvalue reconstruction
Mat3x3 Mat3x3::eigenReconst(const Vec3d& v0, const Vec3d& v1,
           const Vec3d& v2, const Vec3d& l) {
  return Mat3x3(v0,v1,v2,l);
}

// Output functions //
void Mat3x3::print() {
  cout << "[" << x[0] << " " << x[1] << " " << x[2] << "]" << endl
       << "[" << x[3] << " " << x[4] << " " << x[5] << "]" << endl
       << "[" << x[6] << " " << x[7] << " " << x[8] << "]" << endl;
}

// Define << operator to handle output to screen
ostream& operator << (ostream& os, const Mat3x3& m) {
  return os << "[" << m.x[0] << " " << m.x[1] << " " << m.x[2] << "]" << endl
      << "[" << m.x[3] << " " << m.x[4] << " " << m.x[5] << "]" << endl
      << "[" << m.x[6] << " " << m.x[7] << " " << m.x[8] << "]" << endl;
}

} // namespace bmtk

---

bmtk
Author(s): Benjamin Pitzer
autogenerated on Mon Mar 4 11:05:51 2013