The C++ Matrix3x3 Class

Overview

Affine transformations in n-dimensional space

An affine transformation is defined as a mapping of a point in an n-dimensional affine space to another point in the same space. Affine transformations include rotation, scale, mirror, and translation.

If P is an n-dimensional point in an affine space, then we say that X(P) is the affine transformation of the point. Affine transformations can also be applied to vectors. If v is such a vector, then we say that X(v) is the affine transformation of the vector.

An affine transformation can be represented as (M, t) where M is an n x n matrix, and t is an n-dimensional vector in the associated n-dimensional vector space. Assuming this representation, the transformation is applied to points as:

Similarly, it is applied to vectors as:

By embedding our geometry in an (n+1)-dimensional projective space, we can represent the affine transformation as an (n+1)x(n+1) matrix in which M is the upper left nxn portion, the components of t comprise the first n elements of the last column, and the bottom row is (0, ..., 0, 1).

Affine transformations in 3-dimensional space

We restrict our attention to 3D points and vectors by fixing n=3. In terms of the application of an affine transformation as introduced above (i.e., X(P) =MP + t), instances of class Matrix3x3 are used to represent M, and instances of class AffVector are used to represent t.

As noted above, affine transformations on points in a 3D space can be represented using a 4x4 matrix operating on the projective space equivalents of those points. An instance of class Matrix4x4 can be created for this purpose. The normal mode of operation is to create the 3x3 matrix M first using one of the methods described here, then create the 4x4 matrix using a constructor which takes the 3x3 matrix and one of the following: (i) a fixed point, (ii) a point whose pre- and post-images are known, or (iii) an explicit translation vector, t. Among other things, such a 4x4 matrix may then be passed into a graphics API like OpenGL which assumes all transformations are represented as 4x4 matrices.

The methods and overloaded arithmetic operators specified here define common ways to construct 3x3 matrices encoding affine transformations and to apply them to 3D points and vectors. While nothing would prevent it, this class was not designed to provide facilities for representing 2D transformations on points embedded in a 3D projective space.

Constructors

• Matrix3x3();
  Constructs a 3x3 identity matrix.
• Matrix3x3(const Matrix3x3& M);
  The copy constructor.
• Matrix3x3(double m11, double m12, double m13, double m21, double m22, double m23, double m31, double m32, double m33);
  Builds a matrix from the specified matrix elements.

As reviewed above, a 3x3 matrix by itself only describes the effect of an affine transformation on vectors. (Although the complete affine transformation for points is usually easily constructed from such a matrix. See the references for details.) Therefore the remaining constructors here are described in terms of specific types of transformations on vectors. The references provide pictures and other relevant explanatory material.

• Mirror matrix construction
  • Matrix3x3(const AffVector& mirrorPlaneNormal);
    Builds a matrix which will mirror an arbitrary vector about a plane perpendicular to the vector mirrorPlaneNormal. If mirrorPlaneNormal is a zero vector, the matrix will be constructed as an identity matrix.
  • Matrix3x3(char axis);
    Builds a matrix which will mirror an arbitrary vector about a plane perpendicular to the indicated axis. The parameter axis must be one of 'x', 'X', 'y', 'Y', 'z', or 'Z'. Any other value will result in the matrix being constructed as an identity matrix.
• Rotation matrix construction (Angle is in radians. Positive or negative values are allowed. The positive direction of rotation is determined by the right-hand rule as applied to the rotation axis vector.)
  • Matrix3x3(const AffVector& rotationAxis, double angle);
    Builds a matrix which will rotate an arbitrary vector about an axis parallel to the vector rotationAxis. If rotationAxis is a zero vector, the matrix will be constructed as an identity matrix.
  • Matrix3x3(char axis, double angle);
    Builds a matrix which will rotate an arbitrary vector by the given angle about an axis parallel to the indicated axis. The parameter axis must be one of 'x', 'X', 'y', 'Y', 'z', or 'Z'. Any other value will result in the matrix being constructed as an identity matrix.
• Scale matrix construction
  • Matrix3x3(double sx, double sy, double sz);
    Builds a matrix which will scale a vector about the origin by the three scale factors given. No error conditions are possible since any or all of the scale factors are allowed to be positive, negative, or zero.
• Tensor product matrix construction
  • Matrix3x3(const AffVector& u, const AffVector& v);
    Builds the matrix which is the tensor product of the two given vectors. The (i,j)-th element of the matrix is u[i]*v[j]. No error conditions are possible. If either u or v is a zero vector, then the resulting matrix will simply be a zero matrix.

