vnl_matrix.h

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

// This is vxl/vnl/vnl_matrix.h
#ifndef vnl_matrix_h_
#define vnl_matrix_h_
#ifdef VCL_NEEDS_PRAGMA_INTERFACE
#pragma interface
#endif
//:
// \file
// \brief An ordinary mathematical matrix

#include <vcl_iosfwd.h>
#include <vnl/vnl_tag.h>
#include <vnl/vnl_c_vector.h>
#ifndef NDEBUG
# include "vnl_config.h"
# if VNL_CONFIG_CHECK_BOUNDS
#  include <vnl/vnl_error.h>
#  include <vcl_cassert.h>
# endif
#else
# undef VNL_CONFIG_CHECK_BOUNDS
# define VNL_CONFIG_CHECK_BOUNDS 0
# undef ERROR_CHECKING
#endif

#include <vxl_config_itk.h>

export template <class T> class vnl_vector;
export template <class T> class vnl_matrix;

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

#define v vnl_vector<T>
#define m vnl_matrix<T>
template <class T> m operator+ (T const&, m const&);
template <class T> m operator- (T const&, m const&);
template <class T> m operator* (T const&, m const&);
template <class T> m element_product(m const&, m const&);
template <class T> m element_quotient(m const&, m const&);
template <class T> T dot_product (m const&, m const&);
template <class T> T inner_product (m const&, m const&);
template <class T> T cos_angle(m const&, m const& );
template <class T> vcl_ostream& operator<< (vcl_ostream&, m const&);
template <class T> vcl_istream& operator>> (vcl_istream&, m&);
#undef v
#undef m

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

enum vnl_matrix_type
{
  vnl_matrix_null,
  vnl_matrix_identity
};

//:  An ordinary mathematical matrix
// The vnl_matrix<T> class implements two-dimensional arithmetic
// matrices  for  a user-specified numeric data type. Using the
// parameterized types facility of C++,  it  is  possible,  for
// example, for the user to create a matrix of rational numbers
// by parameterizing the vnl_matrix class over the Rational  class.
// The  only  requirement  for the type is that it supports the
// basic arithmetic operators.
//
// Note: Unlike   the   other   sequence   classes,   the
// vnl_matrix<T>  class is fixed-size. It will not grow once the
// size has been specified to the constructor or changed by the
// assignment  or  multiplication  operators.  The vnl_matrix<T>
// class is row-based with addresses of rows being cached,  and
// elements accessed as m[row][col].
//
// Note: The matrix can, however, be resized using the resize(nr,nc) function.
//
// Note: Indexing of the matrix is zero-based, so the top-left element is M(0,0).
//
// Note: Inversion of matrix M, and other operations such as solving systems of linear
// equations are handled by the matrix decomposition classes in vnl/algo, such
// as matrix_inverse, svd, qr etc.
//
// Note: Use a vnl_vector<T> with these matrices.

#ifdef VXL_CXX_SLICED_DESTRUCTOR_BUG
class vnl_matrix_base_own
{
public:
  vnl_matrix_base_own(): vnl_matrix_own_data(1) {}
protected:
  int vnl_matrix_own_data;
};

