vnl_matrix.txx

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

// This is vxl/vnl/vnl_matrix.txx
#ifndef vnl_matrix_txx_
#define vnl_matrix_txx_
//:
// \file
//
// Copyright (C) 1991 Texas Instruments Incorporated.
// Copyright (C) 1992 General Electric Company.
//
// Permission is granted to any individual or institution to use, copy, modify,
// and distribute this software, provided that this complete copyright and
// permission notice is maintained, intact, in all copies and supporting
// documentation.
//
// Texas Instruments Incorporated, General Electric Company,
// provides this software "as is" without express or implied warranty.
//
// Created: MBN 04/21/89  Initial design and implementation
// Updated: MBN 06/22/89  Removed non-destructive methods
// Updated: LGO 08/09/89  Inherit from Generic
// Updated: MBN 08/20/89  Changed template usage to reflect new syntax
// Updated: MBN 09/11/89  Added conditional exception handling and base class
// Updated: LGO 10/05/89  Don't re-allocate data in operator= when same size
// Updated: LGO 10/19/89  Add extra parameter to varargs constructor
// Updated: MBN 10/19/89  Added optional argument to set_compare method
// Updated: LGO 12/08/89  Allocate column data in one chunk
// Updated: LGO 12/08/89  Clean-up get and put, add const everywhere.
// Updated: LGO 12/19/89  Remove the map and reduce methods
// Updated: MBN 02/22/90  Changed size arguments from int to unsigned int
// Updated: MJF 06/30/90  Added base class name to constructor initializer
// Updated: VDN 02/21/92  New lite version
// Updated: VDN 05/05/92  Use envelope to avoid unecessary copying
// Updated: VDN 09/30/92  Matrix inversion with singular value decomposition
// Updated: AWF 08/21/96  set_identity, normalize_rows, scale_row.
// Updated: AWF 09/30/96  set_row/set_column methods.  Const-correct data_block().
// Updated: AWF 14/02/97  get_n_rows, get_n_columns.
// Updated: PVR 20/03/97  get_row, get_column.
//
// The parameterized vnl_matrix<T> class implements two dimensional arithmetic
// matrices of a user specified type. The only constraint placed on the type is
// that it must overload the following operators: +, -,  *,  and /. Thus, it
// will be possible to have a vnl_matrix over vcl_complex<T>. The vnl_matrix<T>
// class is static in size, that is once a vnl_matrix<T> of a particular size
// has been created, there is no dynamic growth or resize method available.
//
// Each matrix contains  a protected  data section  that has a T** slot that
// points to the  physical memory allocated  for the two  dimensional array. In
// addition, two integers  specify   the number  of  rows  and columns  for the
// matrix.  These values  are provided in the  constructors. A single protected
// slot  contains a pointer  to a compare  function  to   be used  in  equality
// operations. The default function used is the built-in == operator.
//
// Four  different constructors are provided.  The  first constructor takes two
// integer arguments  specifying the  row  and column  size.   Enough memory is
// allocated to hold row*column elements  of type Type.  The second constructor
// takes the  same two  first arguments, but  also accepts  an additional third
// argument that is  a reference to  an  object of  the appropriate  type whose
// value is used as an initial fill value.  The third constructor is similar to
// the third, except that it accpets a variable number of initialization values
// for the Matrix.  If there are  fewer values than elements,  the rest are set
// to zero. Finally, the last constructor takes a single argument consisting of
// a reference to a Matrix and duplicates its size and element values.
//
// Methods   are  provided   for destructive   scalar   and Matrix    addition,
// multiplication, check for equality  and inequality, fill, reduce, and access
// and set individual elements.  Finally, both  the  input and output operators
// are overloaded to allow for fomatted input and output of matrix elements.
//
// Good matrix inversion is needed. We choose singular value decomposition,
// since it is general and works great for nearly singular cases. Singular
// value decomposition is preferred to LU decomposition, since the accuracy
// of the pivots is independent from the left->right top->down elimination.
// LU decomposition also does not give eigenvectors and eigenvalues when
// the matrix is symmetric.
//
// Several different constructors are provided. See .h file for brief descriptions.

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

#include "vnl_matrix.h"

#include <vcl_cassert.h>
#include <vcl_cstdlib.h>  // abort()
#include <vcl_cctype.h>   // isspace()
#include <vcl_iostream.h>
#include <vcl_vector.h>
#include <vcl_algorithm.h>

#include <vnl/vnl_math.h>
#include <vnl/vnl_vector.h>
#include <vnl/vnl_c_vector.h>
#include <vnl/vnl_numeric_traits.h>
//--------------------------------------------------------------------------------