C++ Overloaded Operators

In the following, assume that instances of Matrix3x3 called M, N, S1, and S2 have been declared and initialized. Similarly, let us assume that P and Q have been defined as instances of class AffPoint and that u and v have been defined as instances of class AffVector.

• operators implemented as instance methods
  • Matrix3x3 operator=(const Matrix3x3& RHS);
    The normal assignment operator. For example: M = N = S1 + S2;
  • bool operator==(const Matrix3x3& RHS) const;
    The matrix comparison function. Returns true if and only if each corresponding element is equal within BasicDistanceTol.
  • Matrix3x3 operator*=(const Matrix3x3& RHS);
    Analogous to *= for simple numeric types. For example: "M *= N;" performs the 3x3 matrix multiplication "M * N", and the result replaces M.
  • Matrix3x3 operator*=(double f);
    Each element of this matrix is scaled by the scale factor, f. For example: M *= 2.3;.
  • Matrix3x3 operator+=(const Matrix3x3& RHS);
    Corresponding elements of this matrix and the RHS matrix are added. The resulting matrix replaces this matrix. For example: M += N;
  • AffPoint operator*(const AffPoint& P) const;
    The operation "P = M * Q;" multiples the 3D point Q by M, storing the resulting transformed point in P.
  • AffVector operator*(const AffVector& v) const;
    The operation "v = M * u;" multiples the 3D vector u by M, storing the resulting transformed vector in v.
• operators implemented as external functions
  • Matrix3x3 operator*(const Matrix3x3& M1, const Matrix3x3& M2);
    The two indicated 3x3 matrices are multiplied together, and the product is returned as the function return value. For example: Matrix3x3 prod = M * N;.
  • Matrix3x3 operator+(const Matrix3x3& M1, const Matrix3x3& M2);
    The two indicated 3x3 matrices are added element by element. The resulting matrix is returned as the function return value. For example: Matrix3x3 sum = M + N;.
  • Matrix3x3 operator-(const Matrix3x3& M1, const Matrix3x3& M2);
    The two indicated 3x3 matrices are subtracted element by element. The resulting matrix is returned as the function return value. For example: Matrix3x3 diff = M - N;.
  • Matrix3x3 operator*(double f, const Matrix3x3& M);
    Each element of the matrix M is multiplied by f, and the resulting matrix is returned as the function return value. For example: Matrix3x3 scaledM = 10.3 * M;.
  • ostream& operator<<(ostream& os, const Matrix3x3& m);
    The matrix is written to the given output stream. (See the class methods outputFormat and ioFormat below for details.) Example: cout << M;
  • istream& operator>>(istream& is, Matrix3x3& m);
    The matrix is read from the given input stream. (See the class methods inputFormat and ioFormat below for details.) Example: cin >> M;

Matrix3x3 Instance Methods

In the following, assume that an instance of Matrix3x3 called M has been declared and initialized.

• double determinant() const;
  Computes and returns the determinant of the matrix. For example: double det = M.determinant();
• double elementAt(int i, int j) const;
  This routine returns the matrix element at (i,j). For example: double elem = M.elementAt(2,1); The results are unpredictable if either i or j is outside the range 0..2.
• int extractAxisAngle(AffVector& w, double& theta) const;
  A 3x3 matrix which is orthogonal and right-handed describes a rigid transformation on vectors. Any such transformation is equivalent to a single rotation about some axis. This routine verifies that this matrix is orthogonal and right-handed (see methods below). If it is, it computes and returns in w and theta the axis and angle which describe the equivalent rotation. If this matrix is either not orthogonal or not right-handed, the values of w and theta are undefined, and one of the following two error messages will be returned as the function return value: Matrix3x3::NotOrthogonal, Matrix3x3::NotRightHanded. If w and theta are successfully determined, the function return value will be Matrix3x3::Extracted_wTheta. Any other function return value corresponds to an unanticipated internal error, and w and theta should be considered undefined. For example:

      AffVector equivAxis;
      double    rotAngle;
      int     result = M.extractAxisAngle(equivAxis, rotAngle);
      if (result == Matrix3x3::Extracted_wTheta)
      {
          // code which uses 'equivAxis' and/or 'rotAngle'
      }
      else
      {
          // error handling code....
      }
  

