vnl_matrix_fixed_ref.h

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

// This is core/vnl/vnl_matrix_fixed_ref.h
#ifndef vnl_matrix_fixed_ref_h_
#define vnl_matrix_fixed_ref_h_
#ifdef VCL_NEEDS_PRAGMA_INTERFACE
#pragma interface
#endif
//:
//  \file
//  \brief Fixed size stack-stored vnl_matrix
//
//  > > From: Amitha Perera [mailto:[e-mail]]
//  > > Sent: Monday, October 07, 2002 3:18 PM
//  > > Subject: vnl_vector_fixed_ref
//  > >
//  > > I'm working on separating vnl_vector and vnl_vector_fixed in the VXL
//  > > tree, as I mailed a while ago to the vxl-maintainers list. I noticed
//  > > that you'd committed a vnl_vector_fixed_ref class which doesn't seem
//  > > to provide any additional functionality over vnl_vector_ref. May I
//  > > remove it, or is there some use for it?
//  > >
//  > > FYI, the author is listed as "Paul P. Smyth, Vicon Motion
//  > > Systems Ltd."
//  > > and the comment is dated 02 May 2001.
//
//  Paul Smyth <[e-mail]> writes:
//  > The rationale behind it was that I had some (fast) algorithms for
//  > matrix/vector operations that made use of compile-time knowledge of the
//  > vector and matrix sizes.
//  > This was typically appropriate when trying to interpret a fixed-size
//  > subvector within a large vector of parameters as e.g. a translation.
//  >
//  > As I saw it, the various types of vector possible were: (with their current
//  > names)
//  > - pointer to memory, plus compile-time knowledge of vector size ( T*, and enum{size}) = vnl_vector_fixed_ref
//  > - ownership of memory, plus compile-time size = vnl_vector_fixed
//  > - pointer to memory, plus run-time only knowledge of size (T* and size()) = vnl_vector_ref
//  > - ownership of memory, variably sized = vnl_vector
//  >
//  > I had a conversation with Andrew Fitzgibbon, where he reckoned that the best
//  > thing to do with vnl vectors etc. was to create entirely separate types, and
//  > routines for conversion between them (possibly implicitly), rather that
//  > trying to establish a class hierarchy, which may add too many burdens in
//  > terms of object size for small vectors/matrices.
//  >
//  > Sorry - I've now found the debate on the maintaners list!
//  >
//  > Anyway, I believe that vector_fixed_ref is very necessary, and that you
//  > should be able to convert from a vector_fixed to a vector_fixed_ref - say
//  > using an as_ref() member on vector_fixed or standalone function.
//  > And I believe that for the restructured classes, vector_fixed_ref and
//  > vector_fixed should not be related by inheritance, as that would place an
//  > excessive burden on the size of vector_fixed.
//  >
//  > ------
//  > Another issue - do you have a mechanism for dealing with const data safely?
//  > {
//  >   template<typename T, int n>
//  >   vnl_vector_fixed_ref(T* i_Data);
//  >
//  >   void MyFunction(const vnl_vector<double> & Input)
//  >   {
//  >     // take a reference to the first 3 elements of Input
//  >     vnl_vector_fixed_ref<double,3> ref(Input.begin());
//  >     // compiler error - as making vector_fixed_ref from const
//  > double *
//  >   }
//  > }
//  >
//  > The options appear to be
//  > 1) Make a separate class vnl_vector_fixed_ref_const
//  > 2) Make vnl_vector_fixed_ref so it can be instantiated with
//  > vnl_vector_fixed_ref<double,n> AND vnl_vector_fixed_ref<const double,n>, and
//  > gives appropriate behaviour - would probably require a to_const function
//  > which generates vnl_vector_fixed_ref<const T,n> from
//  > vnl_vector_fixed_ref<T,n>
//  >
//  > ------
//  > Another note is that a number of routines that use vector_fixed currently
//  > (e.g. cross_3d) should really use vector_fixed_ref as an input, because they
//  > should be able to operate on fixed vector references as well as fixed
//  > vectors.
//  >
//  > While I'm at it, has it been decided that the vnl_vector and vnl_vector_ref
//  > classes are to remain unchanged? Because having vnl_vector as the base, and
//  > vnl_vector_ref derived from it is a real pain in the backside. A vector
//  > which may or may not own its own memory is a more general type than one
//  > which does own it's own memory, and having vnl_vector as the base means that
//  > all sorts of nastinesses can happen. Simply, a vector_ref Is-not a type of
//  > vector.
//  > If anything, it should be the other way round.
//  >
//  > void DoAssign(vnl_vector<double> & RefToMemoryIDontOwn, const vnl_vector<double> & NewContents)
//  > {
//  >   RefToMemoryIDontOwn = NewContents;
//  > }
//  >
//  > void DeleteTwice()
//  > {
//  >   vnl_vector<double> vec1(3, 0); // size 3 - news 3*double
//  >   vnl_vector<double> vec2(4,1); // size 4 news 4 * double
//  >   vnl_vector_ref<double> ref_to_1(3,vec1.begin()); // copies pointer
//  >   DoAssign(ref_to_1, vec2); // deletes memory owned by 1, news 4 * double
//  >   // vec1 now points to deleted memory, and will crash when goes out of scope
//  > }
//  >
//  > Maybe that issue isn't on your agenda - but it's a bit of a disaster. I know
//  > that fixing this might break some code.
//  >
//  > ---------
//  > Sorry for rolling all these things into one - I'd be interested to know what
//  > you think. But please don't kill my vnl_vector_ref!
//  >
//  > Paul.
//
//    vnl_matrix_fixed_ref is a fixed-size vnl_matrix for which the data space
//    has been supplied externally.  This is useful for two main tasks:
//
//    (a) Treating some row-based "C" matrix as a vnl_matrix in order to
//    perform vnl_matrix operations on it.
//
//    (b) Declaring a vnl_matrix that uses entirely stack-based storage for the
//    matrix.
//
//    The big warning is that returning a vnl_matrix_fixed_ref pointer will free
//    non-heap memory if deleted through a vnl_matrix pointer.  This should be
//    very difficult though, as vnl_matrix_fixed_ref objects may not be constructed
//    using operator new.  This in turn is plausible as the point is to avoid
//    such calls.
//
// \author Andrew W. Fitzgibbon, Oxford RRG
// \date   04 Aug 96
//
// \verbatim
//  Modifications:
//   27-Nov-1996 Peter Vanroose - added default constructor which allocates matrix storage
//    4-Jul-2003 Paul Smyth - general cleanup and rewrite; interface now as vnl_matrix_fixed
//   15-Aug-2003 Peter Vanroose - removed "duplicate" operator=(vnl_matrix_fixed<T,n> const&)
//    8-Dec-2006 Markus Moll - changed operator>> signature (to const& argument)
// \endverbatim
//
//-----------------------------------------------------------------------------

