vnl_matrix.txx

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

    for (unsigned int j = 0; j < this->num_cols; j++)
      this->data[i][j] = value;
}

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

template<class T>
void vnl_matrix<T>::fill_diagonal (T const& value) {
  for (unsigned int i = 0; i < this->num_rows && i < this->num_cols; i++)
    this->data[i][i] = value;
}

//: Assigns value to all elements of a matrix. O(m*n).

//template<class T>
//vnl_matrix<T>& vnl_matrix<T>::operator= (T const& value) {
//  for (unsigned i = 0; i < this->num_rows; i++)    // For each row in Matrix
//    for (unsigned j = 0; j < this->num_cols; j++)  // For each column in Matrix
//      this->data[i][j] = value;                 // Assign value
//  return *this;                                 // Return Matrix reference
//}


//: Copies all elements of rhs matrix into lhs matrix. O(m*n).
// If needed, the arrays in lhs matrix are freed up, and new arrays are
// allocated to match the dimensions of the rhs matrix.

template<class T>
vnl_matrix<T>& vnl_matrix<T>::operator= (vnl_matrix<T> const& rhs) {
  if (this != &rhs) { // make sure *this != m
    if (rhs.data) {
      this->make_size(rhs.num_rows, rhs.num_cols);
      for (unsigned int i = 0; i < this->num_rows; i++)
        for (unsigned int j = 0; j < this->num_cols; j++)
          this->data[i][j] = rhs.data[i][j];
    }
    else {
      // rhs is default-constructed.
      clear();
    }
  }
  return *this;
}

template<class T>
void vnl_matrix<T>::print(vcl_ostream& os) const {
  for (unsigned int i = 0; i < this->rows(); i++) {
    for (unsigned int j = 0; j < this->columns(); j++)
      os << this->data[i][j] << " ";
    os << "\n";
  }
}

//: Prints the 2D array of elements of a matrix out to a stream.
// O(m*n).

template<class T>
vcl_ostream& operator<< (vcl_ostream& os, vnl_matrix<T> const& m) {
  for (unsigned int i = 0; i < m.rows(); ++i) {
    for (unsigned int j = 0; j < m.columns(); ++j)
      os << m(i, j) << " ";
    os << vcl_endl;
  }
  return os;
}

//: Read an vnl_matrix from an ascii vcl_istream.
// Automatically determines file size if the input matrix has zero size.
template<class T>
vcl_istream& operator>>(vcl_istream& s, vnl_matrix<T>& M) {
  M.read_ascii(s);
  return s;
}

template<class T>
void vnl_matrix<T>::inline_function_tickler()
{
  vnl_matrix<T> M;
  // fsm: hack to get 2.96 to instantiate the inline function.
  M = T(1) + T(3) * M;
}

template<class T>
vnl_matrix<T>& vnl_matrix<T>::operator+= (T value) {
  for (unsigned int i = 0; i < this->num_rows; i++)
    for (unsigned int j = 0; j < this->num_cols; j++)
      this->data[i][j] += value;
  return *this;
}

template<class T>
vnl_matrix<T>& vnl_matrix<T>::operator-= (T value) {
  for (unsigned int i = 0; i < this->num_rows; i++)
    for (unsigned int j = 0; j < this->num_cols; j++)
      this->data[i][j] -= value;
  return *this;
}

template<class T>
vnl_matrix<T>& vnl_matrix<T>::operator*= (T value) {
  for (unsigned int i = 0; i < this->num_rows; i++)
    for (unsigned int j = 0; j < this->num_cols; j++)
      this->data[i][j] *= value;
  return *this;
}

template<class T>
vnl_matrix<T>& vnl_matrix<T>::operator/= (T value) {
  for (unsigned int i = 0; i < this->num_rows; i++)
    for (unsigned int j = 0; j < this->num_cols; j++)
      this->data[i][j] /= value;
  return *this;
}


//: Adds lhs matrix with rhs matrix, and stores in place in lhs matrix.
// O(m*n). The dimensions of the two matrices must be identical.

template<class T>
vnl_matrix<T>& vnl_matrix<T>::operator+= (vnl_matrix<T> const& rhs) {
#ifndef NDEBUG
  if (this->num_rows != rhs.num_rows ||
      this->num_cols != rhs.num_cols)           // Size match?
    vnl_error_matrix_dimension ("operator+=",
                                this->num_rows, this->num_cols,
                                rhs.num_rows, rhs.num_cols);
#endif
  for (unsigned int i = 0; i < this->num_rows; i++)    // For each row
    for (unsigned int j = 0; j < this->num_cols; j++)  // For each element in column
      this->data[i][j] += rhs.data[i][j];       // Add elements
  return *this;
}


