Added LinearAlgebra class with vector/matrix operations. References #3969. Fixes #4194

#ifndef LinearAlgebra_h

#define LinearAlgebra_h

/**

 * @file

 * $Revision: 1.6 $

 * $Date: 2009/12/22 02:09:54 $

 *

 *   Unless noted otherwise, the portions of Isis written by the USGS are public

 *   domain. See individual third-party library and package descriptions for

 *   intellectual property information,user agreements, and related information.

 *

 *   Although Isis has been used by the USGS, no warranty, expressed or implied,

 *   is made by the USGS as to the accuracy and functioning of such software

 *   and related material nor shall the fact of distribution constitute any such

 *   warranty, and no responsibility is assumed by the USGS in connection

 *   therewith.

 *

 *   For additional information, launch

 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html in a browser or see

 *   the Privacy & Disclaimers page on the Isis website,

 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on

 *   http://www.usgs.gov/privacy.html.

 */

// std library

#include <iostream>

// boost library

#include <boost/numeric/ublas/matrix.hpp>

#include <boost/numeric/ublas/io.hpp>

// Qt Library

#include <QDebug>

#include <QPair>

namespace Isis {

  class Angle;

  class Distance;

  class Latitude;

  class Longitude;

  /**

   * This class holds all static methods to perform linear algebra 

   * operations on vectors and matrices. 

   *  

   * This class has 4 representations of a 3 dimensional rotation 

   * <ul> 

   *   <li> A 3 X 3 rotation matrix </li>

   *   <li> An angle/axis pair - This pair is made up of a vector (axis)

   *   and angle of rotation about that axis.</li>

   *   <li> A quaternion - This is saved as a 4-valued unit vector of the

   *   form q = (q0, q1, q2, q3) where for the angle/axis pair of the form

   *   (A, (ux, uy, uz)), we have q0 = cos(A/2), q1 = ux*sin(A/2), q2 =

   *   uy*sin(A/2), q3 = uz*sin(A/2).</li>

   *   <li> A set of Euler angles - This is a set of three angles and

   *   their corresponding axis of rotation. In this case, each axis

   *   must be one of the x, y, or z axes for the coordinate system.</li>

   * </ul>

   *  

   * @author 2013-06-01 Jeannie Backer

   *

   * @internal

   *   @history 2013-06-01 Jeannie Backer - Original version.

   *   @history 2013-08-07 Kimberly Oyama - Updated documentation.

   *   @history 2016-07-25 Jeannie Backer - Updated documentation and test.

   *  

   *  

   *   @todo document methods (a) add naif routine names to documentation where appropriate,

   *                          (b) clean up comments within methods,

   *                          (3) use naif documentation where appropriate

   *   @todo add links to naif documentation web pages???

   *   @todo add proper scopes

   *   @todo implement norm2 (L2 norm)?

   *   @todo implement trace? 

   *   @todo common validate errors - vector size match, 3x3 matrix, 3D vector 

   *   @todo what is more expensive - unitize vector and compare to 1 or isOrthogonal?

   *   @todo rename local variables m2q and q2m

   *   @todo implement qdq2av (Quaternion and quaternion derivative to angular velocity)?

   *   @todo implement qxq (Quaternion times quaternion)?

   *   @todo implement qderiv (Quadratic derivative)?

   *   @todo When converting quaternion, axis-angle, Euler angle to rot 

   *         matrix, should we verify that the matrix is valid (is a rotation matrix)???

   *   @todo isQuaternion??? 

   *   @todo should we throw error in q2m if q = 0???

   *   @todo q2m - why check if sharpen when q is not unit if it is required

   *         to be a unit? throw an error?

   *   @todo axisar - throw error if axis is origin? If we rotate about 

   *         about origin, then we get a zero matrix - not a rotation...

   *         except maybe trivially... (or do we get identity??? indicating no rotation)

   *   @todo double check all indexing on code copied from NAIF - indices 

   *         should reverse col/row order and subtract 1

   */

  class LinearAlgebra {

    public:

      /**

       * Definition for an Isis::LinearAlgebra::Matrix of doubles. This is a 

       * typedef for a boost matrix. 

       *  

       * Note: This typedef is used so that we can add functionality to an 

       * existing matrix type and/or change which third party library's matrix

       * we are using without changing all references to this type in the ISIS

       * API. 

       */

      typedef boost::numeric::ublas::matrix<double> Matrix;

      /**

       * Definition for an Isis::LinearAlgebra::Vector of doubles. This is a 

       * typedef for a boost vector. 

       *  

       * Note: This typedef is used so that we can add functionality to an 

       * existing vector type and/or change which third party library's vector

       * we are using without changing all references to this type in the ISIS

       * API. 

       */

      typedef boost::numeric::ublas::vector<double> Vector;

      // define AxisAngle and EulerAngle

