Merging upstream version 1.44.3.

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 1874 insertions, 0 deletions

diff --git a/ml/dlib/dlib/matrix/matrix_utilities_abstract.h b/ml/dlib/dlib/matrix/matrix_utilities_abstract.h
new file mode 100644
index 000000000..ad4c91167
--- /dev/null
+++ b/ml/dlib/dlib/matrix/matrix_utilities_abstract.h
@@ -0,0 +1,1874 @@
+// Copyright (C) 2006  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_MATRIx_UTILITIES_ABSTRACT_
+#ifdef DLIB_MATRIx_UTILITIES_ABSTRACT_
+
+#include "matrix_abstract.h"
+#include <complex>
+#include "../pixel.h"
+#include "../geometry/rectangle.h"
+#inclue <vector>
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+//                                   Simple matrix utilities 
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+    template <typename EXP>
+    constexpr bool is_row_major (
+        const matrix_exp<EXP>&
+    );
+    /*!
+        ensures
+            - returns true if and only if the given matrix expression uses the row_major_layout.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp diag (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns a column vector R that contains the elements from the diagonal 
+              of m in the order R(0)==m(0,0), R(1)==m(1,1), R(2)==m(2,2) and so on.
+    !*/
+
+    template <typename EXP>
+    struct diag_exp
+    {
+        /*!
+            WHAT THIS OBJECT REPRESENTS
+                This struct allows you to determine the type of matrix expression 
+                object returned from the diag() function.  An example makes its
+                use clear:
+
+                template <typename EXP>
+                void do_something( const matrix_exp<EXP>& mat)
+                {
+                    // d is a matrix expression that aliases mat.
+                    typename diag_exp<EXP>::type d = diag(mat);
+
+                    // Print the diagonal of mat.  So we see that by using
+                    // diag_exp we can save the object returned by diag() in
+                    // a local variable.    
+                    cout << d << endl;
+
+                    // Note that you can only save the return value of diag() to
+                    // a local variable if the argument to diag() has a lifetime
+                    // beyond the diag() expression.  The example shown above is
+                    // OK but the following would result in undefined behavior:
+                    typename diag_exp<EXP>::type bad = diag(mat + mat);
+                }
+        !*/
+        typedef type_of_expression_returned_by_diag type;
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp diagm (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_vector(m) == true
+              (i.e. m is a row or column matrix)
+        ensures
+            - returns a square matrix M such that:
+                - diag(M) == m
+                - non diagonal elements of M are 0
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp trans (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns the transpose of the matrix m
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_type::type dot (
+        const matrix_exp& m1,
+        const matrix_exp& m2
+    );
+    /*!
+        requires
+            - is_vector(m1) == true
+            - is_vector(m2) == true
+            - m1.size() == m2.size()
+            - m1.size() > 0
+        ensures
+            - returns the dot product between m1 and m2. That is, this function 
+              computes and returns the sum, for all i, of m1(i)*m2(i).
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp lowerm (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns a matrix M such that:
+                - M::type == the same type that was in m
+                - M has the same dimensions as m
+                - M is the lower triangular part of m.  That is:
+                    - if (r >= c) then
+                        - M(r,c) == m(r,c)
+                    - else
+                        - M(r,c) == 0
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp lowerm (
+        const matrix_exp& m,
+        const matrix_exp::type scalar_value
+    );
+    /*!
+        ensures
+            - returns a matrix M such that:
+                - M::type == the same type that was in m
+                - M has the same dimensions as m
+                - M is the lower triangular part of m except that the diagonal has
+                  been set to scalar_value.  That is:
+                    - if (r > c) then
+                        - M(r,c) == m(r,c)
+                    - else if (r == c) then
+                        - M(r,c) == scalar_value 
+                    - else
+                        - M(r,c) == 0
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp upperm (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns a matrix M such that:
+                - M::type == the same type that was in m
+                - M has the same dimensions as m
+                - M is the upper triangular part of m.  That is:
+                    - if (r <= c) then
+                        - M(r,c) == m(r,c)
+                    - else
+                        - M(r,c) == 0
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp upperm (
+        const matrix_exp& m,
+        const matrix_exp::type scalar_value
+    );
+    /*!
+        ensures
+            - returns a matrix M such that:
+                - M::type == the same type that was in m
+                - M has the same dimensions as m
+                - M is the upper triangular part of m except that the diagonal has
+                  been set to scalar_value.  That is:
+                    - if (r < c) then
+                        - M(r,c) == m(r,c)
+                    - else if (r == c) then
+                        - M(r,c) == scalar_value 
+                    - else
+                        - M(r,c) == 0
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp make_symmetric (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.nr() == m.nc()
+              (i.e. m must be a square matrix)
+        ensures
+            - returns a matrix M such that:
+                - M::type == the same type that was in m
+                - M has the same dimensions as m
+                - M is a symmetric matrix, that is, M == trans(M) and
+                  it is constructed from the lower triangular part of m.  Specifically,
+                  we have:
+                    - lowerm(M) == lowerm(m)
+                    - upperm(M) == trans(lowerm(m))
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T, 
+        long NR, 
+        long NC, 
+        T val
+        >
+    const matrix_exp uniform_matrix (
+    );
+    /*!
+        requires
+            - NR > 0 && NC > 0
+        ensures
+            - returns an NR by NC matrix with elements of type T and all set to val.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T,
+        long NR, 
+        long NC
+        >
+    const matrix_exp uniform_matrix (
+        const T& val
+    );
+    /*!
+        requires
+            - NR > 0 && NC > 0
+        ensures
+            - returns an NR by NC matrix with elements of type T and all set to val.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T
+        >
+    const matrix_exp uniform_matrix (
+        long nr,
+        long nc,
+        const T& val
+    );
+    /*!
+        requires
+            - nr >= 0 && nc >= 0
+        ensures
+            - returns an nr by nc matrix with elements of type T and all set to val.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp ones_matrix (
+        const matrix_exp& mat
+    );
+    /*!
+        requires
+            - mat.nr() >= 0 && mat.nc() >= 0
+        ensures
+            - Let T denote the type of element in mat. Then this function
+              returns uniform_matrix<T>(mat.nr(), mat.nc(), 1)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T
+        >
+    const matrix_exp ones_matrix (
+        long nr,
+        long nc
+    );
+    /*!
+        requires
+            - nr >= 0 && nc >= 0
+        ensures
+            - returns uniform_matrix<T>(nr, nc, 1)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp zeros_matrix (
+        const matrix_exp& mat
+    );
+    /*!
+        requires
+            - mat.nr() >= 0 && mat.nc() >= 0
+        ensures
+            - Let T denote the type of element in mat. Then this function
+              returns uniform_matrix<T>(mat.nr(), mat.nc(), 0)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T
+        >
+    const matrix_exp zeros_matrix (
+        long nr,
+        long nc
+    );
+    /*!
+        requires
+            - nr >= 0 && nc >= 0
+        ensures
+            - returns uniform_matrix<T>(nr, nc, 0)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp identity_matrix (
+        const matrix_exp& mat
+    );
+    /*!
+        requires
+            - mat.nr() == mat.nc()
+        ensures
+            - returns an identity matrix with the same dimensions as mat and
+              containing the same type of elements as mat.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T
+        >
+    const matrix_exp identity_matrix (
+        long N
+    );
+    /*!
+        requires
+            - N > 0
+        ensures
+            - returns an N by N identity matrix with elements of type T.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T, 
+        long N
+        >
+    const matrix_exp identity_matrix (
+    );
+    /*!
+        requires
+            - N > 0
+        ensures
+            - returns an N by N identity matrix with elements of type T.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp linspace (
+        double start,
+        double end,
+        long num
+    );
+    /*!
+        requires
+            - num >= 0
+        ensures
+            - returns a matrix M such that:
+                - M::type == double 
+                - is_row_vector(M) == true
+                - M.size() == num
+                - M == a row vector with num linearly spaced values beginning with start
+                  and stopping with end.  
+                - M(num-1) == end 
+                - if (num > 1) then
+                    - M(0) == start
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp logspace (
+        double start,
+        double end,
+        long num
+    );
+    /*!
+        requires
+            - num >= 0
+        ensures
+            - returns a matrix M such that:
+                - M::type == double 
+                - is_row_vector(M) == true
+                - M.size() == num
+                - M == a row vector with num logarithmically spaced values beginning with 
+                  10^start and stopping with 10^end.  
+                  (i.e. M == pow(10, linspace(start, end, num)))
+                - M(num-1) == 10^end
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp linpiece (
+        const double val,
+        const matrix_exp& joints
+    );
+    /*!
+        requires
+            - is_vector(joints) == true
+            - joints.size() >= 2
+            - for all valid i < j:
+                - joints(i) < joints(j)
+        ensures
+            - linpiece() is useful for creating piecewise linear functions of val.  For
+              example, if w is a parameter vector then you can represent a piecewise linear
+              function of val as: f(val) = dot(w, linpiece(val, linspace(0,100,5))).  In
+              this case, f(val) is piecewise linear on the intervals [0,25], [25,50],
+              [50,75], [75,100].  Moreover, w(i) defines the derivative of f(val) in the
+              i-th interval.  Finally, outside the interval [0,100] f(val) has a derivative
+              of zero and f(0) == 0.
+            - To be precise, this function returns a column vector L such that:
+                - L.size() == joints.size()-1
+                - is_col_vector(L) == true
+                - L contains the same type of elements as joints.
+                - for all valid i:
+                - if (joints(i) < val)
+                    - L(i) == min(val,joints(i+1)) - joints(i)
+                - else
+                    - L(i) == 0
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        long R,
+        long C
+        >
+    const matrix_exp rotate (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m
+                - R has the same dimensions as m
+                - for all valid r and c:
+                  R( (r+R)%m.nr() , (c+C)%m.nc() ) == m(r,c)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp fliplr (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - flips the matrix m from left to right and returns the result.  
+              I.e. reverses the order of the columns.
+            - returns a matrix M such that:
+                - M::type == the same type that was in m
+                - M has the same dimensions as m
+                - for all valid r and c:
+                  M(r,c) == m(r, m.nc()-c-1)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp flipud (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - flips the matrix m from up to down and returns the result.  
+              I.e. reverses the order of the rows.
+            - returns a matrix M such that:
+                - M::type == the same type that was in m
+                - M has the same dimensions as m
+                - for all valid r and c:
+                  M(r,c) == m(m.nr()-r-1, c)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp flip (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - flips the matrix m from up to down and left to right and returns the 
+              result.  I.e. returns flipud(fliplr(m)).
+            - returns a matrix M such that:
+                - M::type == the same type that was in m
+                - M has the same dimensions as m
+                - for all valid r and c:
+                  M(r,c) == m(m.nr()-r-1, m.nc()-c-1)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp reshape (
+        const matrix_exp& m,
+        long rows,
+        long cols
+    );
+    /*!
+        requires
+            - m.size() == rows*cols
+            - rows > 0
+            - cols > 0
+        ensures
+            - returns a matrix M such that: 
+                - M.nr() == rows
+                - M.nc() == cols
+                - M.size() == m.size()
+                - for all valid r and c:
+                    - let IDX = r*cols + c
+                    - M(r,c) == m(IDX/m.nc(), IDX%m.nc())
+
+            - i.e. The matrix m is reshaped into a new matrix of rows by cols
+              dimension.  Additionally, the elements of m are laid into M in row major 
+              order.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp reshape_to_column_vector (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns a matrix M such that: 
+                - is_col_vector(M) == true
+                - M.size() == m.size()
+                - for all valid r and c:
+                    - m(r,c) == M(r*m.nc() + c)
+
+            - i.e. The matrix m is reshaped into a column vector.  Note that
+              the elements are pulled out in row major order.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        long R,
+        long C
+        >
+    const matrix_exp removerc (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.nr() > R >= 0
+            - m.nc() > C >= 0
+        ensures
+            - returns a matrix M such that:
+                - M.nr() == m.nr() - 1
+                - M.nc() == m.nc() - 1
+                - M == m with its R row and C column removed
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp removerc (
+        const matrix_exp& m,
+        long R,
+        long C
+    );
+    /*!
+        requires
+            - m.nr() > R >= 0
+            - m.nc() > C >= 0
+        ensures
+            - returns a matrix M such that:
+                - M.nr() == m.nr() - 1
+                - M.nc() == m.nc() - 1
+                - M == m with its R row and C column removed
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        long R
+        >
+    const matrix_exp remove_row (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.nr() > R >= 0
+        ensures
+            - returns a matrix M such that:
+                - M.nr() == m.nr() - 1
+                - M.nc() == m.nc() 
+                - M == m with its R row removed
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp remove_row (
+        const matrix_exp& m,
+        long R
+    );
+    /*!
+        requires
+            - m.nr() > R >= 0
+        ensures
+            - returns a matrix M such that:
+                - M.nr() == m.nr() - 1
+                - M.nc() == m.nc() 
+                - M == m with its R row removed
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        long C
+        >
+    const matrix_exp remove_col (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.nc() > C >= 0
+        ensures
+            - returns a matrix M such that:
+                - M.nr() == m.nr() 
+                - M.nc() == m.nc() - 1 
+                - M == m with its C column removed
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp remove_col (
+        const matrix_exp& m,
+        long C
+    );
+    /*!
+        requires
+            - m.nc() > C >= 0
+        ensures
+            - returns a matrix M such that:
+                - M.nr() == m.nr() 
+                - M.nc() == m.nc() - 1 
+                - M == m with its C column removed
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+       typename target_type
+       >
+    const matrix_exp matrix_cast (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns a matrix R where for all valid r and c:
+              R(r,c) == static_cast<target_type>(m(r,c))
+              also, R has the same dimensions as m.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T,
+        long NR,
+        long NC,
+        typename MM,
+        typename U,
+        typename L
+        >
+    void set_all_elements (
+        matrix<T,NR,NC,MM,L>& m,
+        U value
+    );
+    /*!
+        ensures
+            - for all valid r and c:
+              m(r,c) == value
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::matrix_type tmp (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns a temporary matrix object that is a copy of m. 
+              (This allows you to easily force a matrix_exp to fully evaluate)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T, 
+        long NR, 
+        long NC, 
+        typename MM, 
+        typename L
+        >
+    uint32 hash (
+        const matrix<T,NR,NC,MM,L>& item,
+        uint32 seed = 0
+    );
+    /*!
+        requires
+            - T is a standard layout type (e.g. a POD type like int, float, 
+              or a simple struct).
+        ensures
+            - returns a 32bit hash of the data stored in item.  
+            - Each value of seed results in a different hash function being used.  
+              (e.g. hash(item,0) should generally not be equal to hash(item,1))
+            - uses the murmur_hash3() routine to compute the actual hash.
+            - Note that if the memory layout of the elements in item change between
+              hardware platforms then hash() will give different outputs.  If you want
+              hash() to always give the same output for the same input then you must 
+              ensure that elements of item always have the same layout in memory.
+              Typically this means using fixed width types and performing byte swapping
+              to account for endianness before passing item to hash().
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    // if matrix_exp contains non-complex types (e.g. float, double)
+    bool equal (
+        const matrix_exp& a,
+        const matrix_exp& b,
+        const matrix_exp::type epsilon = 100*std::numeric_limits<matrix_exp::type>::epsilon()
+    );
+    /*!
+        ensures
+            - if (a and b don't have the same dimensions) then
+                - returns false
+            - else if (there exists an r and c such that abs(a(r,c)-b(r,c)) > epsilon) then
+                - returns false
+            - else
+                - returns true
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    // if matrix_exp contains std::complex types 
+    bool equal (
+        const matrix_exp& a,
+        const matrix_exp& b,
+        const matrix_exp::type::value_type epsilon = 100*std::numeric_limits<matrix_exp::type::value_type>::epsilon()
+    );
+    /*!
+        ensures
+            - if (a and b don't have the same dimensions) then
+                - returns false
+            - else if (there exists an r and c such that abs(real(a(r,c)-b(r,c))) > epsilon 
+              or abs(imag(a(r,c)-b(r,c))) > epsilon) then
+                - returns false
+            - else
+                - returns true
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp pointwise_multiply (
+        const matrix_exp& a,
+        const matrix_exp& b 
+    );
+    /*!
+        requires
+            - a.nr() == b.nr()
+            - a.nc() == b.nc()
+            - a and b both contain the same type of element (one or both
+              can also be of type std::complex so long as the underlying type
+              in them is the same)
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in a and b.
+                - R has the same dimensions as a and b. 
+                - for all valid r and c:
+                  R(r,c) == a(r,c) * b(r,c)
+    !*/
+
+    const matrix_exp pointwise_multiply (
+        const matrix_exp& a,
+        const matrix_exp& b,
+        const matrix_exp& c 
+    );
+    /*!
+        performs pointwise_multiply(a,pointwise_multiply(b,c));
+    !*/
+
+    const matrix_exp pointwise_multiply (
+        const matrix_exp& a,
+        const matrix_exp& b,
+        const matrix_exp& c,
+        const matrix_exp& d 
+    );
+    /*!
+        performs pointwise_multiply(pointwise_multiply(a,b),pointwise_multiply(c,d));
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp join_rows (
+        const matrix_exp& a,
+        const matrix_exp& b 
+    );
+    /*!
+        requires
+            - a.nr() == b.nr() || a.size() == 0 || b.size() == 0
+            - a and b both contain the same type of element
+        ensures
+            - This function joins two matrices together by concatenating their rows.
+            - returns a matrix R such that:
+                - R::type == the same type that was in a and b.
+                - R.nr() == a.nr() == b.nr()
+                - R.nc() == a.nc() + b.nc()
+                - for all valid r and c:
+                    - if (c < a.nc()) then
+                        - R(r,c) == a(r,c) 
+                    - else
+                        - R(r,c) == b(r, c-a.nc()) 
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp join_cols (
+        const matrix_exp& a,
+        const matrix_exp& b 
+    );
+    /*!
+        requires
+            - a.nc() == b.nc() || a.size() == 0 || b.size() == 0
+            - a and b both contain the same type of element
+        ensures
+            - This function joins two matrices together by concatenating their columns.
+            - returns a matrix R such that:
+                - R::type == the same type that was in a and b.
+                - R.nr() == a.nr() + b.nr()
+                - R.nc() == a.nc() == b.nc()
+                - for all valid r and c:
+                    - if (r < a.nr()) then
+                        - R(r,c) == a(r,c) 
+                    - else
+                        - R(r,c) == b(r-a.nr(), c) 
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp tensor_product (
+        const matrix_exp& a,
+        const matrix_exp& b 
+    );
+    /*!
+        requires
+            - a and b both contain the same type of element
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in a and b.
+                - R.nr() == a.nr() * b.nr()  
+                - R.nc() == a.nc() * b.nc()  
+                - for all valid r and c:
+                  R(r,c) == a(r/b.nr(), c/b.nc()) * b(r%b.nr(), c%b.nc())
+                - I.e. R is the tensor product of matrix a with matrix b
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp cartesian_product (
+        const matrix_exp& A,
+        const matrix_exp& B 
+    );
+    /*!
+        requires
+            - A and B both contain the same type of element
+        ensures
+            - Think of A and B as sets of column vectors.  Then this function 
+              returns a matrix that contains a set of column vectors that is
+              the Cartesian product of the sets A and B.  That is, the resulting
+              matrix contains every possible combination of vectors from both A and
+              B.
+            - returns a matrix R such that:
+                - R::type == the same type that was in A and B.
+                - R.nr() == A.nr() + B.nr()  
+                - R.nc() == A.nc() * B.nc()  
+                - Each column of R is the concatenation of a column vector
+                  from A with a column vector from B.  
+                - for all valid r and c:
+                    - if (r < A.nr()) then
+                        - R(r,c) == A(r, c/B.nc())
+                    - else
+                        - R(r,c) == B(r-A.nr(), c%B.nc())
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp scale_columns (
+        const matrix_exp& m,
+        const matrix_exp& v
+    );
+    /*!
+        requires
+            - is_vector(v) == true
+            - v.size() == m.nc()
+            - m and v both contain the same type of element
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m and v.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                  R(r,c) == m(r,c) * v(c)
+                - i.e. R is the result of multiplying each of m's columns by
+                  the corresponding scalar in v.
+
+            - Note that this function is identical to the expression m*diagm(v).  
+              That is, the * operator is overloaded for this case and will invoke
+              scale_columns() automatically as appropriate.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp scale_rows (
+        const matrix_exp& m,
+        const matrix_exp& v
+    );
+    /*!
+        requires
+            - is_vector(v) == true
+            - v.size() == m.nr()
+            - m and v both contain the same type of element
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m and v.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                  R(r,c) == m(r,c) * v(r)
+                - i.e. R is the result of multiplying each of m's rows by
+                  the corresponding scalar in v.
+
+            - Note that this function is identical to the expression diagm(v)*m.  
+              That is, the * operator is overloaded for this case and will invoke
+              scale_rows() automatically as appropriate.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename T>
+    void sort_columns (
+        matrix<T>& m,
+        matrix<T>& v
+    );
+    /*!
+        requires
+            - is_col_vector(v) == true
+            - v.size() == m.nc()
+            - m and v both contain the same type of element
+        ensures
+            - the dimensions for m and v are not changed
+            - sorts the columns of m according to the values in v.
+              i.e. 
+                - #v == the contents of v but in sorted order according to
+                  operator<.  So smaller elements come first.
+                - Let #v(new(i)) == v(i) (i.e. new(i) is the index element i moved to)
+                - colm(#m,new(i)) == colm(m,i) 
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename T>
+    void rsort_columns (
+        matrix<T>& m,
+        matrix<T>& v
+    );
+    /*!
+        requires
+            - is_col_vector(v) == true
+            - v.size() == m.nc()
+            - m and v both contain the same type of element
+        ensures
+            - the dimensions for m and v are not changed
+            - sorts the columns of m according to the values in v.
+              i.e. 
+                - #v == the contents of v but in sorted order according to
+                  operator>.  So larger elements come first.
+                - Let #v(new(i)) == v(i) (i.e. new(i) is the index element i moved to)
+                - colm(#m,new(i)) == colm(m,i) 
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::type length_squared (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_vector(m) == true
+        ensures
+            - returns sum(squared(m))
+              (i.e. returns the square of the length of the vector m)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::type length (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_vector(m) == true
+        ensures
+            - returns sqrt(sum(squared(m)))
+              (i.e. returns the length of the vector m)
+            - if (m contains integer valued elements) then  
+                - The return type is a double that represents the length.  Therefore, the
+                  return value of length() is always represented using a floating point
+                  type. 
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    bool is_row_vector (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - if (m.nr() == 1) then
+                - return true
+            - else
+                - returns false
+    !*/
+
+    bool is_col_vector (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - if (m.nc() == 1) then
+                - return true
+            - else
+                - returns false
+    !*/
+
+    bool is_vector (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - if (is_row_vector(m) || is_col_vector(m)) then
+                - return true
+            - else
+                - returns false
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    bool is_finite (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns true if all the values in m are finite values and also not any kind
+              of NaN value.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+//                      Thresholding relational operators 
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator< (
+        const matrix_exp& m,
+        const S& s
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (m(r,c) < s) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator< (
+        const S& s,
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (s < m(r,c)) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator<= (
+        const matrix_exp& m,
+        const S& s
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (m(r,c) <= s) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator<= (
+        const S& s,
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (s <= m(r,c)) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator> (
+        const matrix_exp& m,
+        const S& s
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (m(r,c) > s) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator> (
+        const S& s,
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (s > m(r,c)) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator>= (
+        const matrix_exp& m,
+        const S& s
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (m(r,c) >= s) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator>= (
+        const S& s,
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (s >= m(r,c)) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator== (
+        const matrix_exp& m,
+        const S& s
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (m(r,c) == s) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator== (
+        const S& s,
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (s == m(r,c)) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator!= (
+        const matrix_exp& m,
+        const S& s
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (m(r,c) != s) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename S>
+    const matrix_exp operator!= (
+        const S& s,
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_built_in_scalar_type<S>::value == true 
+            - is_built_in_scalar_type<matrix_exp::type>::value == true
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m.
+                - R has the same dimensions as m. 
+                - for all valid r and c:
+                    - if (s != m(r,c)) then
+                        - R(r,c) == 1
+                    - else
+                        - R(r,c) == 0
+                - i.e. R is a binary matrix of all 1s or 0s.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+//                              Statistics
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::type min (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.size() > 0
+        ensures
+            - returns the value of the smallest element of m.  If m contains complex
+              elements then the element returned is the one with the smallest norm
+              according to std::norm().
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp min_pointwise (
+        const matrix_exp& a,
+        const matrix_exp& b 
+    );
+    /*!
+        requires
+            - a.nr() == b.nr()
+            - a.nc() == b.nc()
+            - a and b both contain the same type of element 
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in a and b.
+                - R has the same dimensions as a and b. 
+                - for all valid r and c:
+                  R(r,c) == std::min(a(r,c), b(r,c))
+    !*/
+
+    const matrix_exp min_pointwise (
+        const matrix_exp& a,
+        const matrix_exp& b,
+        const matrix_exp& c 
+    );
+    /*!
+        performs min_pointwise(a,min_pointwise(b,c));
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::type max (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.size() > 0
+        ensures
+            - returns the value of the biggest element of m.  If m contains complex
+              elements then the element returned is the one with the largest norm
+              according to std::norm().
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp max_pointwise (
+        const matrix_exp& a,
+        const matrix_exp& b 
+    );
+    /*!
+        requires
+            - a.nr() == b.nr()
+            - a.nc() == b.nc()
+            - a and b both contain the same type of element 
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in a and b.
+                - R has the same dimensions as a and b. 
+                - for all valid r and c:
+                  R(r,c) == std::max(a(r,c), b(r,c))
+    !*/
+
+    const matrix_exp max_pointwise (
+        const matrix_exp& a,
+        const matrix_exp& b,
+        const matrix_exp& c 
+    );
+    /*!
+        performs max_pointwise(a,max_pointwise(b,c));
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    void find_min_and_max (
+        const matrix_exp& m,
+        matrix_exp::type& min_val,
+        matrix_exp::type& max_val
+    );
+    /*!
+        requires
+            - m.size() > 0
+        ensures
+            - #min_val == min(m)
+            - #max_val == max(m)
+            - This function computes both the min and max in just one pass
+              over the elements of the matrix m.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    long index_of_max (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_vector(m) == true
+            - m.size() > 0 
+        ensures
+            - returns the index of the largest element in m.  
+              (i.e. m(index_of_max(m)) == max(m))
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    long index_of_min (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - is_vector(m) == true
+            - m.size() > 0 
+        ensures
+            - returns the index of the smallest element in m.  
+              (i.e. m(index_of_min(m)) == min(m))
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    point max_point (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.size() > 0
+        ensures
+            - returns the location of the maximum element of the array, that is, if the
+              returned point is P then it will be the case that: m(P.y(),P.x()) == max(m).
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    dlib::vector<double,2> max_point_interpolated (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.size() > 0
+        ensures
+            - Like max_point(), this function finds the location in m with the largest
+              value.  However, we additionally use some quadratic interpolation to find the
+              location of the maximum point with sub-pixel accuracy.  Therefore, the
+              returned point is equal to max_point(m) + some small sub-pixel delta.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    point min_point (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.size() > 0
+        ensures
+            - returns the location of the minimum element of the array, that is, if the
+              returned point is P then it will be the case that: m(P.y(),P.x()) == min(m).
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::type sum (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns the sum of all elements in m
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp sum_rows (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.size() > 0
+        ensures
+            - returns a row matrix that contains the sum of all the rows in m. 
+            - returns a matrix M such that
+                - M::type == the same type that was in m
+                - M.nr() == 1
+                - M.nc() == m.nc()
+                - for all valid i:
+                    - M(i) == sum(colm(m,i)) 
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp sum_cols (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - m.size() > 0
+        ensures
+            - returns a column matrix that contains the sum of all the columns in m. 
+            - returns a matrix M such that
+                - M::type == the same type that was in m
+                - M.nr() == m.nr() 
+                - M.nc() == 1
+                - for all valid i:
+                    - M(i) == sum(rowm(m,i)) 
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::type prod (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns the results of multiplying all elements of m together. 
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::type mean (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns the mean of all elements in m. 
+              (i.e. returns sum(m)/(m.nr()*m.nc()))
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::type variance (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns the unbiased sample variance of all elements in m 
+              (i.e. 1.0/(m.nr()*m.nc() - 1)*(sum of all pow(m(i,j) - mean(m),2)))
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp::type stddev (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns sqrt(variance(m))
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix covariance (
+        const matrix_exp& m
+    );
+    /*!
+        requires
+            - matrix_exp::type == a dlib::matrix object
+            - is_col_vector(m) == true
+            - m.size() > 1
+            - for all valid i, j:
+                - is_col_vector(m(i)) == true 
+                - m(i).size() > 0
+                - m(i).size() == m(j).size() 
+                - i.e. m contains only column vectors and all the column vectors
+                  have the same non-zero length
+        ensures
+            - returns the unbiased sample covariance matrix for the set of samples
+              in m.  
+              (i.e. 1.0/(m.nr()-1)*(sum of all (m(i) - mean(m))*trans(m(i) - mean(m))))
+            - the returned matrix will contain elements of type matrix_exp::type::type.
+            - the returned matrix will have m(0).nr() rows and columns.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename rand_gen>
+    const matrix<double> randm( 
+        long nr,
+        long nc,
+        rand_gen& rnd
+    );
+    /*!
+        requires
+            - nr >= 0
+            - nc >= 0
+            - rand_gen == an object that implements the rand/rand_float_abstract.h interface
+        ensures
+            - generates a random matrix using the given rnd random number generator
+            - returns a matrix M such that
+                - M::type == double
+                - M.nr() == nr
+                - M.nc() == nc
+                - for all valid i, j:
+                    - M(i,j) == a random number such that 0 <= M(i,j) < 1
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    inline const matrix<double> randm( 
+        long nr,
+        long nc
+    );
+    /*!
+        requires
+            - nr >= 0
+            - nc >= 0
+        ensures
+            - generates a random matrix using std::rand() 
+            - returns a matrix M such that
+                - M::type == double
+                - M.nr() == nr
+                - M.nc() == nc
+                - for all valid i, j:
+                    - M(i,j) == a random number such that 0 <= M(i,j) < 1
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    inline const matrix_exp gaussian_randm (
+        long nr,
+        long nc,
+        unsigned long seed = 0
+    );
+    /*!
+        requires
+            - nr >= 0
+            - nc >= 0
+        ensures
+            - returns a matrix with its values filled with 0 mean unit variance Gaussian
+              random numbers.  
+            - Each setting of the seed results in a different random matrix.
+            - The returned matrix is lazily evaluated using the expression templates
+              technique.  This means that the returned matrix doesn't take up any memory
+              and is only an expression template.  The values themselves are computed on
+              demand using the gaussian_random_hash() routine.  
+            - returns a matrix M such that
+                - M::type == double
+                - M.nr() == nr
+                - M.nc() == nc
+                - for all valid i, j:
+                    - M(i,j) == gaussian_random_hash(i,j,seed) 
+    !*/
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+//                                 Pixel and Image Utilities
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename T,
+        typename P
+        >
+    const matrix<T,pixel_traits<P>::num,1> pixel_to_vector (
+        const P& pixel
+    );
+    /*!
+        requires
+            - pixel_traits<P> must be defined
+        ensures
+            - returns a matrix M such that:
+                - M::type == T
+                - M::NC == 1 
+                - M::NR == pixel_traits<P>::num
+                - if (pixel_traits<P>::grayscale) then
+                    - M(0) == pixel 
+                - if (pixel_traits<P>::rgb) then
+                    - M(0) == pixel.red 
+                    - M(1) == pixel.green 
+                    - M(2) == pixel.blue 
+                - if (pixel_traits<P>::hsi) then
+                    - M(0) == pixel.h 
+                    - M(1) == pixel.s 
+                    - M(2) == pixel.i 
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename P
+        >
+    void vector_to_pixel (
+        P& pixel,
+        const matrix_exp& vector 
+    );
+    /*!
+        requires
+            - vector::NR == pixel_traits<P>::num
+            - vector::NC == 1 
+              (i.e. you have to use a statically dimensioned vector)
+        ensures
+            - if (pixel_traits<P>::grayscale) then
+                - pixel == M(0) 
+            - if (pixel_traits<P>::rgb) then
+                - pixel.red   == M(0)  
+                - pixel.green == M(1) 
+                - pixel.blue  == M(2)  
+            - if (pixel_traits<P>::hsi) then
+                - pixel.h == M(0)
+                - pixel.s == M(1)
+                - pixel.i == M(2)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        long lower,
+        long upper 
+        >
+    const matrix_exp clamp (
+        const matrix_exp& m
+    );
+    /*!
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m
+                - R has the same dimensions as m
+                - for all valid r and c:
+                    - if (m(r,c) > upper) then
+                        - R(r,c) == upper
+                    - else if (m(r,c) < lower) then
+                        - R(r,c) == lower
+                    - else
+                        - R(r,c) == m(r,c)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp clamp (
+        const matrix_exp& m,
+        const matrix_exp::type& lower,
+        const matrix_exp::type& upper
+    );
+    /*!
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m
+                - R has the same dimensions as m
+                - for all valid r and c:
+                    - if (m(r,c) > upper) then
+                        - R(r,c) == upper
+                    - else if (m(r,c) < lower) then
+                        - R(r,c) == lower
+                    - else
+                        - R(r,c) == m(r,c)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp clamp (
+        const matrix_exp& m,
+        const matrix_exp& lower,
+        const matrix_exp& upper
+    );
+    /*!
+        requires
+            - m.nr() == lower.nr()
+            - m.nc() == lower.nc()
+            - m.nr() == upper.nr()
+            - m.nc() == upper.nc()
+            - m, lower, and upper all contain the same type of elements. 
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m
+                - R has the same dimensions as m
+                - for all valid r and c:
+                    - if (m(r,c) > upper(r,c)) then
+                        - R(r,c) == upper(r,c)
+                    - else if (m(r,c) < lower(r,c)) then
+                        - R(r,c) == lower(r,c)
+                    - else
+                        - R(r,c) == m(r,c)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp lowerbound (
+        const matrix_exp& m,
+        const matrix_exp::type& thresh 
+    );
+    /*!
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m
+                - R has the same dimensions as m
+                - for all valid r and c:
+                    - if (m(r,c) >= thresh) then
+                        - R(r,c) == m(r,c)
+                    - else
+                        - R(r,c) == thresh
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    const matrix_exp upperbound (
+        const matrix_exp& m,
+        const matrix_exp::type& thresh 
+    );
+    /*!
+        ensures
+            - returns a matrix R such that:
+                - R::type == the same type that was in m
+                - R has the same dimensions as m
+                - for all valid r and c:
+                    - if (m(r,c) <= thresh) then
+                        - R(r,c) == m(r,c)
+                    - else
+                        - R(r,c) == thresh
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_MATRIx_UTILITIES_ABSTRACT_
+