//: Substract lhs matrix with rhs matrix and store in place in lhs matrix.
// O(m*n).
// The dimensions of the two matrices must be identical.

template<class T>
vnl_matrix<T>& vnl_matrix<T>::operator-= (vnl_matrix<T> const& rhs) {
#ifndef NDEBUG
  if (this->num_rows != rhs.num_rows ||
      this->num_cols != rhs.num_cols) // Size?
    vnl_error_matrix_dimension ("operator-=",
                                this->num_rows, this->num_cols,
                                rhs.num_rows, rhs.num_cols);
#endif
  for (unsigned int i = 0; i < this->num_rows; i++)
    for (unsigned int j = 0; j < this->num_cols; j++)
      this->data[i][j] -= rhs.data[i][j];
  return *this;
}


template<class T>
vnl_matrix<T> operator- (T const& value, vnl_matrix<T> const& m) {
  vnl_matrix<T> result(m.rows(),m.columns());
  for (unsigned int i = 0; i < m.rows(); i++)  // For each row
    for (unsigned int j = 0; j < m.columns(); j++) // For each elmt. in column
      result.put(i,j, value - m.get(i,j) );     // subtract from value elmt.
  return result;
}


#if 0 // commented out
//: Returns new matrix which is the product of m1 with m2, m1 * m2.
// O(n^3). Number of columns of first matrix must match number of rows
// of second matrix.

template<class T>
vnl_matrix<T> vnl_matrix<T>::operator* (vnl_matrix<T> const& rhs) const {
#ifndef NDEBUG
  if (this->num_cols != rhs.num_rows)           // dimensions do not match?
    vnl_error_matrix_dimension("operator*",
                               this->num_rows, this->num_cols,
                               rhs.num_rows, rhs.num_cols);
#endif
  vnl_matrix<T> result(this->num_rows, rhs.num_cols); // Temp to store product
  for (unsigned i = 0; i < this->num_rows; i++) {  // For each row
    for (unsigned j = 0; j < rhs.num_cols; j++) {  // For each element in column
      T sum = 0;
      for (unsigned k = 0; k < this->num_cols; k++) // Loop over column values
        sum += (this->data[i][k] * rhs.data[k][j]);     // Multiply
      result(i,j) = sum;
    }
  }
  return result;
}
#endif

//: Returns new matrix which is the negation of THIS matrix.
// O(m*n).

template<class T>
vnl_matrix<T> vnl_matrix<T>::operator- () const {
  vnl_matrix<T> result(this->num_rows, this->num_cols);
  for (unsigned int i = 0; i < this->num_rows; i++)
    for (unsigned int j = 0; j < this->num_cols; j++)
      result.data[i][j] = - this->data[i][j];
  return result;
}

#if 0 // commented out
//: Returns new matrix with elements of lhs matrix added with value.
// O(m*n).

template<class T>
vnl_matrix<T> vnl_matrix<T>::operator+ (T const& value) const {
  vnl_matrix<T> result(this->num_rows, this->num_cols);
  for (unsigned i = 0; i < this->num_rows; i++)    // For each row
    for (unsigned j = 0; j < this->num_cols; j++)  // For each element in column
      result.data[i][j] = (this->data[i][j] + value);   // Add scalar
  return result;
}


//: Returns new matrix with elements of lhs matrix multiplied with value.
// O(m*n).

template<class T>
vnl_matrix<T> vnl_matrix<T>::operator* (T const& value) const {
  vnl_matrix<T> result(this->num_rows, this->num_cols);
  for (unsigned i = 0; i < this->num_rows; i++)    // For each row
    for (unsigned j = 0; j < this->num_cols; j++)  // For each element in column
      result.data[i][j] = (this->data[i][j] * value);   // Multiply
  return result;
}


//: Returns new matrix with elements of lhs matrix divided by value. O(m*n).
template<class T>
vnl_matrix<T> vnl_matrix<T>::operator/ (T const& value) const {
  vnl_matrix<T> result(this->num_rows, this->num_cols);
  for (unsigned i = 0; i < this->num_rows; i++)    // For each row
    for (unsigned j = 0; j < this->num_cols; j++)  // For each element in column
      result.data[i][j] = (this->data[i][j] / value);   // Divide
  return result;
}
#endif