• void extractRows(AffVector& row1, AffVector& row2, AffVector& row3) const;
  This routine extracts the three rows of a matrix, placing each in the indicated output vectors. One could, for example, use this routine while checking whether a matrix is right-handed:

      AffVector r1, r2, r3;
      M.extractRows(r1,r2,r3);
      if (AffVector::dot(r1.cross(r2),r3))
      {
          // the matrix is right-handed....
      }
      else
      {
          // it isn't....
      }
  

• bool isOrthogonal() const;
  This routine returns true if this matrix is orthogonal. For example:

      if (M.isOrthogonal())
          cout << M;
  

• bool isRightHanded() const;
  This routine returns true if this matrix is right-handed. For example:

      if (M.isRightHanded())
          cout << M;
  

• double largestDiagonalElement(int& pos) const;
  This routine will return the position of the largest element along the diagonal of the matrix. For example:

      int  pos = -1;
      M.largestDiagonalElement(pos);
      cout << "Largest diagonal element is: " << M.elementAt(pos,pos) << '\n';
  

• void multiply(const double a[ ], double b[ ]) const;
  This routine multiplies the array in 'a' (which is assumed to have length 3) by the matrix, placing the result in the output array 'b' (which is also assumed to be of length 3). Do not pass the same array as the actual parameter for 'a' and 'b'.
• double setElementAt(int i, int j, double newValue);
  This routine sets the matrix element at (i,j) to newValue. For example: M.setElementAt(2,1,-10.329);

  The matrix is unchanged if either i or j is outside the range 0..2.

• double trace() const;
  This routine computes and returns the trace of the matrix. The trace is the sum of the diagonal elements. For example: double t = M.trace();

Matrix3x3 Public Constants

The following constants are Matrix3x3 class variables that can be used in your programs.

• Special Matrices
  • const Matrix3x3 IdentityMatrix;
    The identity matrix.
  • const Matrix3x3 ZeroMatrix;
    The zero matrix.

For example, the mirror matrix constructor described above does in part the equivalent of the following:

		AffVector n(1.0, 2.0, -1.3);
		n.normalize();
		Matrix3x3 N = Matrix3x3::tensorProductMatrix(n,n);
		Matrix3x3 Mirror = Matrix3x3::Identity - 2.0 * N;

• Values for the parameter ioAspect in the class methods ioFormat, inputFormat, and outputFormat
  • const int elementFieldWidth;
  • const int openMatrixDelimiter;
  • const int openRowDelimiter;
  • const int elementSeparator;
  • const int rowSeparator;

See the description of the indicated class methods for instructions and examples on the use of these constants.

Matrix3x3 Class Methods

Most of these class methods provide alternative ways to construct 3x3 matrices. If any of these methods cannot create the matrix due to some error in supplied parameters (e.g., a zero-length vector was provided), an identity matrix will be returned. Angles are assumed to be measured in radians.

• Matrix3x3 crossProductMatrix(const AffVector& u);
  This class method computes and returns a matrix, U such that, for any arbitrary vector, v:
  u x v = Uv

• char* getInputFormat(int iAspect);
  
• char* getOutputFormat(int oAspect);
  
These methods return the character strings currently being used as indicated. (See the descriptions of inputFormat, ioFormat, and outputFormat below.) A copy dynamically allocated by copyString (defined in Basic.h) is returned. It is the responsibility of the caller to use delete [] when the string is no longer needed.
• int getOutputElementFieldWidth();
  
This method returns the field width currently being used to output matrix elements.
• void inputFormat(int iAspect, char iValue);
• void inputFormat(int iAspect, const char* iValue);
• void ioFormat(int ioAspect, char ioValue);
• void ioFormat(int ioAspect, const char* ioValue);
• void outputFormat(int oAspect, char oValue);
• void outputFormat(int oAspect, const char* oValue);
• void outputFormat(int oAspect, int oValue);
  Use of these class methods allows you to control whether the >> and << operators defined above use parentheses, commas, and so forth. By default, the input operator assumes that no opening delimiters or separators are used; the output operator assumes that braces are used as opening delimiters, and blanks are used as element separators. Consider the following code:

        cin >> M;
        cout << "M = " << M << '\n';
  

  Assuming the default conditions specified, if the input is:

  13.2 -1.2 -0.334 1.9 21.5 97 1.2 18.2 0.32
  

  then the output will be:

  M = {{        13.2         -1.2       -0.334} {         1.9         21.5           97} {         1.2         18.2         0.32}}
  

  The methods inputFormat and outputFormat allow you to change these assumptions for the input and output operators, respectively. The ioFormat method is equivalent to invoking both inputFormat and outputFormat with the given parameters.

  The valid values for "ioAspect" are listed in the "public constants" section above. Using these, we can format the preceding more reasonably as:

      cin >> M;
  
      Matrix3x3::outputFormat(Matrix3x3::elementFieldWidth, 5);
      Matrix3x3::outputFormat(Matrix3x3::openMatrixDelimiter, "{\n");
      Matrix3x3::outputFormat(Matrix3x3::rowSeparator, ",\n");
  
      cout << "M = " << M << '\n';
  

  If the input is the same as indicated above, the output will be:

  M = {
  {   13.2    -1.2  -0.334},
  {    1.9    21.5      97},
  {    1.2    18.2    0.32}
  }
  

  Any character string can be specified for any of these formatting options (subject to a maximum length). For the "open" strings, a matching closing character string will be generated by matching each input character with an "inverse". If, for example, the string "(*\n" is specified as the openDelimiter, then "\n*)" is generated as the closing delimiter.

  The input operator does not actually check for an exact match on delimiters and separators. Instead it knows the number of nonblank characters in each, and it simply skips that many nonblank characters on input.

