vpHomogeneousMatrix.cpp

/****************************************************************************
 *
 * ViSP, open source Visual Servoing Platform software.
 * Copyright (C) 2005 - 2022 by Inria. All rights reserved.
 *
 * This software is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 * See the file LICENSE.txt at the root directory of this source
 * distribution for additional information about the GNU GPL.
 *
 * For using ViSP with software that can not be combined with the GNU
 * GPL, please contact Inria about acquiring a ViSP Professional
 * Edition License.
 *
 * See http://visp.inria.fr for more information.
 *
 * This software was developed at:
 * Inria Rennes - Bretagne Atlantique
 * Campus Universitaire de Beaulieu
 * 35042 Rennes Cedex
 * France
 *
 * If you have questions regarding the use of this file, please contact
 * Inria at visp@inria.fr
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 * Description:
 * Homogeneous matrix.
 *
 * Authors:
 * Eric Marchand
 *
 *****************************************************************************/
 
#include <visp3/core/vpDebug.h>
#include <visp3/core/vpException.h>
#include <visp3/core/vpHomogeneousMatrix.h>
#include <visp3/core/vpMatrix.h>
#include <visp3/core/vpPoint.h>
#include <visp3/core/vpQuaternionVector.h>
 
vpHomogeneousMatrix::vpHomogeneousMatrix(const vpTranslationVector &t, const vpQuaternionVector &q)
  : vpArray2D<double>(4, 4)
{
  buildFrom(t, q);
  (*this)[3][3] = 1.;
}
 
vpHomogeneousMatrix::vpHomogeneousMatrix() : vpArray2D<double>(4, 4), m_index(0) { eye(); }
 
vpHomogeneousMatrix::vpHomogeneousMatrix(const vpHomogeneousMatrix &M) : vpArray2D<double>(4, 4), m_index(0)
{
  *this = M;
}
 
vpHomogeneousMatrix::vpHomogeneousMatrix(const vpTranslationVector &t, const vpThetaUVector &tu)
  : vpArray2D<double>(4, 4), m_index(0)
{
  buildFrom(t, tu);
  (*this)[3][3] = 1.;
}
 
vpHomogeneousMatrix::vpHomogeneousMatrix(const vpTranslationVector &t, const vpRotationMatrix &R)
  : vpArray2D<double>(4, 4), m_index(0)
{
  insert(R);
  insert(t);
  (*this)[3][3] = 1.;
}
 
vpHomogeneousMatrix::vpHomogeneousMatrix(const vpPoseVector &p) : vpArray2D<double>(4, 4), m_index(0)
{
  buildFrom(p[0], p[1], p[2], p[3], p[4], p[5]);
  (*this)[3][3] = 1.;
}
 
vpHomogeneousMatrix::vpHomogeneousMatrix(const std::vector<float> &v) : vpArray2D<double>(4, 4), m_index(0)
{
  buildFrom(v);
  (*this)[3][3] = 1.;
}
 
#if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11)
vpHomogeneousMatrix::vpHomogeneousMatrix(const std::initializer_list<double> &list)
  : vpArray2D<double>(4, 4), m_index(0)
{
  if (list.size() == 12) {
    std::copy(list.begin(), list.end(), data);
    data[12] = 0.;
    data[13] = 0.;
    data[14] = 0.;
    data[15] = 1.;
  } else if (list.size() == 16) {
    std::copy(list.begin(), list.end(), data);
    for (size_t i = 12; i < 15; i++) {
      if (std::fabs(data[i]) > std::numeric_limits<double>::epsilon()) {
        throw(vpException(vpException::fatalError,
                          "Cannot initialize homogeneous matrix. "
                          "List element %d (%f) should be 0.",
                          i, data[i]));
      }
    }
    if (std::fabs(data[15] - 1.) > std::numeric_limits<double>::epsilon()) {
      throw(vpException(vpException::fatalError,
                        "Cannot initialize homogeneous matrix. "
                        "List element 15 (%f) should be 1.",
                        data[15]));
    }
  } else {
    throw(vpException(vpException::fatalError,
                      "Cannot initialize homogeneous matrix from a list (%d elements) that has not 12 or 16 elements",
                      list.size()));
  }
 
  if (!isAnHomogeneousMatrix()) {
    if (isAnHomogeneousMatrix(1e-3)) {
      // re-orthogonalize rotation matrix since the input is close to a valid rotation matrix
      vpRotationMatrix R(*this);
      R.orthogonalize();
 
      data[0] = R[0][0];
      data[1] = R[0][1];
      data[2] = R[0][2];
      data[4] = R[1][0];
      data[5] = R[1][1];
      data[6] = R[1][2];
      data[8] = R[2][0];
      data[9] = R[2][1];
      data[10] = R[2][2];
    } else {
      throw(vpException(
          vpException::fatalError,
          "Homogeneous matrix initialization fails since its elements are not valid (rotation part or last row)"));
    }
  }
}
#endif
 