template<class T>
class vnl_matrix: public vnl_matrix_base_own
#else
template<class T>
class vnl_matrix
#endif
{
 public:
  //: Default constructor creates an empty matrix of size 0,0.
  vnl_matrix () :
    num_rows(0),
    num_cols(0),
    data(0)
  {
  }

  //: Construct a matrix of size r rows by c columns
  // Contents are unspecified.
  // Complexity $O(1)$
  vnl_matrix(unsigned r, unsigned c);                           // r rows, c cols.

  //: Construct a matrix of size r rows by c columns, and all emelemnts equal to v0
  // Complexity $O(r.c)$
  vnl_matrix(unsigned r, unsigned c, T const& v0);              // r rows, c cols, value v0.

  //: Construct a matrix of size r rows by c columns, with a special type
  // Contents are specified by t
  // Complexity $O(r.c)$
  vnl_matrix(unsigned r, unsigned c, vnl_matrix_type t);        // r rows, c cols, special type

  //: Construct a matrix of size r rows by c columns, initialised by an automatic array
  // The first n elements, are initialised row-wise, to values.
  // Complexity $O(n)$
  vnl_matrix(unsigned r, unsigned c, unsigned n, T const values[]);  // use automatic arrays.

  //: Construct a matrix of size r rows by c columns, initialised by a memory block
  // The values are initialise row wise from the data.
  // Complexity $O(r.c)$
  vnl_matrix(T const* data_block, unsigned r, unsigned c);      // fill row-wise.

  //: Copy construct a matrix
  // Complexity $O(r.c)$
  vnl_matrix(vnl_matrix<T> const&);                             // from another matrix.

#ifndef VXL_DOXYGEN_SHOULD_SKIP_THIS
// <internal>
  // These constructors are here so that operator* etc can take
  // advantage of the C++ return value optimization.
  vnl_matrix (vnl_matrix<T> const &, vnl_matrix<T> const &, vnl_tag_add); // M + M
  vnl_matrix (vnl_matrix<T> const &, vnl_matrix<T> const &, vnl_tag_sub); // M - M
  vnl_matrix (vnl_matrix<T> const &, T,                     vnl_tag_mul); // M * s
  vnl_matrix (vnl_matrix<T> const &, T,                     vnl_tag_div); // M / s
  vnl_matrix (vnl_matrix<T> const &, T,                     vnl_tag_add); // M + s
  vnl_matrix (vnl_matrix<T> const &, T,                     vnl_tag_sub); // M - s
  vnl_matrix (vnl_matrix<T> const &, vnl_matrix<T> const &, vnl_tag_mul); // M * M
  vnl_matrix (vnl_matrix<T> &that, vnl_tag_grab)
    : num_rows(that.num_rows), num_cols(that.num_cols), data(that.data)
  { that.num_cols=that.num_rows=0; that.data=0; } // "*this" now uses "that"'s data.
// </internal>
#endif

  //: Matrix destructor
  ~vnl_matrix()
  {
    // save some fcalls if data is 0 (i.e. in matrix_fixed)
#ifdef VXL_CXX_SLICED_DESTRUCTOR_BUG
    if (data && vnl_matrix_own_data) destroy();
#else
    if (data) destroy();
#endif
  }

// Basic 2D-Array functionality-------------------------------------------

  //: Return number of rows
  unsigned rows ()    const { return num_rows; }

  //: Return number of columns
   // A synonym for cols()
 unsigned columns ()  const { return num_cols; }

  //: Return number of columns
  // A synonym for columns()
  unsigned cols ()    const { return num_cols; }

  //: Return number of elements
  // This equals rows() * cols()
  unsigned size ()    const { return rows()*cols(); }

  //: set element with boundary checks if error checking is on.
  void put (unsigned r, unsigned c, T const&);

  //: get element with boundary checks if error checking is on.
  T    get (unsigned r, unsigned c) const;

  //: return pointer to given row
  // No boundary checking here.
  T       * operator[] (unsigned r) { return data[r]; }

  //: return pointer to given row
  // No boundary checking here.
  T const * operator[] (unsigned r) const { return data[r]; }

  //: Access an element for reading or writing
  // There are assert style boundary checks - #define NDEBUG to turn them off.
  T       & operator() (unsigned r, unsigned c)
  {
#if VNL_CONFIG_CHECK_BOUNDS  && (!defined NDEBUG)
    assert(r<rows());   // Check the row index is valid
    assert(c<cols());   // Check the column index is valid
#endif
    return this->data[r][c];
  }

  //: Access an element for reading
  // There are assert style boundary checks - #define NDEBUG to turn them off.
  T const & operator() (unsigned r, unsigned c) const
  {
#if VNL_CONFIG_CHECK_BOUNDS  && (!defined NDEBUG)
    assert(r<rows());   // Check the row index is valid
    assert(c<cols());   // Check the column index is valid
#endif
    return this->data[r][c];
  }


// Filling and copying------------------------------------------------

  //: Set all elements of matrix to specified value.
  // Complexity $O(r.c)$
  void fill (T const&);

  //: Set all diagonal elements of matrix to specified value.
  // Complexity $O(\min(r,c))$
  void fill_diagonal (T const&);

  //: Fill (laminate) this matrix with the given data.
  // We assume that p points to a contiguous rows*cols array, stored rowwise.
  void copy_in(T const *);

  //: Fill (laminate) this matrix with the given data.
  // A synonym for copy_in()
  void set(T const *d) { copy_in(d); }

  //: Fill the given array with this matrix.
  // We assume that p points to
  // a contiguous rows*cols array, stored rowwise.
  // No bounds checking on the array
  void copy_out(T *) const;


  //: Set all elements to value v
  // Complexity $O(r.c)$
  vnl_matrix<T>& operator= (T const&v) { fill(v); return *this; }

  //: Copies all elements of rhs matrix into lhs matrix.
  // Complexity $O(\min(r,c))$
  vnl_matrix<T>& operator= (vnl_matrix<T> const&);

// Arithmetic ----------------------------------------------------
  // note that these functions should not pass scalar as a const&.
  // Look what would happen to A /= A(0,0).

  //: Add rhs to each element of lhs matrix in situ
  vnl_matrix<T>& operator+= (T value);

  //: Subtract rhs from each element of lhs matrix in situ
  vnl_matrix<T>& operator-= (T value);

  //: Scalar multiplication in situ of lhs matrix  by rhs
  vnl_matrix<T>& operator*= (T value);

  //: Scalar division of lhs matrix  in situ by rhs
  vnl_matrix<T>& operator/= (T value);

  //: Add rhs to lhs  matrix in situ
  vnl_matrix<T>& operator+= (vnl_matrix<T> const&);
  //: Subtract rhs from lhs matrix in situ
  vnl_matrix<T>& operator-= (vnl_matrix<T> const&);
  //: Multiply lhs matrix in situ by rhs
  vnl_matrix<T>& operator*= (vnl_matrix<T> const&rhs) { *this = (*this) * rhs; return *this; }

  //: Negate all elements of matrix
  vnl_matrix<T> operator- () const;


  //: Add rhs to each element of lhs matrix and return result in new matrix
  vnl_matrix<T> operator+ (T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_add()); }

  //: Subtract rhs from each element of lhs matrix and return result in new matrix
  vnl_matrix<T> operator- (T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_sub()); }

  //: Scalar multiplication of lhs matrix by rhs  and return result in new matrix
  vnl_matrix<T> operator* (T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_mul()); }

  //: Scalar division of lhs matrix by rhs and return result in new matrix
  vnl_matrix<T> operator/ (T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_div()); }

  //: Matrix add rhs to lhs matrix and return result in new matrix
  vnl_matrix<T> operator+ (vnl_matrix<T> const& rhs) const { return vnl_matrix<T>(*this, rhs, vnl_tag_add()); }
  //: Matrix subtract rhs from lhs and return result in new matrix
  vnl_matrix<T> operator- (vnl_matrix<T> const& rhs) const { return vnl_matrix<T>(*this, rhs, vnl_tag_sub()); }
  //: Matrix multiply lhs by rhs matrix and return result in new matrix
  vnl_matrix<T> operator* (vnl_matrix<T> const& rhs) const { return vnl_matrix<T>(*this, rhs, vnl_tag_mul()); }

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

  //: Make a new matrix by applying function to each element.
  vnl_matrix<T> apply(T (*f)(T)) const;

  //: Make a new matrix by applying function to each element.
  vnl_matrix<T> apply(T (*f)(T const&)) const;

  //: Return transpose
  vnl_matrix<T> transpose () const;

  //: Return conjugate transpose
  vnl_matrix<T> conjugate_transpose () const;

  //: Set values of this matrix to those of M, starting at [top,left]
  vnl_matrix<T>& update (vnl_matrix<T> const&, unsigned top=0, unsigned left=0);

  //: Set the elements of the i'th column to v[j]  (No bounds checking)
  void set_column(unsigned i, T const * v);

  //: Set the elements of the i'th column to value
  void set_column(unsigned i, T value );

  //: Set j-th colum to v
  void set_column(unsigned j, vnl_vector<T> const& v);

  //: Set columns to those in M, starting at starting_column
  void set_columns(unsigned starting_column, vnl_matrix<T> const& M);

  //: Set the elements of the i'th row to v[j]  (No bounds checking)