• Matrix3x3 mirrorMatrix(const AffVector& mirrorPlaneNormal);
  This class method creates and returns a 3x3 matrix which will mirror vectors about a plane perpendicular to the given mirror plane normal vector. If the given vector has zero length, a 3x3 identity matrix is returned.
• Matrix3x3 rotationMatrix(const AffVector& rotationAxis, double angle);
  This class method creates and returns a 3x3 matrix which will rotate vectors about an axis parallel to the given rotation axis vector. If the given vector has zero length, a 3x3 identity matrix is returned.
• Matrix3x3 rotationMatrix(const AffVector& uFrom, const AffVector& uTo);
  This class method creates and returns a 3x3 matrix which will rotate the "from" vector onto the "to" vector. If either vector has zero length, a 3x3 identity matrix is returned.
• Matrix3x3 rotationMatrix(const AffVector& uFrom, const AffVector& vFrom, const AffVector& uTo, const AffVector& vTo);
  This class method creates and returns a 3x3 matrix which rotates the orientation vectors of one local coordinate system onto those of another. The orientation of the two coordinate systems is defined by a (u, v) vector pair. Specifically, the u vector defines the direction of a local x axis, and the component of v perpendicular to u defines the direction of a local y axis. The matrix returned rotates the "from" coordinate system directions onto those of the "to" coordinate system. More precisely, it is equivalent to first rotating the "uFrom" vector onto the "uTo" vector, and then rotating about this common "u" direction until (i) "vFrom" is in the plane of "uTo-vTo", and (ii) the component of "vFrom" perpendicular to "uTo" is parallel to the component of "vTo" perpendicular to "uTo". If it is not possible to create this matrix for any reason (e.g., if any of the four vectors has zero length or if the component of v perpendicular to u is a zero vector for either pair), then an identity matrix is returned.
• Matrix3x3 scaleMatrix(double sx, double sy, double sz);
  This class method returns a 3x3 matrix which will scale vectors by the indicated scale factors in the three coordinate directions. No error conditions are possible; scale factors can be positive, negative, or zero.
• Matrix3x3 shearMatrix(const AffVector& n, const AffVector& u, double f);
  This class method returns a 3x3 matrix which will shear vectors by an amount that is proportional to the length of their component parallel to n. Imagine that n defines the orientation of a "shear plane". Define n_hat as a unit vector in the direction of n, and define u_shear as a unit vector in the direction of the component of u perpendicular to n_hat. The u_shear vector defines the shear direction in the shear plane. Vectors in the shear plane are not modified at all. An arbitrary vector v not in the shear plane is sheared in the +u_shear or -u_shear direction, according to the side of the plane to which it points. The amount of the shear is proportional to the length of v parallel to n. The parameter f is the constant of proportionality. Specifically, the matrix returned will transform v into v' as follows: v' = v + f * (v ^. n_hat)u_shear If any of the vectors mentioned in the paragraph above have zero length and/or if f=0, an identity matrix is returned.
• Matrix3x3 tensorProductMatrix(const AffVector& u, const AffVector& v);
  This class method returns the 3x3 matrix determined by the tensor product of the two given vectors.
• Matrix3x3 xRotationMatrix(double angle);
• Matrix3x3 yRotationMatrix(double angle);
• Matrix3x3 zRotationMatrix(double angle);
  These three class methods return a 3x3 matrix which will rotate vectors about the indicated axis by the given angle in radians.