vpHomogeneousMatrix::vpHomogeneousMatrix(const std::vector<double> &v) : vpArray2D<double>(4, 4), m_index(0)
{
  buildFrom(v);
  (*this)[3][3] = 1.;
}
 
vpHomogeneousMatrix::vpHomogeneousMatrix(double tx, double ty, double tz, double tux, double tuy, double tuz)
  : vpArray2D<double>(4, 4), m_index(0)
{
  buildFrom(tx, ty, tz, tux, tuy, tuz);
  (*this)[3][3] = 1.;
}
 
void vpHomogeneousMatrix::buildFrom(const vpTranslationVector &t, const vpThetaUVector &tu)
{
  insert(tu);
  insert(t);
}
 
void vpHomogeneousMatrix::buildFrom(const vpTranslationVector &t, const vpRotationMatrix &R)
{
  insert(R);
  insert(t);
}
 
void vpHomogeneousMatrix::buildFrom(const vpPoseVector &p)
{
  vpTranslationVector tv(p[0], p[1], p[2]);
  vpThetaUVector tu(p[3], p[4], p[5]);
 
  insert(tu);
  insert(tv);
}
 
void vpHomogeneousMatrix::buildFrom(const vpTranslationVector &t, const vpQuaternionVector &q)
{
  insert(t);
  insert(q);
}
 
void vpHomogeneousMatrix::buildFrom(double tx, double ty, double tz, double tux, double tuy, double tuz)
{
  vpRotationMatrix R(tux, tuy, tuz);
  vpTranslationVector t(tx, ty, tz);
 
  insert(R);
  insert(t);
}
 
void vpHomogeneousMatrix::buildFrom(const std::vector<float> &v)
{
  if (v.size() != 12 && v.size() != 16) {
    throw(vpException(vpException::dimensionError, "Cannot convert std::vector<float> to vpHomogeneousMatrix"));
  }
 
  for (unsigned int i = 0; i < 12; i++)
    this->data[i] = (double)v[i];
}
 
void vpHomogeneousMatrix::buildFrom(const std::vector<double> &v)
{
  if (v.size() != 12 && v.size() != 16) {
    throw(vpException(vpException::dimensionError, "Cannot convert std::vector<double> to vpHomogeneousMatrix"));
  }
 
  for (unsigned int i = 0; i < 12; i++)
    this->data[i] = v[i];
}
 
vpHomogeneousMatrix &vpHomogeneousMatrix::operator=(const vpHomogeneousMatrix &M)
{
  for (int i = 0; i < 4; i++) {
    for (int j = 0; j < 4; j++) {
      rowPtrs[i][j] = M.rowPtrs[i][j];
    }
  }
  return *this;
}
 
vpHomogeneousMatrix vpHomogeneousMatrix::operator*(const vpHomogeneousMatrix &M) const
{
  vpHomogeneousMatrix p;
 
  vpRotationMatrix R1, R2, R;
  vpTranslationVector T1, T2, T;
 
  extract(T1);
  M.extract(T2);
 
  extract(R1);
  M.extract(R2);
 
  R = R1 * R2;
 
  T = R1 * T2 + T1;
 
  p.insert(T);
  p.insert(R);
 
  return p;
}
 
vpHomogeneousMatrix &vpHomogeneousMatrix::operator*=(const vpHomogeneousMatrix &M)
{
  (*this) = (*this) * M;
  return (*this);
}
 
