itkEigenDecompositionSignConvention.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 itkEigenDecompositionSignConvention_h
#define itkEigenDecompositionSignConvention_h
 
#include "vnl/vnl_matrix.h"
#include <cmath>
 
namespace itk::detail
{

template <typename T>
void
CanonicalizeEigenvectorColumnSigns(vnl_matrix<T> & V)
{
  const unsigned int rows = V.rows();
  for (unsigned int j = 0; j < V.cols(); ++j)
  {
    unsigned int pivot = 0;
    for (unsigned int i = 1; i < rows; ++i)
    {
      if (std::abs(V(i, j)) > std::abs(V(pivot, j)))
      {
        pivot = i;
      }
    }
    if (V(pivot, j) < T{ 0 })
    {
      V.scale_column(j, T{ -1 });
    }
  }
}

template <typename TMatrixU, typename TMatrixV>
void
CanonicalizeColumnSignsPaired(TMatrixU & u, TMatrixV & paired)
{
  // u and paired share a column count but may differ in row count (an SVD's
  // thin U is rows x k, V is cols x k), so each is flipped over its own rows.
  const unsigned int uRows = u.rows();
  const unsigned int pairedRows = paired.rows();
  for (unsigned int j = 0; j < u.cols(); ++j)
  {
    unsigned int pivot = 0;
    for (unsigned int i = 1; i < uRows; ++i)
    {
      if (std::abs(u(i, j)) > std::abs(u(pivot, j)))
      {
        pivot = i;
      }
    }
    if (u(pivot, j) < typename TMatrixU::element_type{ 0 })
    {
      for (unsigned int i = 0; i < uRows; ++i)
      {
        u(i, j) = -u(i, j);
      }
      for (unsigned int i = 0; i < pairedRows; ++i)
      {
        paired(i, j) = -paired(i, j);
      }
    }
  }
}
 
} // namespace itk::detail
 
#endif // itkEigenDecompositionSignConvention_h