#include <vcl_cassert.h>
#include <vcl_iosfwd.h>
#include <vcl_cstring.h> // memcpy()
#include <vnl/vnl_matrix_fixed.h>
#include <vnl/vnl_vector_fixed.h>
#include <vnl/vnl_vector_fixed_ref.h>
#include <vnl/vnl_c_vector.h>

//: Fixed size stack-stored vnl_matrix
// vnl_matrix_fixed_ref is a fixed-size vnl_matrix for which the data space
// has been supplied externally.  This is useful for two main tasks:
//
// (a) Treating some row-based "C" matrix as a vnl_matrix in order to
// perform vnl_matrix operations on it.
//
// (b) Declaring a vnl_matrix that uses entirely stack-based storage for the
// matrix.
//
// The big warning is that returning a vnl_matrix_fixed_ref pointer will free
// non-heap memory if deleted through a vnl_matrix pointer.  This should be
// very difficult though, as vnl_matrix_fixed_ref objects may not be constructed
// using operator new.  This in turn is plausible as the point is to avoid
// such calls.
//

template <class T, unsigned num_rows, unsigned num_cols>
class vnl_matrix_fixed_ref_const
{
 protected:
  const T* data_;
 public:
  vnl_matrix_fixed_ref_const(const vnl_matrix_fixed<T,num_rows,num_cols>& rhs)
    : data_(rhs.data_block())
  {
  }
  explicit vnl_matrix_fixed_ref_const(const T * dataptr)
    : data_(dataptr)
  {
  }
  vnl_matrix_fixed_ref_const(const vnl_matrix_fixed_ref_const<T,num_rows,num_cols> & rhs)
    : data_(rhs.data_)
  {
  }
  //: Get j-th row
  vnl_vector_fixed<T,num_rows> get_row(unsigned row_index) const
  {
    vnl_vector<T> v(num_cols);
    for (unsigned int j = 0; j < num_cols; j++)    // For each element in row
      v[j] = (*this)(row_index,j);
    return v;
  }

  //: Get j-th column
  vnl_vector_fixed<T,num_cols> get_column(unsigned column_index) const
  {
    vnl_vector<T> v(num_rows);
    for (unsigned int j = 0; j < num_rows; j++)
      v[j] = (*this)(j,column_index);
    return v;
  }
  const T * data_block() const { return data_; }

  //: Const iterators
  typedef T const *const_iterator;
  //: Iterator pointing to start of data
  const_iterator begin() const { return data_; }
  //: Iterator pointing to element beyond end of data
  const_iterator end() const { return begin() + this->size(); }

  //: Type defs for iterators
  typedef const T element_type;
  //: Type defs for iterators
  typedef const T       *iterator;

  T const & operator() (unsigned r, unsigned c) const
  {
#if VNL_CONFIG_CHECK_BOUNDS  && (!defined NDEBUG)
    assert(r<num_rows);   // Check the row index is valid
    assert(c<num_cols);   // Check the column index is valid
#endif
    return *(data_ + num_cols * r + c);
  }

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

  //: 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 num_rows*num_cols; }


  //: Print matrix to os in some hopefully sensible format
  void print(vcl_ostream& os) const;

  void copy_out(T *) const;

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

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

  //: Make a new matrix by applying function to each element.
  vnl_matrix_fixed<T,num_rows,num_cols> 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;

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

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

  //: Return largest absolute value
  abs_t absolute_value_max() const { return array_inf_norm(); }

  // $ || M ||_1 := \max_j \sum_i | M_{ij} | $
  abs_t operator_one_norm() const;

  // $ || M ||_\inf := \max_i \sum_j | M_{ij} | $
  abs_t operator_inf_norm() const;

  //: Return Frobenius norm of matrix (sqrt of sum of squares of its elements)
  abs_t frobenius_norm() const { return vnl_c_vector<T>::two_norm(begin(), size()); }

  //: Return Frobenius norm of matrix (sqrt of sum of squares of its elements)
  abs_t fro_norm() const { return frobenius_norm(); }

  //: Return RMS of all elements
  abs_t rms() const { return vnl_c_vector<T>::rms_norm(begin(), size()); }