vnl_matrix_fixed.h

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

// This is vxl/vnl/vnl_matrix_fixed.h
#ifndef vnl_matrix_fixed_h_
#define vnl_matrix_fixed_h_
#ifdef VCL_NEEDS_PRAGMA_INTERFACE
#pragma interface
#endif

//:
// \file
// \brief fixed size matrix
//
// \author Andrew W. Fitzgibbon, Oxford RRG
// \date   04 Aug 96
//
// \verbatim
// Modifications:
// Peter Vanroose, 23 Nov 1996:  added explicit copy constructor
// LSB (Manchester) 15/03/2001:  added Binary I/O and tidied up the documentation
//   Feb.2002 - Peter Vanroose - brief doxygen comment placed on single line
//   Oct.2002 - Amitha Perera  - separated vnl_matrix and vnl_matrix_fixed,
//                               removed necessity for vnl_matrix_fixed_ref
// 26Oct.2002 - Peter Vanroose - added inplace_transpose() method
// \endverbatim
//-----------------------------------------------------------------------------

#include <vcl_cassert.h>

#include "vnl_matrix.h"
#include "vnl_matrix_ref.h"
#include "vnl_vector.h"
#include "vnl_vector_fixed.h"

// This mess is for a MSVC6 workaround.
//
// The problem: the matrix-matrix operator* should be written as a
// non-member function since vxl (currently) forbits the use of member
// templates. However, when declared as
//
//     template<typename T, unsigned m, unsigned n, unsigned o>
//     matrix<T,m,o> operator*( matrix<T,m,n>, matrix<T,n,o> );
//
// MSVC6 does not find it. A solution is to declare it as a member
// template. However, the obvious
//
//     template<unsigned o>
//     matrix<T,num_rows,o> operator*( matrix<T,num_cols,o> );
//
// causes an internal compiler error. It turns out that if the new
// template parameter "o" comes _first_, then all is okay. Now, we
// can't change the signature of vnl_matrix_fixed to <unsigned cols,
// unsigned rows, type>, so we use a "hidden" helper matrix. Except
// that user defined conversion operators and conversion constructors
// are not called for templated functions. So we have to use a helper
// base class. The base class is empty, which means that there is no
// loss in space or time efficiency. Finally, we have:
//
//   template<unsigned cols, unsigned rows, typename T>
//   class fake_base { };
//
//   template<typename T, unsigned rows, unsigned cols>
//   class matrix : public fake_base<cols,rows,T>
//   {
//      template<unsigned o>
//      matrix<T,rows,o>  operator*( fake_base<o,cols,T> );
//   };
//
// Notice how "o" is first in the list of template parameters. Since
// base class conversions _are_ performed during template matching,
// matrix<T,m,n> is matched as fake_base<n,m,T>, and all is good. For
// some values of good.
//
// Of course, all this trickery is pre-processed away for conforming
// compilers.
//
template <class T, unsigned int num_rows, unsigned int num_cols>
class vnl_matrix_fixed;
template <class T, unsigned M, unsigned N>
inline
vnl_vector_fixed<T, M> vnl_matrix_fixed_mat_vec_mult(const vnl_matrix_fixed<T, M, N>& a, const vnl_vector_fixed<T, N>& b);
template <class T, unsigned M, unsigned N, unsigned O>
inline
vnl_matrix_fixed<T, M, O> vnl_matrix_fixed_mat_mat_mult(const vnl_matrix_fixed<T, M, N>& a, const vnl_matrix_fixed<T, N, O>& b);
#if VCL_VC60
template<unsigned cols, unsigned rows, typename T>
class vnl_matrix_fixed_fake_base
{
};
#define VNL_MATRIX_FIXED_VCL60_WORKAROUND : public vnl_matrix_fixed_fake_base<num_cols,num_rows,T>
#else
#define VNL_MATRIX_FIXED_VCL60_WORKAROUND /* no workaround. Phew. */
#endif