// This macro allocates and initializes the dynamic storage used by a vnl_matrix.
#define vnl_matrix_alloc_blah(rowz_, colz_) \
do { \
  this->num_rows = (rowz_); \
  this->num_cols = (colz_); \
  if (this->num_rows && this->num_cols) { \
    /* Allocate memory to hold the row pointers */ \
    this->data = vnl_c_vector<T>::allocate_Tptr(this->num_rows); \
    /* Allocate memory to hold the elements of the matrix */ \
    T* elmns = vnl_c_vector<T>::allocate_T(this->num_rows * this->num_cols); \
    /* Fill in the array of row pointers */ \
    for (unsigned int i = 0; i < this->num_rows; ++ i) \
      this->data[i] = elmns + i*this->num_cols; \
  } \
  else { \
   /* This is to make sure .begin() and .end() work for 0xN matrices: */ \
   (this->data = vnl_c_vector<T>::allocate_Tptr(1))[0] = 0; \
  } \
} while (false)

// This macro releases the dynamic storage used by a vnl_matrix.
#define vnl_matrix_free_blah \
do { \
  if (this->data) { \
    if (this->num_cols && this->num_rows) { \
      vnl_c_vector<T>::deallocate(this->data[0], this->num_cols * this->num_rows); \
      vnl_c_vector<T>::deallocate(this->data, this->num_rows); \
    } else { \
      vnl_c_vector<T>::deallocate(this->data, 1); \
    } \
  } \
} while (false)

//: Creates a matrix with given number of rows and columns.
// Elements are not initialized. O(m*n).

template<class T>
vnl_matrix<T>::vnl_matrix (unsigned rowz, unsigned colz)
{
  vnl_matrix_alloc_blah(rowz, colz);
}

//: Creates a matrix with given number of rows and columns, and initialize all elements to value. O(m*n).

template<class T>
vnl_matrix<T>::vnl_matrix (unsigned rowz, unsigned colz, T const& value)
{
  vnl_matrix_alloc_blah(rowz, colz);
  for (unsigned int i = 0; i < rowz; ++ i)
    for (unsigned int j = 0; j < colz; ++ j)
      this->data[i][j] = value;
}

//: r rows, c cols, special type.  Currently implements "identity".
template <class T>
vnl_matrix<T>::vnl_matrix(unsigned r, unsigned c, vnl_matrix_type t)
{
  vnl_matrix_alloc_blah(r, c);
  if (t == vnl_matrix_identity) {
    assert(r == c);
    for (unsigned int i = 0; i < r; ++ i)
      for (unsigned int j = 0; j < c; ++ j)
        this->data[i][j] = (i==j) ? T(1) : T(0);
  }
}

#if 1 // fsm: who uses this?
//: Creates a matrix with given dimension (rows, cols) and initialize first n elements, row-wise, to values. O(m*n).

template<class T>
vnl_matrix<T>::vnl_matrix (unsigned rowz, unsigned colz, unsigned n, T const values[])
{
  vnl_matrix_alloc_blah(rowz, colz);
  if (n > rowz*colz)
    n = rowz*colz;
  T *dst = this->data[0];
  for (unsigned k=0; k<n; ++k)
    dst[k] = values[k];
}
#endif

//: Creates a matrix from a block array of data, stored row-wise.
// O(m*n).

template<class T>
vnl_matrix<T>::vnl_matrix (T const* datablck, unsigned rowz, unsigned colz)
{
  vnl_matrix_alloc_blah(rowz, colz);
  unsigned int n = rowz*colz;
  T *dst = this->data[0];
  for (unsigned int k=0; k<n; ++k)
    dst[k] = datablck[k];
}


//: Creates a new matrix and copies all the elements.
// O(m*n).

template<class T>
vnl_matrix<T>::vnl_matrix (vnl_matrix<T> const& from)
{
  if (from.data) {
    vnl_matrix_alloc_blah(from.num_rows, from.num_cols);
    unsigned int n = this->num_rows * this->num_cols;
    T *dst = this->data[0];
    T const *src = from.data[0];
    for (unsigned int k=0; k<n; ++k)
      dst[k] = src[k];
  }
  else {
    num_rows = 0;
    num_cols = 0;
    data = 0;
  }
}

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

template<class T>
vnl_matrix<T>::vnl_matrix (vnl_matrix<T> const &A, vnl_matrix<T> const &B, vnl_tag_add)
{
#ifndef NDEBUG
  if (A.num_rows != B.num_rows || A.num_cols != B.num_cols)
    vnl_error_matrix_dimension ("vnl_tag_add", A.num_rows, A.num_cols, B.num_rows, B.num_cols);
#endif

  vnl_matrix_alloc_blah(A.num_rows, A.num_cols);

  unsigned int n = A.num_rows * A.num_cols;
  T const *a = A.data[0];
  T const *b = B.data[0];
  T *dst = this->data[0];

  for (unsigned int i=0; i<n; ++i)
    dst[i] = a[i] + b[i];
}