//: Return the matrix made by applying "f" to each element.
template <class T>
vnl_matrix<T> vnl_matrix<T>::apply(T (*f)(T const&)) const {
  vnl_matrix<T> ret(num_rows, num_cols);
  vnl_c_vector<T>::apply(this->data[0], num_rows * num_cols, f, ret.data_block());
  return ret;
}

//: Return the matrix made by applying "f" to each element.
template <class T>
vnl_matrix<T> vnl_matrix<T>::apply(T (*f)(T)) const {
  vnl_matrix<T> ret(num_rows, num_cols);
  vnl_c_vector<T>::apply(this->data[0], num_rows * num_cols, f, ret.data_block());
  return ret;
}

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

//: Returns new matrix with rows and columns transposed.
// O(m*n).

template<class T>
vnl_matrix<T> vnl_matrix<T>::transpose() const {
  vnl_matrix<T> result(this->num_cols, this->num_rows);
  for (unsigned int i = 0; i < this->num_cols; i++)
    for (unsigned int j = 0; j < this->num_rows; j++)
      result.data[i][j] = this->data[j][i];
  return result;
}

// adjoint/hermitian transpose

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

//: Replaces the submatrix of THIS matrix, starting at top left corner, by the elements of matrix m. O(m*n).
// This is the reverse of extract().

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


//: Returns a copy of submatrix of THIS matrix, specified by the top-left corner and size in rows, cols. O(m*n).
// Use update() to copy new values of this submatrix back into THIS matrix.

template<class T>
vnl_matrix<T> vnl_matrix<T>::extract (unsigned rowz, unsigned colz,
                                      unsigned top, unsigned left) const{
#ifndef NDEBUG
  unsigned int bottom = top + rowz;
  unsigned int right = left + colz;
  if ((this->num_rows < bottom) || (this->num_cols < right))
    vnl_error_matrix_dimension ("extract",
                                this->num_rows, this->num_cols, 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.data[i][j] = data[top+i][left+j];
  return result;
}

//: Returns the dot product of the two matrices. O(m*n).
// This is the sum of all pairwise products of the elements m1[i,j]*m2[i,j].

template<class T>
T dot_product (vnl_matrix<T> const& m1, vnl_matrix<T> const& m2) {
#ifndef NDEBUG
  if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) // Size?
    vnl_error_matrix_dimension ("dot_product",
                                m1.rows(), m1.columns(),
                                m2.rows(), m2.columns());
#endif
  return vnl_c_vector<T>::dot_product(m1.begin(), m2.begin(), m1.rows()*m1.cols());
}

//: Hermitian inner product.
// O(mn).

template<class T>
T inner_product (vnl_matrix<T> const& m1, vnl_matrix<T> const& m2) {
#ifndef NDEBUG
  if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) // Size?
    vnl_error_matrix_dimension ("inner_product",
                                m1.rows(), m1.columns(),
                                m2.rows(), m2.columns());
#endif
  return vnl_c_vector<T>::inner_product(m1.begin(), m2.begin(), m1.rows()*m1.cols());
}

// cos_angle. O(mn).

template<class T>
T cos_angle (vnl_matrix<T> const& a, vnl_matrix<T> const& b) {
  typedef typename vnl_numeric_traits<T>::real_t real_t;
  typedef typename vnl_numeric_traits<T>::abs_t abs_t;
  typedef typename vnl_numeric_traits<abs_t>::real_t abs_r;

  T ab = inner_product(a,b);
  abs_t a_b = (abs_t)vcl_sqrt( (abs_r)vnl_math_abs(inner_product(a,a) * inner_product(b,b)) );

  return T( ab / a_b);
}

//: Returns new matrix whose elements are the products m1[ij]*m2[ij].
// O(m*n).

template<class T>
vnl_matrix<T> element_product (vnl_matrix<T> const& m1,
                               vnl_matrix<T> const& m2) {
#ifndef NDEBUG
  if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) // Size?
    vnl_error_matrix_dimension ("element_product",
                                m1.rows(), m1.columns(), m2.rows(), m2.columns());
#endif
  vnl_matrix<T> result(m1.rows(), m1.columns());
  for (unsigned int i = 0; i < m1.rows(); i++)
    for (unsigned int j = 0; j < m1.columns(); j++)
      result.put(i,j, m1.get(i,j) * m2.get(i,j) );
  return result;
}