vpColVector vpHomogeneousMatrix::operator*(const vpColVector &v) const
{
  if (v.getRows() != 4) {
    throw(vpException(vpException::dimensionError,
                      "Cannot multiply a (4x4) homogeneous matrix by a "
                      "(%dx1) column vector",
                      v.getRows()));
  }
  vpColVector p(rowNum);
 
  p = 0.0;
 
  for (unsigned int j = 0; j < 4; j++) {
    for (unsigned int i = 0; i < 4; i++) {
      p[i] += rowPtrs[i][j] * v[j];
    }
  }
 
  return p;
}
 
vpPoint vpHomogeneousMatrix::operator*(const vpPoint &bP) const
{
  vpPoint aP;
 
  vpColVector v(4), v1(4);
 
  v[0] = bP.get_X();
  v[1] = bP.get_Y();
  v[2] = bP.get_Z();
  v[3] = bP.get_W();
 
  v1[0] = (*this)[0][0] * v[0] + (*this)[0][1] * v[1] + (*this)[0][2] * v[2] + (*this)[0][3] * v[3];
  v1[1] = (*this)[1][0] * v[0] + (*this)[1][1] * v[1] + (*this)[1][2] * v[2] + (*this)[1][3] * v[3];
  v1[2] = (*this)[2][0] * v[0] + (*this)[2][1] * v[1] + (*this)[2][2] * v[2] + (*this)[2][3] * v[3];
  v1[3] = (*this)[3][0] * v[0] + (*this)[3][1] * v[1] + (*this)[3][2] * v[2] + (*this)[3][3] * v[3];
 
  v1 /= v1[3];
 
  //  v1 = M*v ;
  aP.set_X(v1[0]);
  aP.set_Y(v1[1]);
  aP.set_Z(v1[2]);
  aP.set_W(v1[3]);
 
  aP.set_oX(v1[0]);
  aP.set_oY(v1[1]);
  aP.set_oZ(v1[2]);
  aP.set_oW(v1[3]);
 
  return aP;
}
 
vpTranslationVector vpHomogeneousMatrix::operator*(const vpTranslationVector &t) const
{
  vpTranslationVector t_out;
  t_out[0] = (*this)[0][0] * t[0] + (*this)[0][1] * t[1] + (*this)[0][2] * t[2] + (*this)[0][3];
  t_out[1] = (*this)[1][0] * t[0] + (*this)[1][1] * t[1] + (*this)[1][2] * t[2] + (*this)[1][3];
  t_out[2] = (*this)[2][0] * t[0] + (*this)[2][1] * t[1] + (*this)[2][2] * t[2] + (*this)[2][3];
 
  return t_out;
}
 
vpHomogeneousMatrix &vpHomogeneousMatrix::operator<<(double val)
{
  m_index = 0;
  data[m_index] = val;
  return *this;
}
 
vpHomogeneousMatrix &vpHomogeneousMatrix::operator,(double val)
{
  m_index++;
  if (m_index >= size()) {
    throw(vpException(vpException::dimensionError,
                      "Cannot set homogenous matrix out of bounds. It has only %d elements while you try to initialize "
                      "with %d elements",
                      size(), m_index + 1));
  }
  data[m_index] = val;
  return *this;
}
 
/*********************************************************************/
 
bool vpHomogeneousMatrix::isAnHomogeneousMatrix(double threshold) const
{
  vpRotationMatrix R;
  extract(R);
 
  const double epsilon = std::numeric_limits<double>::epsilon();
  return R.isARotationMatrix(threshold) && vpMath::nul((*this)[3][0], epsilon) && vpMath::nul((*this)[3][1], epsilon) &&
         vpMath::nul((*this)[3][2], epsilon) && vpMath::equal((*this)[3][3], 1.0, epsilon);
}
 
bool vpHomogeneousMatrix::isValid() const
{
  for (unsigned int i = 0; i < size(); i++) {
    if (vpMath::isNaN(data[i])) {
      return false;
    }
  }
  return true;
}
 
void vpHomogeneousMatrix::extract(vpRotationMatrix &R) const
{
  for (unsigned int i = 0; i < 3; i++)
    for (unsigned int j = 0; j < 3; j++)
      R[i][j] = (*this)[i][j];
}
 