//: Fixed size, stack-stored, space-efficient matrix.
// vnl_matrix_fixed is a fixed-length, stack storage vector. It has
// the same storage size as a C-style array. It is not related via
// inheritance to vnl_matrix. However, it can be converted cheaply to
// a vnl_matrix_ref.
//
// Read the overview documentation of vnl_vector_fixed. The text there
// applies here.
template <class T, unsigned int num_rows, unsigned int num_cols>
class vnl_matrix_fixed  VNL_MATRIX_FIXED_VCL60_WORKAROUND
{
public:
  typedef unsigned int size_type;

private:
  T data_[num_rows][num_cols]; // Local storage

public:

  //: Construct an empty num_rows*num_cols matrix
  vnl_matrix_fixed() {}

  //: Construct an m*n matrix and fill with value
  explicit vnl_matrix_fixed(T value)
  {
    T* p = data_[0];
    unsigned int n = num_rows * num_cols;
    while (n--)
      *p++ = value;
  }

  //: Construct an m*n Matrix and copy data into it row-wise.
  explicit vnl_matrix_fixed(const T* datablck)
  {
    vcl_memcpy(data_[0], datablck, num_rows*num_cols*sizeof(T));
  }

  //: Construct an m*n Matrix and copy rhs into it.
  //  Abort if rhs is not the same size.
  vnl_matrix_fixed(const vnl_matrix_fixed& rhs)
  {
    vcl_memcpy(data_[0], rhs.data_block(), num_rows*num_cols*sizeof(T));
  }

  //: Construct an m*n Matrix and copy rhs into it.
  //  Abort if rhs is not the same size.
  vnl_matrix_fixed(const vnl_matrix<T>& rhs)
  {
    assert(rhs.rows() == num_rows && rhs.columns() == num_cols);
    vcl_memcpy(data_[0], rhs.data_block(), num_rows*num_cols*sizeof(T));
  }

  //  Destruct the m*n matrix.
  // An explicit destructor seems to be necessary, at least for gcc 3.0.0,
  // to avoid the compiler generating multiple versions of it.
  // (This way, a weak symbol is generated; otherwise not.  A bug of gcc 3.0.)
  ~vnl_matrix_fixed() {}

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

  //: Copy a vnl_matrix into this.
  //  Abort if rhs is not the same size.
  vnl_matrix_fixed& operator=(const vnl_matrix<T>& rhs)
  {
    assert(rhs.rows() == num_rows && rhs.columns() == num_cols);
    vcl_memcpy(data_[0], rhs.data_block(), num_rows*num_cols*sizeof(T));
    return *this;
  }

  //: Copy another vnl_matrix_fixed<T,m,n> into this.
  vnl_matrix_fixed& operator=(const vnl_matrix_fixed& rhs)
  {
    vcl_memcpy(data_[0], rhs.data_block(), num_rows*num_cols*sizeof(T));
    return *this;
  }

// 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
  void put (unsigned r, unsigned c, T const& v) { (*this)(r,c) = v; }

  //: get element
  T    get (unsigned r, unsigned c) const { return (*this)(r,c); }

  //: 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);

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

  //: 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;

  //: Transpose this matrix efficiently, if it is a square matrix
  void inplace_transpose();


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

  //: Add \a s to each element of lhs matrix in situ
  vnl_matrix_fixed& operator+= (T s)
  {
    add( data_block(), s, data_block() ); return *this;
  }

  //: Subtract \a s from each element of lhs matrix in situ
  vnl_matrix_fixed& operator-= (T s)
  {
    sub( data_block(), s, data_block() ); return *this;
  }

  //:
  vnl_matrix_fixed& operator*= (T s)
  {
    mul( data_block(), s, data_block() ); return *this;
  }

  //:
  vnl_matrix_fixed& operator/= (T s)
  {
    div( data_block(), s, data_block() ); return *this;
  }

  //:
  vnl_matrix_fixed& operator+= (vnl_matrix_fixed const& m)
  {
    add( data_block(), m.data_block(), data_block() ); return *this;
  }

