vnl_matrix_fixed.txx

InsightToolkit-1.4.0(有大量的优化算法程序)

// This is vxl/vnl/vnl_matrix_fixed.txx
#ifndef vnl_matrix_fixed_txx_
#define vnl_matrix_fixed_txx_
//:
// \file
#include "vnl_matrix_fixed.h"

#include <vcl_cmath.h>
#include <vcl_iostream.h>
#include <vcl_cstdlib.h> // for abort
#include <vcl_cassert.h>

#include "vnl_error.h"
#include "vnl_math.h"

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::add( const T* a, const T* b, T* r )
{
  unsigned int count = nrows*ncols;
  while ( count-- )
    *(r++) = *(a++) + *(b++);
}


template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::add( const T* a, T b, T* r )
{
  unsigned int count = nrows*ncols;
  while ( count-- )
    *(r++) = *(a++) + b;
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::sub( const T* a, const T* b, T* r )
{
  unsigned int count = nrows*ncols;
  while ( count-- )
    *(r++) = *(a++) - *(b++);
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::sub( const T* a, T b, T* r )
{
  unsigned int count = nrows*ncols;
  while ( count-- )
    *(r++) = *(a++) - b;
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::sub( T a, const T* b, T* r )
{
  unsigned int count = nrows*ncols;
  while ( count-- )
    *(r++) = a - *(b++);
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::mul( const T* a, const T* b, T* r )
{
  unsigned int count = nrows*ncols;
  while ( count-- )
    *(r++) = *(a++) * *(b++);
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::mul( const T* a, T b, T* r )
{
  unsigned int count = nrows*ncols;
  while ( count-- )
    *(r++) = *(a++) * b;
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::div( const T* a, const T* b, T* r )
{
  unsigned int count = nrows*ncols;
  while ( count-- )
    *(r++) = *(a++) / *(b++);
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::div( const T* a, T b, T* r )
{
  unsigned int count = nrows*ncols;
  while ( count-- )
    *(r++) = *(a++) / b;
}

template<class T, unsigned nrows, unsigned ncols>
bool
vnl_matrix_fixed<T,nrows,ncols>::equal( const T* a, const T* b )
{
  unsigned int count = nrows*ncols;
  while ( count-- )
    if ( *(a++) != *(b++) )  return false;
  return true;
}

//------------------------------------------------------------


template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::fill (T value)
{
  for (unsigned int i = 0; i < nrows; i++)
    for (unsigned int j = 0; j < ncols; j++)
      (*this)(i,j) = value;
}


template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::fill_diagonal (T value)
{
  for (unsigned int i = 0; i < nrows && i < ncols; i++)
    (*this)(i,i) = value;
}


template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::print(vcl_ostream& os) const
{
  for (unsigned int i = 0; i < nrows; i++)
  {
    os << (*this)(i,0);
    for (unsigned int j = 1; j < ncols; j++)
      os << " " << (*this)(i,j);
    os << "\n";
  }
}


template <class T, unsigned nrows, unsigned ncols>
vnl_matrix_fixed<T,nrows,ncols>
vnl_matrix_fixed<T,nrows,ncols>::apply(T (*f)(T const&)) const
{
  vnl_matrix_fixed<T,nrows,ncols> ret;
  vnl_c_vector<T>::apply(this->data_[0], rows()*cols(), f, ret.data_block());
  return ret;
}

template <class T, unsigned nrows, unsigned ncols>
vnl_matrix_fixed<T,nrows,ncols>
vnl_matrix_fixed<T,nrows,ncols>::apply(T (*f)(T)) const
{
  vnl_matrix_fixed<T,nrows,ncols> ret;
  vnl_c_vector<T>::apply(this->data_[0], rows()*cols(), f, ret.data_block());
  return ret;
}

////--------------------------- Additions------------------------------------


template<class T, unsigned nrows, unsigned ncols>
vnl_matrix_fixed<T,ncols,nrows>
vnl_matrix_fixed<T,nrows,ncols>::transpose() const
{
  vnl_matrix_fixed<T,ncols,nrows> result;
  for (unsigned int i = 0; i < cols(); i++)
    for (unsigned int j = 0; j < rows(); j++)
      result(i,j) = (*this)(j,i);
  return result;
}

template<class T, unsigned nrows, unsigned ncols>
vnl_matrix_fixed<T,ncols,nrows>
vnl_matrix_fixed<T,nrows,ncols>::conjugate_transpose() const
{
  vnl_matrix_fixed<T,ncols,nrows> result(transpose());
  vnl_c_vector<T>::conjugate(result.begin(),  // src
                             result.begin(),  // dst
                             result.size());  // size of block
  return result;
}

template<class T, unsigned nrows, unsigned ncols>
vnl_matrix_fixed<T,nrows,ncols>&
vnl_matrix_fixed<T,nrows,ncols>::update (vnl_matrix<T> const& m,
                                         unsigned top, unsigned left)
{
  const unsigned int bottom = top + m.rows();
  const unsigned int right = left + m.cols();
#ifndef NDEBUG
  if (nrows < bottom || ncols < right)
    vnl_error_matrix_dimension ("update",
                                bottom, right, m.rows(), m.cols());
#endif
  for (unsigned int i = top; i < bottom; i++)
    for (unsigned int j = left; j < right; j++)
      (*this)(i,j) = m(i-top,j-left);
  return *this;
}


template<class T, unsigned nrows, unsigned ncols>
vnl_matrix<T>
vnl_matrix_fixed<T,nrows,ncols>::extract (unsigned rowz, unsigned colz,
                                          unsigned top, unsigned left) const
{
#ifndef NDEBUG
  unsigned int bottom = top + rowz;
  unsigned int right = left + colz;
  if ((nrows < bottom) || (ncols < right))
    vnl_error_matrix_dimension ("extract",
                                nrows, ncols, bottom, right);
#endif
  vnl_matrix<T> result(rowz, colz);
  for (unsigned int i = 0; i < rowz; i++)      // actual copy of all elements
    for (unsigned int j = 0; j < colz; j++)    // in submatrix
      result(i,j) = (*this)(top+i,left+j);
  return result;
}


template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::copy_in(T const *p)
{
  T* dp = this->data_block();
  unsigned int i = nrows * ncols;
  while (i--)
    *dp++ = *p++;
}

template<class T, unsigned nrows, unsigned ncols>
void vnl_matrix_fixed<T,nrows,ncols>::copy_out(T *p) const
{
  T const* dp = this->data_block();
  unsigned int i = nrows*ncols;
  while (i--)
    *p++ = *dp++;
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::set_identity()
{
#ifndef NDEBUG
  if (nrows != ncols) // Size?
    vnl_error_matrix_nonsquare ("set_identity");
#endif
  // Two simple loops are generally better than having a branch inside
  // the loop. Probably worth the O(n) extra writes.
  for (unsigned int i = 0; i < nrows; i++)
    for (unsigned int j = 0; j < ncols; j++)
        (*this)(i,j) = T(0);
  for (unsigned int i = 0; i < nrows; i++)
    (*this)(i,i) = T(1);
}

//: Make each row of the matrix have unit norm.
// All-zero rows are ignored.
template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::normalize_rows()
{
  typedef typename vnl_numeric_traits<T>::abs_t abs_t;
  for (unsigned int i = 0; i < nrows; i++)
  {
    abs_t norm(0); // double will not do for all types.
    for (unsigned int j = 0; j < ncols; j++)
      norm += vnl_math_squared_magnitude( (*this)(i,j) );

    if (norm != 0)
    {
      typedef typename vnl_numeric_traits<abs_t>::real_t real_t;
      real_t scale = real_t(1)/vcl_sqrt((real_t)norm);
      for (unsigned int j = 0; j < ncols; j++)
      {
        // FIXME need correct rounding here
        // There is e.g. no *standard* operator*=(complex<float>, double), hence the T() cast.
        (*this)(i,j) *= (T)(scale);
      }
    }
  }
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::normalize_columns()
{
  typedef typename vnl_numeric_traits<T>::abs_t abs_t;
  for (unsigned int j = 0; j < ncols; j++) {  // For each column in the Matrix
    abs_t norm(0); // double will not do for all types.
    for (unsigned int i = 0; i < nrows; i++)
      norm += vnl_math_squared_magnitude( (*this)(i,j) );

    if (norm != 0)
    {
      typedef typename vnl_numeric_traits<abs_t>::real_t real_t;
      real_t scale = real_t(1)/vcl_sqrt((real_t)norm);
      for (unsigned int i = 0; i < nrows; i++)
      {
        // FIXME need correct rounding here
        // There is e.g. no *standard* operator*=(complex<float>, double), hence the T() cast.
        (*this)(i,j) *= (T)(scale);
      }
    }
  }
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::scale_row(unsigned row_index, T value)
{
#ifndef NDEBUG
  if (row_index >= nrows)
    vnl_error_matrix_row_index("scale_row", row_index);
#endif
  for (unsigned int j = 0; j < ncols; j++)
    (*this)(row_index,j) *= value;
}

template<class T, unsigned nrows, unsigned ncols>
void
vnl_matrix_fixed<T,nrows,ncols>::scale_column(unsigned column_index, T value)
{
#ifndef NDEBUG
  if (column_index >= ncols)
    vnl_error_matrix_col_index("scale_column", column_index);
#endif
  for (unsigned int j = 0; j < nrows; j++)
    (*this)(j,column_index) *= value;
}

//: Returns a copy of n rows, starting from "row"
template<class T, unsigned nrows, unsigned ncols>
vnl_matrix<T>
vnl_matrix_fixed<T,nrows,ncols>::get_n_rows (unsigned row, unsigned n) const
{
#ifndef NDEBUG
  if (row + n > nrows)
    vnl_error_matrix_row_index ("get_n_rows", row);
#endif

  // Extract data rowwise.
  return vnl_matrix<T>(data_[row], n, ncols);
}

template<class T, unsigned nrows, unsigned ncols>