void vpHomogeneousMatrix::extract(vpTranslationVector &t) const
{
  t[0] = (*this)[0][3];
  t[1] = (*this)[1][3];
  t[2] = (*this)[2][3];
}
void vpHomogeneousMatrix::extract(vpThetaUVector &tu) const
{
  vpRotationMatrix R;
  (*this).extract(R);
  tu.buildFrom(R);
}
 
void vpHomogeneousMatrix::extract(vpQuaternionVector &q) const
{
  vpRotationMatrix R;
  (*this).extract(R);
  q.buildFrom(R);
}
 
void vpHomogeneousMatrix::insert(const vpRotationMatrix &R)
{
  for (unsigned int i = 0; i < 3; i++)
    for (unsigned int j = 0; j < 3; j++)
      (*this)[i][j] = R[i][j];
}
 
void vpHomogeneousMatrix::insert(const vpThetaUVector &tu)
{
  vpRotationMatrix R(tu);
  insert(R);
}
 
void vpHomogeneousMatrix::insert(const vpTranslationVector &t)
{
  (*this)[0][3] = t[0];
  (*this)[1][3] = t[1];
  (*this)[2][3] = t[2];
}
 
void vpHomogeneousMatrix::insert(const vpQuaternionVector &q) { insert(vpRotationMatrix(q)); }
 
vpHomogeneousMatrix vpHomogeneousMatrix::inverse() const
{
  vpHomogeneousMatrix Mi;
 
  vpRotationMatrix R;
  extract(R);
  vpTranslationVector T;
  extract(T);
 
  vpTranslationVector RtT;
  RtT = -(R.t() * T);
 
  Mi.insert(R.t());
  Mi.insert(RtT);
 
  return Mi;
}
 
void vpHomogeneousMatrix::eye()
{
  (*this)[0][0] = 1;
  (*this)[1][1] = 1;
  (*this)[2][2] = 1;
  (*this)[3][3] = 1;
 
  (*this)[0][1] = (*this)[0][2] = (*this)[0][3] = 0;
  (*this)[1][0] = (*this)[1][2] = (*this)[1][3] = 0;
  (*this)[2][0] = (*this)[2][1] = (*this)[2][3] = 0;
  (*this)[3][0] = (*this)[3][1] = (*this)[3][2] = 0;
}
 
void vpHomogeneousMatrix::inverse(vpHomogeneousMatrix &M) const { M = inverse(); }
 
void vpHomogeneousMatrix::save(std::ofstream &f) const
{
  if (!f.fail()) {
    f << *this;
  } else {
    throw(vpException(vpException::ioError, "Cannot save homogeneous matrix: ostream not open"));
  }
}
 
void vpHomogeneousMatrix::load(std::ifstream &f)
{
  if (!f.fail()) {
    for (unsigned int i = 0; i < 4; i++) {
      for (unsigned int j = 0; j < 4; j++) {
        f >> (*this)[i][j];
      }
    }
  } else {
    throw(vpException(vpException::ioError, "Cannot load homogeneous matrix: ifstream not open"));
  }
}
 
void vpHomogeneousMatrix::orthogonalizeRotation()
{
  vpRotationMatrix R(*this);
  R.orthogonalize();
 
  data[0] = R[0][0];
  data[1] = R[0][1];
  data[2] = R[0][2];
  data[4] = R[1][0];
  data[5] = R[1][1];
  data[6] = R[1][2];
  data[8] = R[2][0];
  data[9] = R[2][1];
  data[10] = R[2][2];
}
 
void vpHomogeneousMatrix::print() const
{
  vpPoseVector r(*this);
  std::cout << r.t();
}
 
void vpHomogeneousMatrix::convert(std::vector<float> &M)
{
  M.resize(12);
  for (unsigned int i = 0; i < 12; i++)
    M[i] = (float)(this->data[i]);
}
 
void vpHomogeneousMatrix::convert(std::vector<double> &M)
{
  M.resize(12);
  for (unsigned int i = 0; i < 12; i++)
    M[i] = this->data[i];
}
 