  //:
  vnl_matrix_fixed& operator+= (vnl_matrix<T> const& m)
  {
    assert( m.rows() == rows() && m.cols() == cols() );
    add( data_block(), m.data_block(), data_block() ); return *this;
  }

  //:
  vnl_matrix_fixed& operator-= (vnl_matrix_fixed const& m)
  {
    sub( data_block(), m.data_block(), data_block() ); return *this;
  }

  //:
  vnl_matrix_fixed& operator-= (vnl_matrix<T> const& m)
  {
    assert( m.rows() == rows() && m.cols() == cols() );
    sub( data_block(), m.data_block(), data_block() ); return *this;
  }

  //: Negate all elements of matrix
  vnl_matrix_fixed operator- () const
  {
    vnl_matrix_fixed r;
    sub( T(0), data_block(), r.data_block() );
    return r;
  }

#if VCL_VC60
  template<unsigned o>
  vnl_matrix_fixed<T,num_rows,o> operator*( vnl_matrix_fixed_fake_base<o,num_cols,T> const& mat ) const
  {
    vnl_matrix_fixed<T,num_cols,o> const& b = static_cast<vnl_matrix_fixed<T,num_cols,o> const&>(mat);
    return vnl_matrix_fixed_mat_mat_mult<T,num_rows,num_cols,o>( *this, b );
  }
  vnl_vector_fixed<T, num_rows> operator*( vnl_vector_fixed<T, num_cols> const& b) const
  {
    return vnl_matrix_fixed_mat_vec_mult<T,num_rows,num_cols>(*this,b);
  }
#endif

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

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

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

  //: Return transpose
  vnl_matrix_fixed<T,num_cols,num_rows> transpose () const;

  //: Return conjugate transpose
  vnl_matrix_fixed<T,num_cols,num_rows> conjugate_transpose () const;

  //: Set values of this matrix to those of M, starting at [top,left]
  vnl_matrix_fixed& 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)
  void set_row   (unsigned i, T const * v);

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

  //: Set the i-th row
  void set_row   (unsigned i, vnl_vector<T> const&);

  //: Extract a sub-matrix of size rows x cols, starting at (top,left)
  //  Thus it contains elements  [top,top+rows-1][left,left+cols-1]
  vnl_matrix<T> extract (unsigned rows,  unsigned cols,
                         unsigned top=0, unsigned left=0) const;

  //: Get a vector equal to the given row
  vnl_vector<T> get_row   (unsigned row) const;

  //: Get a vector equal to the given column
  vnl_vector<T> get_column(unsigned col) const;

  //: Get n rows beginning at rowstart
  vnl_matrix<T> get_n_rows   (unsigned rowstart, unsigned n) const;

  //: Get n columns beginning at colstart
  vnl_matrix<T> get_n_columns(unsigned colstart, unsigned n) const;


  // mutators

  //: Set this matrix to an identity matrix
  //  Abort if the matrix is not square
  void set_identity();

  //: Reverse order of rows.
  void flipud();

  //: Reverse order of columns.
  void fliplr();

  //: Normalize each row so it is a unit vector
  //  Zero rows are ignored
  void normalize_rows();

  //: Normalize each column so it is a unit vector
  //  Zero columns are ignored
  void normalize_columns();

  //: Scale elements in given row by a factor of T
  void scale_row   (unsigned row, T value);

  //: Scale elements in given column by a factor of T
  void scale_column(unsigned col, T value);

  //: Type def for norms.
  typedef typename vnl_c_vector<T>::abs_t abs_t;

  //: Return sum of absolute values of elements
  abs_t array_one_norm() const { return vnl_c_vector<T>::one_norm(begin(), size()); }

  //: Return square root of sum of squared absolute element values
  abs_t array_two_norm() const { return vnl_c_vector<T>::two_norm(begin(), size()); }

  //: Return largest absolute element value
  abs_t array_inf_norm() const { return vnl_c_vector<T>::inf_norm(begin(), size()); }

  //: Return sum of absolute values of elements
  abs_t absolute_value_sum() const { return array_one_norm(); }