vpHomogeneousMatrix.cpp

/*
 * ViSP, open source Visual Servoing Platform software.
 * Copyright (C) 2005 - 2025 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 https://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.
 */

 
#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>
 
BEGIN_VISP_NAMESPACE
const unsigned int vpHomogeneousMatrix::constr_value_4 = 4;
vpHomogeneousMatrix::vpHomogeneousMatrix(const vpTranslationVector &t, const vpQuaternionVector &q)
  : vpArray2D<double>(constr_value_4, constr_value_4)
{
  const unsigned int index_3 = 3;
  buildFrom(t, q);
  (*this)[index_3][index_3] = 1.;
}

vpHomogeneousMatrix::vpHomogeneousMatrix() : vpArray2D<double>(constr_value_4, constr_value_4), m_index(0) { eye(); }

vpHomogeneousMatrix::vpHomogeneousMatrix(const vpHomogeneousMatrix &M) : vpArray2D<double>(constr_value_4, constr_value_4), m_index(0)
{
  *this = M;
}

vpHomogeneousMatrix::vpHomogeneousMatrix(const vpTranslationVector &t, const vpThetaUVector &tu)
  : vpArray2D<double>(constr_value_4, constr_value_4), m_index(0)
{
  const unsigned int index_3 = 3;
  buildFrom(t, tu);
  (*this)[index_3][index_3] = 1.;
}

vpHomogeneousMatrix::vpHomogeneousMatrix(const vpTranslationVector &t, const vpRotationMatrix &R)
  : vpArray2D<double>(constr_value_4, constr_value_4), m_index(0)
{
  const unsigned int index_3 = 3;
  insert(R);
  insert(t);
  (*this)[index_3][index_3] = 1.;
}

vpHomogeneousMatrix::vpHomogeneousMatrix(const vpPoseVector &p) : vpArray2D<double>(constr_value_4, constr_value_4), m_index(0)
{
  const unsigned int index_0 = 0;
  const unsigned int index_1 = 1;
  const unsigned int index_2 = 2;
  const unsigned int index_3 = 3;
  const unsigned int index_4 = 4;
  const unsigned int index_5 = 5;
  buildFrom(p[index_0], p[index_1], p[index_2], p[index_3], p[index_4], p[index_5]);
  (*this)[index_3][index_3] = 1.;
}