template<class T>
vnl_matrix<T>::vnl_matrix (vnl_matrix<T> const &A, vnl_matrix<T> const &B, vnl_tag_sub)
{
#ifndef NDEBUG
  if (A.num_rows != B.num_rows || A.num_cols != B.num_cols)
    vnl_error_matrix_dimension ("vnl_tag_sub", A.num_rows, A.num_cols, B.num_rows, B.num_cols);
#endif

  vnl_matrix_alloc_blah(A.num_rows, A.num_cols);

  unsigned int n = A.num_rows * A.num_cols;
  T const *a = A.data[0];
  T const *b = B.data[0];
  T *dst = this->data[0];

  for (unsigned int i=0; i<n; ++i)
    dst[i] = a[i] - b[i];
}

template<class T>
vnl_matrix<T>::vnl_matrix (vnl_matrix<T> const &M, T s, vnl_tag_mul)
{
  vnl_matrix_alloc_blah(M.num_rows, M.num_cols);

  unsigned int n = M.num_rows * M.num_cols;
  T const *m = M.data[0];
  T *dst = this->data[0];

  for (unsigned int i=0; i<n; ++i)
    dst[i] = m[i] * s;
}

template<class T>
vnl_matrix<T>::vnl_matrix (vnl_matrix<T> const &M, T s, vnl_tag_div)
{
  vnl_matrix_alloc_blah(M.num_rows, M.num_cols);

  unsigned int n = M.num_rows * M.num_cols;
  T const *m = M.data[0];
  T *dst = this->data[0];

  for (unsigned int i=0; i<n; ++i)
    dst[i] = m[i] / s;
}

template<class T>
vnl_matrix<T>::vnl_matrix (vnl_matrix<T> const &M, T s, vnl_tag_add)
{
  vnl_matrix_alloc_blah(M.num_rows, M.num_cols);

  unsigned int n = M.num_rows * M.num_cols;
  T const *m = M.data[0];
  T *dst = this->data[0];

  for (unsigned int i=0; i<n; ++i)
    dst[i] = m[i] + s;
}

template<class T>
vnl_matrix<T>::vnl_matrix (vnl_matrix<T> const &M, T s, vnl_tag_sub)
{
  vnl_matrix_alloc_blah(M.num_rows, M.num_cols);

  unsigned int n = M.num_rows * M.num_cols;
  T const *m = M.data[0];
  T *dst = this->data[0];

  for (unsigned int i=0; i<n; ++i)
    dst[i] = m[i] - s;
}

template<class T>
vnl_matrix<T>::vnl_matrix (vnl_matrix<T> const &A, vnl_matrix<T> const &B, vnl_tag_mul)
{
#ifndef NDEBUG
  if (A.num_cols != B.num_rows)
    vnl_error_matrix_dimension("vnl_tag_mul", A.num_rows, A.num_cols, B.num_rows, B.num_cols);
#endif

  unsigned int l = A.num_rows;
  unsigned int m = A.num_cols; // == B.num_rows
  unsigned int n = B.num_cols;

  vnl_matrix_alloc_blah(l, n);

  for (unsigned int i=0; i<l; ++i) {
    for (unsigned int k=0; k<n; ++k) {
      T sum(0);
      for (unsigned int j=0; j<m; ++j)
        sum += A.data[i][j] * B.data[j][k];
      this->data[i][k] = sum;
    }
  }
}

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

//: Frees up the dynamic storage used by matrix.
// O(m*n).

template<class T>
void vnl_matrix<T>::destroy()
{
  vnl_matrix_free_blah;
}

template<class T>
void vnl_matrix<T>::clear()
{
  if (data) {
    destroy();
    num_rows = 0;
    num_cols = 0;
    data = 0;
  }
}

// resize // Resizes the data arrays of THIS matrix to (rows x cols). O(m*n).
// Elements are not initialized. Returns true if size is changed.

template<class T>
bool vnl_matrix<T>::make_size (unsigned rowz, unsigned colz)
{
  if (this->data) {
    // if no change in size, do not reallocate.
    if (this->num_rows == rowz && this->num_cols == colz)
      return false;

    // else, simply release old storage and allocate new.
    vnl_matrix_free_blah;
    vnl_matrix_alloc_blah(rowz, colz);
  }
  else {
    // This happens if the matrix is default constructed.
    vnl_matrix_alloc_blah(rowz, colz);
  }

  return true;
}

#undef vnl_matrix_alloc_blah
#undef vnl_matrix_free_blah

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

//: Sets all elements of matrix to specified value. O(m*n).

template<class T>
void vnl_matrix<T>::fill (T const& value) {
  for (unsigned int i = 0; i < this->num_rows; i++)