matrix3.h

Go to the documentation of this file.

/*
 * Copyright (c) 2011, Mårten Björkman (celle@csc.kth.se) 
 * All rights reserved.
 *
 * 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.The name of Mårten Björkman may not 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 MATRIX3_H
#define MATRIX3_H

#include <iostream>

class Vector2 {
public:
  double x[2];

  Vector2() {
    x[0] = x[1] = 0.0;
  }


  inline double &operator[](int idx) {
    return x[idx];
  } 


  inline double operator[](int idx) const {
    return x[idx];
  } 


  inline double &operator()(int idx) {
    return x[idx];
  } 


  Vector2 operator*(double v) {
    Vector2 res;
    for (int i=0;i<2;i++)
      res.x[i] = x[i]*v;
    return res;
  }


  Vector2 &operator=(double v) {
    for (int i=0;i<2;i++)
      x[i] = v;
    return *this;
  }
  

  Vector2 &operator*=(double v) {
    for (int i=0;i<2;i++)
      x[i] *= v;
    return *this;
  }
  

  Vector2 &operator+=(double v) {
    for (int i=0;i<2;i++)
      x[i] += v;
    return *this;
  }
  

  Vector2 &operator+=(const Vector2 &v) {
    for (int i=0;i<2;i++)
      x[i] += v[i];
    return *this;
  }
  

  friend std::ostream &operator<<(std::ostream &os, const Vector2 &v) {
    os << "[ " << v[0] << " " << v[1] << " ]";
    return os;
  }
};

class Vector3 {
public:
  double x[3];

  Vector3() {
    x[0] = x[1] = x[2] = 0.0;
  }


  inline double &operator[](int idx) {
    return x[idx];
  } 
  

  inline double operator[](int idx) const {
    return x[idx];
  } 


  inline double &operator()(int idx) {
    return x[idx];
  } 
  

  Vector3 operator*(double v) {
    Vector3 res;
    for (int i=0;i<3;i++)
      res.x[i] = x[i]*v;
    return res;
  }


  Vector3 &operator=(double v) {
    for (int i=0;i<3;i++)
      x[i] = v;
    return *this;
  }


  Vector3 &operator*=(double v) {
    for (int i=0;i<3;i++)
      x[i] *= v;
    return *this;
  }


  Vector3 &operator+=(double v) {
    for (int i=0;i<3;i++)
      x[i] += v;
    return *this;
  }


  Vector3 &operator+=(const Vector3 &v) {
    for (int i=0;i<3;i++)
      x[i] += v[i];
    return *this;
  }


  friend std::ostream &operator<<(std::ostream &os, const Vector3 &v) {
    os << "[ " << v[0] << " " << v[1] << " " << v[2] << " ]";
    return os;
  }
};

class Matrix2 {
public:
  Vector2 x[2];

  Matrix2() {  
    identity();
  }


  Matrix2(const Matrix2 &m) {
    for (int i=0;i<2;i++) 
      x[i] = m.x[i];
  }


  Matrix2 &operator=(const Matrix2 &m) {
    for (int i=0;i<2;i++) 
      x[i] = m.x[i];
    return *this;
  }

  void identity() {
    x[0][0] = x[1][1] = 1.0;
    x[0][1] = x[1][0] = 0.0;
  }


  inline Vector2 &operator[](int idx) {
    return x[idx];
  } 


  inline Vector2 operator[](int idx) const {
    return x[idx];
  } 


  inline double &operator()(int r, int c) {
    return x[r][c];
  }


  double determinant() {
    return x[0][0]*x[1][1] - x[0][1]*x[1][0];
  }


  Matrix2 invert() {
    double det = determinant();
    Matrix2 res;
    if (det!=0.0) {
      double idet = 1.0/det;
      res.x[0][0] =  idet*x[1][1];
      res.x[0][1] = -idet*x[0][1];
      res.x[1][0] = -idet*x[1][0];
      res.x[1][1] =  idet*x[0][0];
    }
    return res;
  }


  Matrix2 operator*(const Matrix2 &m) {
    Matrix2 res;
    for (int j=0;j<2;j++)
      for (int i=0;i<2;i++)
        res.x[j][i] = x[j][0]*m.x[0][i] + x[j][1]*m.x[1][i];
    return res;
  }


  Vector2 operator*(const Vector2 &v) {
    Vector2 res;
    for (int i=0;i<2;i++)
      res.x[i] = x[i][0]*v.x[0] + x[i][1]*v.x[1];
    return res;
  }


  Matrix2 operator*(double v) {
    Matrix2 res;
    for (int i=0;i<2;i++)
      res.x[i] = x[i]*v;
    return res;
  }


  Matrix2 &operator=(double v) {
    for (int i=0;i<2;i++)
      x[i] = v;
    return *this;
  }


  Matrix2 &operator*=(double v) {
    for (int i=0;i<2;i++)
      x[i] *= v;
    return *this;
  }
  

  Matrix2 &operator+=(const Matrix2 &m) {
    for (int i=0;i<2;i++)
      x[i] += m[i];
    return *this;
  }


  Matrix2 &operator+=(double v) {
    for (int i=0;i<2;i++)
      x[i] += v;
    return *this;
  }


  friend std::ostream &operator<<(std::ostream &os, const Matrix2 &m) {
    os << "[ " << m[0][0] << " " << m[0][1] << " ;";
    os << "  " << m[1][0] << " " << m[1][1] << " ]";
    return os;
  }
}; 

class Matrix3 {
public:
  Vector3 x[3];
  
  Matrix3() {  
    identity();
  }


  Matrix3(const Matrix3 &m) {
    for (int i=0;i<3;i++) 
      x[i] = m.x[i];
  }
  

  Matrix3 &operator=(const Matrix3 &m) {
    for (int i=0;i<3;i++) 
      x[i] = m.x[i];
    return *this;
  }

  void identity() {
    x[0][0] = x[1][1] = x[2][2] = 1.0;
    x[0][1] = x[0][2] = x[1][0] = 0.0;
    x[1][0] = x[2][0] = x[2][1] = 0.0;
  }
  

  inline Vector3 &operator[](int idx) {
    return x[idx];
  } 


  inline Vector3 operator[](int idx) const {
    return x[idx];
  } 
  

  inline double &operator()(int r, int c) {
    return x[r][c];
  }


  double determinant() {
    return (x[0][0]*(x[1][1]*x[2][2] - x[2][1]*x[1][2]) - 
            x[1][0]*(x[2][1]*x[0][2] - x[0][1]*x[2][2]) +
            x[2][0]*(x[0][1]*x[1][2] - x[1][1]*x[0][2]));         
  }


  Matrix3 invert() { 
    double m00 = x[1][1]*x[2][2] - x[2][1]*x[1][2];
    double m01 = x[2][1]*x[0][2] - x[0][1]*x[2][2];
    double m02 = x[0][1]*x[1][2] - x[1][1]*x[0][2];
    double det = x[0][0]*m00 + x[1][0]*m01 + x[2][0]*m02;
    Matrix3 res;
    if (det!=0.0) {
      double idet = 1.0/det;
      res.x[0][0] = idet*m00;
      res.x[0][1] = idet*m01;
      res.x[0][2] = idet*m02;         
      res.x[1][0] = idet*(x[2][0]*x[1][2] - x[1][0]*x[2][2]);
      res.x[1][1] = idet*(x[0][0]*x[2][2] - x[2][0]*x[0][2]);        
      res.x[1][2] = idet*(x[1][0]*x[0][2] - x[0][0]*x[1][2]);
      res.x[2][0] = idet*(x[1][0]*x[2][1] - x[2][0]*x[1][1]);
      res.x[2][1] = idet*(x[2][0]*x[0][1] - x[0][0]*x[2][1]);
      res.x[2][2] = idet*(x[0][0]*x[1][1] - x[1][0]*x[0][1]);         
    }
    return res;
  }
  

  Matrix3 operator*(const Matrix3 &m) {
    Matrix3 res;
    for (int j=0;j<3;j++)
      for (int i=0;i<3;i++)
        res.x[j][i] = x[j][0]*m.x[0][i] + x[j][1]*m.x[1][i] + x[j][2]*m.x[2][i];
    return res;
  }
  

  Vector3 operator*(const Vector3 &v) {
    Vector3 res;
    for (int i=0;i<3;i++)
      res.x[i] = x[i][0]*v.x[0] + x[i][1]*v.x[1] + x[i][2]*v.x[2];
    return res;
  }
  

  Matrix3 operator*(double v) {
    Matrix3 res;
    for (int i=0;i<3;i++)
      res.x[i] = x[i]*v;
    return res;
  }
  

  Matrix3 &operator=(double v) {
    for (int i=0;i<3;i++)
      x[i] = v;
    return *this;
  }
  

  Matrix3 &operator*=(double v) {
    for (int i=0;i<3;i++)
      x[i] *= v;
    return *this;
  }


  Matrix3 &operator+=(const Matrix3 &m) {
    for (int i=0;i<3;i++)
      x[i] += m[i];
    return *this;
  }
  

  Matrix3 &operator+=(double v) {
    for (int i=0;i<3;i++)
      x[i] += v;
    return *this;
  }


  Matrix2 submatrix(int r, int c) {
    Matrix2 res;
    res.x[0][0] = x[r][c];
    res.x[0][1] = x[r][c+1];
    res.x[1][0] = x[r+1][c];
    res.x[1][1] = x[r+1][c+1];
    return res;
  }



  friend std::ostream &operator<<(std::ostream &os, const Matrix3 &m) {
    os << "[ " << m[0][0] << " " << m[0][1] << " " << m[0][2] << " ;";
    os << "  " << m[1][0] << " " << m[1][1] << " " << m[1][2] << " ;";
    os << "  " << m[2][0] << " " << m[2][1] << " " << m[2][2] << " ]";
    return os;
  }
}; 

#endif