itkMathDeterminant.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkMathDeterminant_h
#define itkMathDeterminant_h
 
#include "itkMacro.h"
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_matrix_fixed.h"
 
#include "itk_eigen.h"
#include ITK_EIGEN(Dense)
 
namespace itk
{
// Forward declaration lets the Matrix overload parse under the circular include
// with itkMatrix.h; its body instantiates only where itk::Matrix is complete.
template <typename T, unsigned int VRows, unsigned int VColumns>
class Matrix;
 
namespace Math
{
namespace detail
{
// Eigen chooses direct cofactor formulas for VDim <= 4 and PartialPivLU beyond;
// the determinant is transpose-invariant, so the row-major map matches ITK
// storage without affecting the value.
template <unsigned int VDim, typename TReal>
TReal
DeterminantEigen(const TReal * inData)
{
  using RowMajor = Eigen::Matrix<TReal, VDim, VDim, Eigen::RowMajor>;
  return Eigen::Map<const RowMajor>(inData).determinant();
}
 
template <typename TReal>
TReal
DynamicDeterminantEigen(const TReal * inData, unsigned int n)
{
  // Eigen's Dynamic-sized determinant always uses LU; dispatch small sizes to
  // the compile-time direct formulas so runtime small matrices stay fast.
  switch (n)
  {
    case 1:
      return inData[0];
    case 2:
      return DeterminantEigen<2, TReal>(inData);
    case 3:
      return DeterminantEigen<3, TReal>(inData);
    case 4:
      return DeterminantEigen<4, TReal>(inData);
    default:
    {
      using RowMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
      return Eigen::Map<const RowMajor>(inData, n, n).determinant();
    }
  }
}
} // namespace detail

template <typename TReal, unsigned int VDim>
TReal
Determinant(const vnl_matrix_fixed<TReal, VDim, VDim> & A)
{
  return detail::DeterminantEigen<VDim, TReal>(A.data_block());
}

template <typename TReal, unsigned int VDim>
TReal
Determinant(const Matrix<TReal, VDim, VDim> & A)
{
  return Determinant<TReal, VDim>(A.GetVnlMatrix());
}

template <typename TReal>
TReal
Determinant(const vnl_matrix<TReal> & A)
{
  const unsigned int rows = A.rows();
  if (rows != A.cols())
  {
    itkGenericExceptionMacro("itk::Math::Determinant requires a square matrix.");
  }
  return detail::DynamicDeterminantEigen<TReal>(A.data_block(), rows);
}

template <typename TReal, unsigned int VRows, unsigned int VColumns>
TReal
Determinant(const vnl_matrix_fixed<TReal, VRows, VColumns> & A)
{
  return Determinant(A.as_ref());
}
 
} // namespace Math
} // namespace itk
 
#endif // itkMathDeterminant_h