      /**

       * Definition for an Axis-Angle pair. This is a three dimensional rotation

       * represented as an axis of rotation and a corresponding rotation angle. 

       * AxisAngle is a typedef for a QPair of an Isis::LinearAlgebra::Vector 

       * and an Isis::Angle. 

       */

      typedef QPair<Vector, Angle> AxisAngle;

      /**

       * Definition for an EulerAngle pair. This is a three dimensional rotation

       * represented as an Euler angle and the number corresponding to its 

       * rotation axis. EulerAngle is a typedef for a QPair of an Isis::Angle 

       * and an integer (1, 2, or 3).

       */

      typedef QPair<Angle, int> EulerAngle;

      // check type of matrix/vector

      static bool isIdentity(const Matrix &matrix);

      static bool isOrthogonal(const Matrix &matrix);

      static bool isRotationMatrix(const Matrix &matrix);

      static bool isZero(const Matrix &matrix);

      static bool isZero(const Vector &vector);

      static bool isEmpty(const Vector &vector);

      static bool isUnit(const Vector &vector);

      // create special matrices

      static Matrix identity(int size);

      static Matrix transpose(const Matrix &matrix);

      static Matrix inverse(const Matrix &matrix);

      static Matrix zeroMatrix(int rows, int columns);

      static Vector zeroVector(int size);

      // 

      static double determinant(const Matrix &matrix);

      // magnitude/norm based calculations

      static Vector normalize(const Vector &vector);

      static double magnitude(const Vector &vector);

      static double absoluteMaximum(const Vector &vector);

      // arithmetic operations

      static Matrix multiply(const Matrix &matrix1, const Matrix &matrix2);

      static Vector multiply(const Matrix &matrix, const Vector &vector);

      static Vector multiply(double scalar, const Vector &vector);

      static Matrix multiply(double scalar, const Matrix &matrix);

      static Vector add(const Vector &vector1, const Vector &vector2);

      static Vector subtract(const Vector &vector1, const Vector &vector2);

      // vector products

      static Vector crossProduct(const Vector &vector1, const Vector &vector2);

      static Vector normalizedCrossProduct(const Vector &vector1, const Vector &vector2);

      static Matrix outerProduct(const Vector &vector1, const Vector &vector2);

      static double dotProduct(const Vector &vector1, const Vector &vector2);

      static double innerProduct(const Vector &vector1, const Vector &vector2);

      // vector movements (projections, rotations)

      static Vector project(const Vector &vector1, const Vector &vector2);

      static Vector rotate(const Vector &vector, const Vector &axis, Angle angle);

      // 

      static Vector perpendicular(const Vector &vector1, const Vector &vector2);

      // converters

      static AxisAngle toAxisAngle(const Matrix &rotationMatrix);//raxisa

      static Matrix toMatrix(const AxisAngle &axisAngle); // axisar

      static Matrix toMatrix(const Vector &axis, Angle angle); // axisar

      static QList<EulerAngle> toEulerAngles(const Matrix &rotationMatrix, const QList<int> axes);// m2eul

      static Matrix toMatrix(const QList<EulerAngle> &eulerAngles); // eul2m

      static Matrix toMatrix(const EulerAngle &angle3, const EulerAngle &angle2, const EulerAngle &angle1); // eul2m

      static Vector toQuaternion(const Matrix &rotationMatrix);// m2q

      static Matrix toMatrix(const Vector &quaternion);// q2m

      //Matrix row/col to vector and vector to row/columns

      static void setRow(Matrix &matrix, const Vector &vector, int rowIndex);

      static void setColumn(Matrix &matrix, const Vector &vector, int columnIndex);

      static Vector row(const Matrix &matrix, int rowIndex);

      static Vector column(const Matrix &matrix, int columnIndex);

      static void setVec3(Vector *v, double v0, double v1, double v2);

      static void setVec4(Vector *v, double v0, double v1, double v2, double v3);

      static Vector vector(double v0, double v1, double v2);

      static Vector vector(double v0, double v1, double v2, double v3);

      static Vector subVector(const Vector &v, int start, int size);

      // No need to implement the std::cout operator<< since this is already

      // done in boost/numeric/ublas/io.hpp, which is included above

      // The print format for a size 3 vector is

      // [3] (1, 2, 3)

      // 

      // friend ostream &operator<<(ostream &os, LinearAlgebra::Vector &vector);

    protected:

      LinearAlgebra();

      ~LinearAlgebra();

    private:

  };

  // these must be declared outside of the class (at the end since they must be

  // declared after the typedefs)

  QDebug operator<<(QDebug dbg, const LinearAlgebra::Vector &vector);

  QDebug operator<<(QDebug dbg, const LinearAlgebra::Matrix &matrix);

  QString toString(const LinearAlgebra::Vector &vector, int precision=15);

};

#endif

ifeq ($(ISISROOT), $(BLANK))

.SILENT:

error:

	echo "Please set ISISROOT";

else

	include $(ISISROOT)/make/isismake.objs

endif