SimpleMatrixSpecializations.ih

#include <cstdlib>
// IMPLEMENTATION of inline methods.

// ----------------------- Standard services ------------------------------
template <typename M, DGtal::Dimension TM, DGtal::Dimension TN>
inline 
typename DGtal::SimpleMatrixSpecializations<M,TM,TN>::Component 
DGtal::SimpleMatrixSpecializations<M,TM,TN>::minorDeterminant(const Matrix &aM,
                                                       const DGtal::Dimension ai, 
                                                       const DGtal::Dimension aj)
{
  BOOST_STATIC_ASSERT(TM == TN);
  ASSERT(ai<TM);
  ASSERT(aj<TN);
  DGtal::SimpleMatrix<Component,TM-1,TN-1> mat;
  DGtal::Dimension indexR=0;
  DGtal::Dimension indexC=0;
  for (DGtal::Dimension i=0; i<TM; i++)
    for (DGtal::Dimension j=0; j<TN; j++)
      {
        if (i!=ai && j!=aj)
          {
            ASSERT(indexR < TM -1);
            ASSERT(indexC < TN -1);
            mat.setComponent(indexR,indexC, aM(i,j));
            indexC++;
          }
        if (indexC==TM-1)
          {
            indexC=0;
            indexR++;
          }

        if (indexR==TM-1)
          return mat.determinant();
      }
 
  return mat.determinant();
}

template <typename M, DGtal::Dimension TM, DGtal::Dimension TN>
inline 
typename DGtal::SimpleMatrixSpecializations<M,TM,TN>::Component 
DGtal::SimpleMatrixSpecializations<M,TM,TN>::determinant(const Matrix &aM)
{
  BOOST_STATIC_ASSERT(TM == TN);
  
  Component d = DGtal::NumberTraits<Component>::ZERO;
  for(DGtal::Dimension i= 0; i< TM; ++i)
    d += aM(i,0) * aM.cofactor(i,0);
  return d;
}
// ----------------------- Specialization 1x1 ------------------------------
    
template <typename M>
inline 
typename DGtal::SimpleMatrixSpecializations<M,1,1>::Component 
DGtal::SimpleMatrixSpecializations<M,1,1>::minorDeterminant(const Matrix &aM,
                                                            const DGtal::Dimension i, 
                                                            const DGtal::Dimension j)
{
  ASSERT(false && "Not defined for 1x1 matrices");
  return NumberTraits<Component>::ZERO;
}

template <typename M>
inline 
typename DGtal::SimpleMatrixSpecializations<M,1,1>::Component 
DGtal::SimpleMatrixSpecializations<M,1,1>::determinant(const Matrix &aM)
{
 
  return aM(0,0);
}
// ----------------------- Specialization 2x2 ------------------------------
    
template <typename M>
inline 
typename DGtal::SimpleMatrixSpecializations<M,2,2>::Component 
DGtal::SimpleMatrixSpecializations<M,2,2>::minorDeterminant(const Matrix &aM,
                                                       const DGtal::Dimension i, 
                                                       const DGtal::Dimension j)
{
  ASSERT(i<2);
  ASSERT(j<2);
 return aM((i+1) % 2,(j+1) % 2);
}

template <typename M>
inline 
typename DGtal::SimpleMatrixSpecializations<M,2,2>::Component 
DGtal::SimpleMatrixSpecializations<M,2,2>::determinant(const Matrix &aM)
{
 
  return aM(0,0)*aM(1,1) - aM(0,1)*aM(1,0);
}
// ----------------------- Specialization 3x3 ------------------------------

template <typename M>
inline 
typename DGtal::SimpleMatrixSpecializations<M,3,3>::Component 
DGtal::SimpleMatrixSpecializations<M,3,3>::minorDeterminant(const Matrix &aM,
                                                       const DGtal::Dimension ai, 
                                                       const DGtal::Dimension aj)
{
  ASSERT(ai<3);
  ASSERT(aj<3);
  DGtal::SimpleMatrix<Component,2,2> mat;
  DGtal::Dimension indexR=0;
  DGtal::Dimension indexC=0;
  for (DGtal::Dimension i=0; i<3; i++)
    for (DGtal::Dimension j=0; j<3; j++)
      {
        if (i!=ai && j!=aj)
          {
            ASSERT(indexR < 3 -1);
            ASSERT(indexC < 3 -1);
            mat.setComponent(indexR,indexC, aM(i,j));
            indexC++;
          }
        if (indexC==3-1)
          {
            indexC=0;
            indexR++;
          }

        if (indexR==3-1)
          return determinant(mat);
      }
 
  return determinant(mat);
}

template <typename M>
inline 
typename DGtal::SimpleMatrixSpecializations<M,3,3>::Component 
DGtal::SimpleMatrixSpecializations<M,3,3>::determinant(const Matrix &aM)
{
  return aM(0,0) * ( (aM(1,1)*aM(2,2))-
                     (aM(1,2)*aM(2,1)) )-
    aM(1,0) * ( (aM(0,1)*aM(2,2))-
                (aM(0,2)*aM(2,1)) )
    + aM(2,0) * ( (aM(0,1)*aM(1,2))-
                  (aM(0,2)*aM(1,1)) );
}

//                                                                           //