vpTranslationVector vpHomogeneousMatrix::getTranslationVector() const
{
  vpTranslationVector tr;
  this->extract(tr);
  return tr;
}
 
vpRotationMatrix vpHomogeneousMatrix::getRotationMatrix() const
{
  vpRotationMatrix R;
  this->extract(R);
  return R;
}
 
vpThetaUVector vpHomogeneousMatrix::getThetaUVector() const
{
  vpThetaUVector tu;
  this->extract(tu);
  return tu;
}
 
vpColVector vpHomogeneousMatrix::getCol(unsigned int j) const
{
  if (j >= getCols())
    throw(vpException(vpException::dimensionError, "Unable to extract a column vector from the homogeneous matrix"));
  unsigned int nb_rows = getRows();
  vpColVector c(nb_rows);
  for (unsigned int i = 0; i < nb_rows; i++)
    c[i] = (*this)[i][j];
  return c;
}
 
vpHomogeneousMatrix vpHomogeneousMatrix::compute3d3dTransformation(const std::vector<vpPoint> &p, const std::vector<vpPoint> &q)
{
  const double N = static_cast<double>(p.size());
 
  vpColVector p_bar(3, 0.0);
  vpColVector q_bar(3, 0.0);
  for (size_t i = 0; i < p.size(); i++) {
    for (unsigned int j = 0; j < 3; j++) {
      p_bar[j] += p.at(i).oP[j];
      q_bar[j] += q.at(i).oP[j];
    }
  }
 
  for (unsigned int j = 0; j < 3; j++) {
    p_bar[j] /= N;
    q_bar[j] /= N;
  }
 
  vpMatrix pc(static_cast<unsigned int>(p.size()), 3);
  vpMatrix qc(static_cast<unsigned int>(q.size()), 3);
 
  for (unsigned int i = 0; i < static_cast<unsigned int>(p.size()); i++) {
    for (unsigned int j = 0; j < 3; j++) {
      pc[i][j] = p.at(i).oP[j] - p_bar[j];
      qc[i][j] = q.at(i).oP[j] - q_bar[j];
    }
  }
 
  const vpMatrix pct_qc = pc.t() * qc;
  vpMatrix U = pct_qc, V;
  vpColVector W;
  U.svd(W, V);
 
  vpMatrix Vt = V.t();
  vpMatrix R = U * Vt;
  if (R.det() < 0) {
    Vt[2][0] *= -1;
    Vt[2][1] *= -1;
    Vt[2][2] *= -1;
 
    R = U * Vt;
  }
 
  const vpColVector t = p_bar - R * q_bar;
 
  return vpHomogeneousMatrix(vpTranslationVector(t), vpRotationMatrix(R));
}
 
vpHomogeneousMatrix vpHomogeneousMatrix::mean(const std::vector<vpHomogeneousMatrix> &vec_M)
{
  vpMatrix meanR(3, 3);
  vpColVector meanT(3);
  vpRotationMatrix R;
  for (size_t i = 0; i < vec_M.size(); i++) {
    R = vec_M[i].getRotationMatrix();
    meanR += (vpMatrix)R;
    meanT += (vpColVector)vec_M[i].getTranslationVector();
  }
  meanR /= static_cast<double>(vec_M.size());
  meanT /= static_cast<double>(vec_M.size());
 
  // Euclidean mean of the rotation matrix following Moakher's method (SIAM 2002)
  vpMatrix M, U, V;
  vpColVector sv;
  meanR.pseudoInverse(M, sv, 1e-6, U, V);
  double det = sv[0] * sv[1] * sv[2];
  if (det > 0) {
    meanR = U * V.t();
  } else {
    vpMatrix D(3, 3);
    D = 0.0;
    D[0][0] = D[1][1] = 1.0;
    D[2][2] = -1;
    meanR = U * D * V.t();
  }
 
  R = meanR;
 
  vpTranslationVector t(meanT);
  vpHomogeneousMatrix mean(t, R);
  return mean;
}
 
#if defined(VISP_BUILD_DEPRECATED_FUNCTIONS)
 
void vpHomogeneousMatrix::setIdentity() { eye(); }
 
#endif //#if defined(VISP_BUILD_DEPRECATED_FUNCTIONS)