vpHomogeneousMatrix::vpHomogeneousMatrix(const std::vector<float> &v) : vpArray2D<double>(constr_value_4, constr_value_4), m_index(0)
{
  const unsigned int index_3 = 3;
  buildFrom(v);
  (*this)[index_3][index_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
      orthogonalizeRotation();
    }
    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>(constr_value_4, constr_value_4), m_index(0)
{
  const unsigned int index_3 = 3;
  buildFrom(v);
  (*this)[index_3][index_3] = 1.;
}

vpHomogeneousMatrix::vpHomogeneousMatrix(double tx, double ty, double tz, double tux, double tuy, double tuz)
  : vpArray2D<double>(constr_value_4, constr_value_4), m_index(0)
{
  const unsigned int index_3 = 3;
  buildFrom(tx, ty, tz, tux, tuy, tuz);
  (*this)[index_3][index_3] = 1.;
}

vpHomogeneousMatrix &vpHomogeneousMatrix::buildFrom(const vpTranslationVector &t, const vpThetaUVector &tu)
{
  insert(tu);
  insert(t);
  return *this;
}

vpHomogeneousMatrix &vpHomogeneousMatrix::buildFrom(const vpTranslationVector &t, const vpRotationMatrix &R)
{
  insert(R);
  insert(t);
  return *this;
}

vpHomogeneousMatrix &vpHomogeneousMatrix::buildFrom(const vpPoseVector &p)
{
  const unsigned int index_0 = 0;
  const unsigned int index_1 = 1;
  const unsigned int index_2 = 2;
  const unsigned int index_3 = 3;
  const unsigned int index_4 = 4;
  const unsigned int index_5 = 5;
  vpTranslationVector tv(p[index_0], p[index_1], p[index_2]);
  vpThetaUVector tu(p[index_3], p[index_4], p[index_5]);
 
  insert(tu);
  insert(tv);
  return *this;
}

vpHomogeneousMatrix &vpHomogeneousMatrix::buildFrom(const vpTranslationVector &t, const vpQuaternionVector &q)
{
  insert(t);
  insert(q);
  return *this;
}

vpHomogeneousMatrix &vpHomogeneousMatrix::buildFrom(const double &tx, const double &ty, const double &tz, const double &tux, const double &tuy, const double &tuz)
{
  vpRotationMatrix R(tux, tuy, tuz);
  vpTranslationVector t(tx, ty, tz);
 
  insert(R);
  insert(t);
  return *this;
}

vpHomogeneousMatrix &vpHomogeneousMatrix::buildFrom(const std::vector<float> &v)
{
  const std::size_t val_12 = 12;
  const std::size_t val_16 = 16;
  const unsigned int val_12ui = 12, val_16ui = 16;
  if ((v.size() != val_12) && (v.size() != val_16)) {
    throw(vpException(vpException::dimensionError, "Cannot convert std::vector<float> to vpHomogeneousMatrix"));
  }
 
  for (unsigned int i = 0; i < val_12ui; ++i) {
    this->data[i] = static_cast<double>(v[i]);
  }
 
  if (v.size() == val_12) {
    data[12] = 0.;
    data[13] = 0.;
    data[14] = 0.;
    data[15] = 1.;
  }
  else {
    for (unsigned int i = val_12ui; i < val_16ui; ++i) {
      this->data[i] = static_cast<double>(v[i]);
    }
  }
 
  if (!isAnHomogeneousMatrix()) {
    if (isAnHomogeneousMatrix(1e-3)) {
      // re-orthogonalize rotation matrix since the input is close to a valid rotation matrix
      orthogonalizeRotation();
    }
    else {
      throw(vpException(
        vpException::fatalError,
        "Homogeneous matrix initialization failed since its elements are not valid (rotation part or last row)"));
    }
  }
  return *this;
}

vpHomogeneousMatrix &vpHomogeneousMatrix::buildFrom(const std::vector<double> &v)
{
  const std::size_t val_12 = 12;
  const std::size_t val_16 = 16;
  const unsigned int val_12ui = 12, val_16ui = 16;
 
 
  if ((v.size() != val_12) && (v.size() != val_16)) {
    throw(vpException(vpException::dimensionError, "Cannot convert std::vector<double> to vpHomogeneousMatrix"));
  }
 
  for (unsigned int i = 0; i < val_12ui; ++i) {
    this->data[i] = static_cast<double>(v[i]);
  }
 
  if (v.size() == val_12) {
    data[12] = 0.;
    data[13] = 0.;
    data[14] = 0.;
    data[15] = 1.;
  }
  else {
    for (unsigned int i = val_12ui; i < val_16ui; ++i) {
      this->data[i] = static_cast<double>(v[i]);
    }
  }
 
  if (!isAnHomogeneousMatrix()) {
    if (isAnHomogeneousMatrix(1e-3)) {
      // re-orthogonalize rotation matrix since the input is close to a valid rotation matrix
      orthogonalizeRotation();
    }
    else {
      throw(vpException(
        vpException::fatalError,
        "Homogeneous matrix initialization fails since its elements are not valid (rotation part or last row)"));
    }
  }
 
  return *this;
}

vpHomogeneousMatrix &vpHomogeneousMatrix::operator=(const vpHomogeneousMatrix &M)
{
  const int val_4 = 4;
  for (int i = 0; i < val_4; ++i) {
    for (int j = 0; j < val_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
{
  const unsigned int val_4 = 4;
  if (v.getRows() != val_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 < val_4; ++j) {
    for (unsigned int i = 0; i < val_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);
 
  const unsigned int index_0 = 0;
  const unsigned int index_1 = 1;
  const unsigned int index_2 = 2;
  const unsigned int index_3 = 3;
 
  v[index_0] = bP.get_X();
  v[index_1] = bP.get_Y();
  v[index_2] = bP.get_Z();
  v[index_3] = bP.get_W();
 
  v1[index_0] = ((*this)[index_0][index_0] * v[index_0]) + ((*this)[index_0][index_1] * v[index_1]) + ((*this)[index_0][index_2] * v[index_2]) + ((*this)[index_0][index_3] * v[index_3]);
  v1[index_1] = ((*this)[index_1][index_0] * v[index_0]) + ((*this)[index_1][index_1] * v[index_1]) + ((*this)[index_1][index_2] * v[index_2]) + ((*this)[index_1][index_3] * v[index_3]);
  v1[index_2] = ((*this)[index_2][index_0] * v[index_0]) + ((*this)[index_2][index_1] * v[index_1]) + ((*this)[index_2][index_2] * v[index_2]) + ((*this)[index_2][index_3] * v[index_3]);
  v1[index_3] = ((*this)[index_3][index_0] * v[index_0]) + ((*this)[index_3][index_1] * v[index_1]) + ((*this)[index_3][index_2] * v[index_2]) + ((*this)[index_3][index_3] * v[index_3]);
 
  v1 /= v1[index_3];
 
  //  --comment: v1 equals M multiplied by v
  aP.set_X(v1[index_0]);
  aP.set_Y(v1[index_1]);
  aP.set_Z(v1[index_2]);
  aP.set_W(v1[index_3]);
 
  aP.set_oX(v1[index_0]);
  aP.set_oY(v1[index_1]);
  aP.set_oZ(v1[index_2]);
  aP.set_oW(v1[index_3]);
 
  return aP;
}

vpTranslationVector vpHomogeneousMatrix::operator*(const vpTranslationVector &t) const
{
  vpTranslationVector t_out;
  const unsigned int index_0 = 0;
  const unsigned int index_1 = 1;
  const unsigned int index_2 = 2;
  const unsigned int index_3 = 3;
  t_out[index_0] = (((*this)[index_0][0] * t[0]) + ((*this)[index_0][1] * t[1]) + ((*this)[index_0][index_2] * t[index_2])) + (*this)[index_0][index_3];
  t_out[index_1] = (((*this)[index_1][0] * t[0]) + ((*this)[index_1][1] * t[1]) + ((*this)[index_1][index_2] * t[index_2])) + (*this)[index_1][index_3];
  t_out[index_2] = (((*this)[index_2][0] * t[0]) + ((*this)[index_2][1] * t[1]) + ((*this)[index_2][index_2] * t[index_2])) + (*this)[index_2][index_3];
 
  return t_out;
}

vpHomogeneousMatrix vpHomogeneousMatrix::operator*(const vpRotationMatrix &R) const
{
  return (vpHomogeneousMatrix((*this).getTranslationVector(), (*this).getRotationMatrix() * R));
}

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 unsigned int index_0 = 0;
  const unsigned int index_1 = 1;
  const unsigned int index_2 = 2;
  const unsigned int index_3 = 3;
  const double epsilon = std::numeric_limits<double>::epsilon();
  bool isLastRowOK = vpMath::nul((*this)[index_3][index_0], epsilon) && vpMath::nul((*this)[index_3][index_1], epsilon) &&
    vpMath::nul((*this)[index_3][index_2], epsilon);
  return R.isARotationMatrix(threshold) && isLastRowOK && vpMath::equal((*this)[index_3][index_3], 1.0, epsilon);
}

bool vpHomogeneousMatrix::isValid() const
{
  unsigned int l_size = size();
  for (unsigned int i = 0; i < l_size; ++i) {
    if (vpMath::isNaN(data[i])) {
      return false;
    }
  }
  return true;
}

void vpHomogeneousMatrix::extract(vpRotationMatrix &R) const
{
  const unsigned int val_3 = 3;
  for (unsigned int i = 0; i < val_3; ++i) {
    for (unsigned int j = 0; j < val_3; ++j) {
      R[i][j] = (*this)[i][j];
    }
  }
}

void vpHomogeneousMatrix::extract(vpTranslationVector &t) const
{
  const unsigned int index_0 = 0;
  const unsigned int index_1 = 1;
  const unsigned int index_2 = 2;
  const unsigned int index_3 = 3;
  t[index_0] = (*this)[index_0][index_3];
  t[index_1] = (*this)[index_1][index_3];
  t[index_2] = (*this)[index_2][index_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)
{
  const unsigned int val_3 = 3;
  for (unsigned int i = 0; i < val_3; ++i) {
    for (unsigned int j = 0; j < val_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)
{
  const unsigned int index_0 = 0;
  const unsigned int index_1 = 1;
  const unsigned int index_2 = 2;
  const unsigned int index_3 = 3;
  (*this)[index_0][index_3] = t[index_0];
  (*this)[index_1][index_3] = t[index_1];
  (*this)[index_2][index_3] = t[index_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()
{
  const unsigned int index_0 = 0;
  const unsigned int index_1 = 1;
  const unsigned int index_2 = 2;
  const unsigned int index_3 = 3;
  (*this)[index_0][index_0] = 1;
  (*this)[index_1][index_1] = 1;
  (*this)[index_2][index_2] = 1;
  (*this)[index_3][index_3] = 1;
 
  (*this)[index_0][index_1] = 0;
  (*this)[index_0][index_2] = 0;
  (*this)[index_0][index_3] = 0;
  (*this)[index_1][index_0] = 0;
  (*this)[index_1][index_2] = 0;
  (*this)[index_1][index_3] = 0;
  (*this)[index_2][index_0] = 0;
  (*this)[index_2][index_1] = 0;
  (*this)[index_2][index_3] = 0;
  (*this)[index_3][index_0] = 0;
  (*this)[index_3][index_1] = 0;
  (*this)[index_3][index_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::save(const std::string &filename) const
{
  std::ofstream f;
  f.open(filename.c_str());
  save(f);
  f.close();
}
 
void vpHomogeneousMatrix::load(std::ifstream &f)
{
  if (!f.fail()) {
    const unsigned int val_4 = 4;
    for (unsigned int i = 0; i < val_4; ++i) {
      for (unsigned int j = 0; j < val_4; ++j) {
        f >> (*this)[i][j];
      }
    }
  }
  else {
    throw(vpException(vpException::ioError, "Cannot load homogeneous matrix: ifstream not open"));
  }
}
 
void vpHomogeneousMatrix::load(const std::string &filename)
{
  std::ifstream f;
  f.open(filename.c_str());
  load(f);
  f.close();
}

void vpHomogeneousMatrix::orthogonalizeRotation()
{
  vpRotationMatrix R(*this);
  R.orthogonalize();
  insert(R);
}

void vpHomogeneousMatrix::print() const
{
  vpPoseVector r(*this);
  std::cout << r.t();
}

void vpHomogeneousMatrix::convert(std::vector<float> &M)
{
  const std::size_t val_12 = 12;
  const unsigned int val_12ui = 12;
  M.resize(val_12);
  for (unsigned int i = 0; i < val_12ui; ++i) {
    M[i] = static_cast<float>(this->data[i]);
  }
}

void vpHomogeneousMatrix::convert(std::vector<double> &M)
{
  const std::size_t val_12 = 12;
  const unsigned int val_12ui = 12;
  M.resize(val_12);
  for (unsigned int i = 0; i < val_12ui; ++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);
  size_t p_size = p.size();
  const unsigned int val_3 = 3;
  for (size_t i = 0; i < p_size; ++i) {
    for (unsigned int j = 0; j < val_3; ++j) {
      p_bar[j] += p.at(i).oP[j];
      q_bar[j] += q.at(i).oP[j];
    }
  }
 
  for (unsigned int j = 0; j < val_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 < val_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) {
    const unsigned int index_0 = 0;
    const unsigned int index_1 = 1;
    const unsigned int index_2 = 2;
    Vt[index_2][index_0] *= -1.;
    Vt[index_2][index_1] *= -1.;
    Vt[index_2][index_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;
  const unsigned int index_0 = 0;
  const unsigned int index_1 = 1;
  const unsigned int index_2 = 2;
  size_t vec_m_size = vec_M.size();
  for (size_t i = 0; i < vec_m_size; ++i) {
    R = vec_M[i].getRotationMatrix();
    meanR += static_cast<vpMatrix>(R);
    meanT += static_cast<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[index_0] * sv[index_1] * sv[index_2];
  if (det > 0) {
    meanR = U * V.t();
  }
  else {
    vpMatrix D(3, 3);
    D = 0.0;
    D[index_0][index_0] = 1.0;
    D[index_1][index_1] = 1.0;
    D[index_2][index_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)
 
#ifdef VISP_HAVE_NLOHMANN_JSON
VP_ATTRIBUTE_NO_DESTROY const std::string vpHomogeneousMatrix::jsonTypeName = "vpHomogeneousMatrix";
#include <visp3/core/vpJsonParsing.h>
void vpHomogeneousMatrix::convert_to_json(nlohmann::json &j) const
{
  const vpArray2D<double> *asArray = (vpArray2D<double>*) this;
  to_json(j, *asArray);
  j["type"] = vpHomogeneousMatrix::jsonTypeName;
}
 
void vpHomogeneousMatrix::parse_json(const nlohmann::json &j)
{
#ifdef ENABLE_VISP_NAMESPACE
  using namespace VISP_NAMESPACE_NAME;
#endif
  vpArray2D<double> *asArray = (vpArray2D<double>*) this;
  if (j.is_object() && j.contains("type")) { // Specific conversions
    const bool converted = convertFromTypeAndBuildFrom<vpHomogeneousMatrix, vpPoseVector>(j, *this);
    if (!converted) {
      from_json(j, *asArray);
    }
  }
  else { // Generic 2D array conversion
    from_json(j, *asArray);
  }
 
  if ((getCols() != 4) && (getRows() != 4)) {
    throw vpException(vpException::badValue, "From JSON, tried to read something that is not a 4x4 matrix");
  }
  if (!isAnHomogeneousMatrix()) {
    throw vpException(vpException::badValue, "From JSON read a non homogeneous matrix into a vpHomogeneousMatrix");
  }
}
#endif
END_VISP_NAMESPACE