// Copyright (C) 2006 Davis E. King (davisking@users.sourceforge.net)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_MATRIx_UTILITIES_
#define DLIB_MATRIx_UTILITIES_
#include "matrix_utilities_abstract.h"
#include "matrix.h"
#include <cmath>
#include <complex>
#include <limits>
#include "../pixel.h"
#include "../geometry/rectangle.h"
#include "../stl_checked.h"
#include <vector>
namespace dlib
{
// ----------------------------------------------------------------------------------------
/*
templates for finding the max of two matrix expressions' dimensions
*/
template <typename EXP1, typename EXP2 = void, typename EXP3 = void, typename EXP4 = void>
struct max_nr;
template <typename EXP1>
struct max_nr<EXP1,void,void,void>
{
const static long val = EXP1::NR;
};
template <typename EXP1, typename EXP2>
struct max_nr<EXP1,EXP2,void,void>
{
const static long val = (EXP1::NR > EXP2::NR) ? (EXP1::NR) : (EXP2::NR);
};
template <typename EXP1, typename EXP2, typename EXP3>
struct max_nr<EXP1,EXP2,EXP3,void>
{
private:
const static long max12 = (EXP1::NR > EXP2::NR) ? (EXP1::NR) : (EXP2::NR);
public:
const static long val = (max12 > EXP3::NR) ? (max12) : (EXP3::NR);
};
template <typename EXP1, typename EXP2, typename EXP3, typename EXP4>
struct max_nr
{
private:
const static long max12 = (EXP1::NR > EXP2::NR) ? (EXP1::NR) : (EXP2::NR);
const static long max34 = (EXP3::NR > EXP4::NR) ? (EXP3::NR) : (EXP4::NR);
public:
const static long val = (max12 > max34) ? (max12) : (max34);
};
template <typename EXP1, typename EXP2 = void, typename EXP3 = void, typename EXP4 = void>
struct max_nc;
template <typename EXP1>
struct max_nc<EXP1,void,void,void>
{
const static long val = EXP1::NC;
};
template <typename EXP1, typename EXP2>
struct max_nc<EXP1,EXP2,void,void>
{
const static long val = (EXP1::NC > EXP2::NC) ? (EXP1::NC) : (EXP2::NC);
};
template <typename EXP1, typename EXP2, typename EXP3>
struct max_nc<EXP1,EXP2,EXP3,void>
{
private:
const static long max12 = (EXP1::NC > EXP2::NC) ? (EXP1::NC) : (EXP2::NC);
public:
const static long val = (max12 > EXP3::NC) ? (max12) : (EXP3::NC);
};
template <typename EXP1, typename EXP2, typename EXP3, typename EXP4>
struct max_nc
{
private:
const static long max12 = (EXP1::NC > EXP2::NC) ? (EXP1::NC) : (EXP2::NC);
const static long max34 = (EXP3::NC > EXP4::NC) ? (EXP3::NC) : (EXP4::NC);
public:
const static long val = (max12 > max34) ? (max12) : (max34);
};
// ----------------------------------------------------------------------------------------
template <
typename OP
>
class matrix_zeroary_exp;
template <
typename M,
typename OP
>
class matrix_unary_exp;
template <
typename M1,
typename M2,
typename OP
>
class matrix_binary_exp;
struct has_destructive_aliasing
{
template <typename M, typename U, long iNR, long iNC, typename MM >
static bool destructively_aliases (
const M& m,
const matrix<U,iNR,iNC,MM>& item
) { return m.aliases(item); }
template <typename M1, typename M2, typename U, long iNR, long iNC, typename MM >
static bool destructively_aliases (
const M1& m1,
const M2& m2,
const matrix<U,iNR,iNC,MM>& item
) { return m1.aliases(item) || m2.aliases(item) ; }
};
struct has_nondestructive_aliasing
{
template <typename M, typename U, long iNR, long iNC, typename MM >
static bool destructively_aliases (
const M& m,
const matrix<U,iNR,iNC,MM>& item
) { return m.destructively_aliases(item); }
template <typename M1, typename M2, typename U, long iNR, long iNC, typename MM >
static bool destructively_aliases (
const M1& m1,
const M2& m2,
const matrix<U,iNR,iNC, MM>& item
) { return m1.destructively_aliases(item) || m2.destructively_aliases(item) ; }
};
template <typename EXP1, typename EXP2 = void, typename EXP3 = void, typename EXP4 = void>
struct preserves_dimensions
{
const static long NR = max_nr<EXP1,EXP2,EXP3,EXP4>::val;
const static long NC = max_nc<EXP1,EXP2,EXP3,EXP4>::val;
typedef typename EXP1::mem_manager_type mem_manager_type;
template <typename M>
static long nr (const M& m) { return m.nr(); }
template <typename M>
static long nc (const M& m) { return m.nc(); }
template <typename M1, typename M2>
static long nr (const M1& m1, const M2& ) { return m1.nr(); }
template <typename M1, typename M2>
static long nc (const M1& m1, const M2& ) { return m1.nc(); }
};
// ----------------------------------------------------------------------------------------
template <
typename EXP
>
const typename matrix_exp<EXP>::type max (
const matrix_exp<EXP>& m
)
{
DLIB_ASSERT(m.size() > 0,
"\ttype max(const matrix_exp& m)"
<< "\n\tYou can't ask for the max() of an empty matrix"
<< "\n\tm.size(): " << m.size()
);
typedef typename matrix_exp<EXP>::type type;
type val = m(0,0);
for (long r = 0; r < m.nr(); ++r)
{
for (long c = 0; c < m.nc(); ++c)
{
type temp = m(r,c);
if (temp > val)
val = temp;
}
}
return val;
}
// ----------------------------------------------------------------------------------------
template <
typename EXP
>
const typename matrix_exp<EXP>::type min (
const matrix_exp<EXP>& m
)
{
DLIB_ASSERT(m.size() > 0,
"\ttype min(const matrix_exp& m)"
<< "\n\tYou can't ask for the min() of an empty matrix"
<< "\n\tm.size(): " << m.size()
);
typedef typename matrix_exp<EXP>::type type;
type val = m(0,0);
for (long r = 0; r < m.nr(); ++r)
{
for (long c = 0; c < m.nc(); ++c)
{
type temp = m(r,c);
if (temp < val)
val = temp;
}
}
return val;
}
// ----------------------------------------------------------------------------------------
template <
typename EXP
>
const typename matrix_exp<EXP>::type length (
const matrix_exp<EXP>& m
)
{
DLIB_ASSERT(m.nr() == 1 || m.nc() == 1,
"\ttype length(const matrix_exp& m)"
<< "\n\tm must be a row or column vector"
<< "\n\tm.nr(): " << m.nr()
<< "\n\tm.nc(): " << m.nc()
);
return std::sqrt(sum(squared(m)));
}
// ----------------------------------------------------------------------------------------
template <
typename EXP
>
const typename matrix_exp<EXP>::type length_squared (
const matrix_exp<EXP>& m
)
{
DLIB_ASSERT(m.nr() == 1 || m.nc() == 1,
"\ttype length_squared(const matrix_exp& m)"
<< "\n\tm must be a row or column vector"
<< "\n\tm.nr(): " << m.nr()
<< "\n\tm.nc(): " << m.nc()
);
return sum(squared(m));
}
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
namespace nric
{
// This namespace contains stuff from Numerical Recipes in C
template <typename T>
inline T pythag(const T& a, const T& b)
{
T absa,absb;
absa=std::abs(a);
absb=std::abs(b);
if (absa > absb)
{
T val = absb/absa;
val *= val;
return absa*std::sqrt(1.0+val);
}
else
{
if (absb == 0.0)
{
return 0.0;
}
else
{
T val = absa/absb;
val *= val;
return absb*std::sqrt(1.0+val);
}
}
}
template <typename T>
inline T sign(const T& a, const T& b)
{
if (b < 0)
{
return -std::abs(a);
}
else
{
return std::abs(a);
}
}
template <
typename T,
long M, long N,
long wN, long wX,
long vN,
long rN, long rX,
typename MM1,
typename MM2,
typename MM3,
typename MM4
>
bool svdcmp(
matrix<T,M,N,MM1>& a,
matrix<T,wN,wX,MM2>& w,
matrix<T,vN,vN,MM3>& v,
matrix<T,rN,rX,MM4>& rv1
)
/*! ( this function is derived from the one in numerical recipes in C chapter 2.6)
requires
- w.nr() == a.nc()
- w.nc() == 1
- v.nr() == a.nc()
- v.nc() == a.nc()
- rv1.nr() == a.nc()
- rv1.nc() == 1
ensures
- computes the singular value decomposition of a
- let W be the matrix such that diag(W) == #w then:
- a == #a*W*trans(#v)
- trans(#a)*#a == identity matrix
- trans(#v)*#v == identity matrix
- #rv1 == some undefined value
- returns true for success and false for failure
!*/
{
DLIB_ASSERT(
w.nr() == a.nc() &&
w.nc() == 1 &&
v.nr() == a.nc() &&
v.nc() == a.nc() &&
rv1.nr() == a.nc() &&
rv1.nc() == 1, "");
COMPILE_TIME_ASSERT(wX == 0 || wX == 1);
COMPILE_TIME_ASSERT(rX == 0 || rX == 1);
const T one = 1.0;
const long max_iter = 30;
const long n = a.nc();
const long m = a.nr();
const T eps = std::numeric_limits<T>::epsilon();
long nm = 0, l = 0;
bool flag;
T anorm,c,f,g,h,s,scale,x,y,z;
g = 0.0;
scale = 0.0;
anorm = 0.0;
for (long i = 0; i < n; ++i)
{
l = i+1;
rv1(i) = scale*g;
g = s = scale = 0.0;
if (i < m)
{
for (long k = i; k < m; ++k)
scale += std::abs(a(k,i));
if (scale)
{
for (long k = i; k < m; ++k)
{
a(k,i) /= scale;
s += a(k,i)*a(k,i);
}
f = a(i,i);
g = -sign(std::sqrt(s),f);
h = f*g - s;
a(i,i) = f - g;
for (long j = l; j < n; ++j)
{
s = 0.0;
for (long k = i; k < m; ++k)
s += a(k,i)*a(k,j);
f = s/h;
for (long k = i; k < m; ++k)
a(k,j) += f*a(k,i);
}
for (long k = i; k < m; ++k)
a(k,i) *= scale;
}
}
w(i) = scale *g;
g=s=scale=0.0;
if (i < m && i < n-1)
{
for (long k = l; k < n; ++k)
scale += std::abs(a(i,k));
if (scale)
{
for (long k = l; k < n; ++k)
{
a(i,k) /= scale;
s += a(i,k)*a(i,k);
}
f = a(i,l);
g = -sign(std::sqrt(s),f);
h = f*g - s;
a(i,l) = f - g;
for (long k = l; k < n; ++k)
rv1(k) = a(i,k)/h;
for (long j = l; j < m; ++j)
{
s = 0.0;
for (long k = l; k < n; ++k)
s += a(j,k)*a(i,k);
for (long k = l; k < n; ++k)
a(j,k) += s*rv1(k);
}
for (long k = l; k < n; ++k)
a(i,k) *= scale;
}
}
anorm = std::max(anorm,(std::abs(w(i))+std::abs(rv1(i))));
}
for (long i = n-1; i >= 0; --i)
{
if (i < n-1)
{
if (g != 0)
{
for (long j = l; j < n ; ++j)
v(j,i) = (a(i,j)/a(i,l))/g;
for (long j = l; j < n; ++j)
{
s = 0.0;
for (long k = l; k < n; ++k)
s += a(i,k)*v(k,j);
for (long k = l; k < n; ++k)
v(k,j) += s*v(k,i);
}
}
for (long j = l; j < n; ++j)
v(i,j) = v(j,i) = 0.0;
}
v(i,i) = 1.0;
g = rv1(i);
l = i;
}
for (long i = std::min(m,n)-1; i >= 0; --i)
{
l = i + 1;
g = w(i);
for (long j = l; j < n; ++j)
a(i,j) = 0.0;
if (g != 0)
{
g = 1.0/g;
for (long j = l; j < n; ++j)
{
s = 0.0;
for (long k = l; k < m; ++k)
s += a(k,i)*a(k,j);
f=(s/a(i,i))*g;
for (long k = i; k < m; ++k)
a(k,j) += f*a(k,i);
}
for (long j = i; j < m; ++j)
a(j,i) *= g;
}
else
{
for (long j = i; j < m; ++j)
a(j,i) = 0.0;
}
++a(i,i);
}
for (long k = n-1; k >= 0; --k)
{
for (long its = 1; its <= max_iter; ++its)
{
flag = true;
for (l = k; l >= 1; --l)
{
nm = l - 1;
if (std::abs(rv1(l)) <= eps*anorm)
{
flag = false;
break;
}
if (std::abs(w(nm)) <= eps*anorm)
{
break;
}
}
if (flag)
{
c = 0.0;
s = 1.0;
for (long i = l; i <= k; ++i)
{
f = s*rv1(i);
rv1(i) = c*rv1(i);
if (std::abs(f) <= eps*anorm)
break;
g = w(i);
h = pythag(f,g);
w(i) = h;
h = 1.0/h;
c = g*h;
s = -f*h;
for (long j = 0; j < m; ++j)
{
y = a(j,nm);
z = a(j,i);
a(j,nm) = y*c + z*s;
a(j,i) = z*c - y*s;
}
}
}
z = w(k);
if (l == k)
{
if (z < 0.0)
{
w(k) = -z;
for (long j = 0; j < n; ++j)
v(j,k) = -v(j,k);
}
break;
}
if (its == max_iter)
return false;
x = w(l);
nm = k - 1;
y = w(nm);
g = rv1(nm);
h = rv1(k);
f = ((y-z)*(y+z) + (g-h)*(g+h))/(2.0*h*y);
g = pythag(f,one);
f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x;
c = s = 1.0;
for (long j = l; j <= nm; ++j)
{
long i = j + 1;
g = rv1(i);
y = w(i);
h = s*g;
g = c*g;
z = pythag(f,h);
rv1(j) = z;
c = f/z;
s = h/z;
f = x*c + g*s;
g = g*c - x*s;
h = y*s;
y *= c;
for (long jj = 0; jj < n; ++jj)
{
x = v(jj,j);
z = v(jj,i);
v(jj,j) = x*c + z*s;
v(jj,i) = z*c - x*s;
}
z = pythag(f,h);
w(j) = z;
if (z != 0)
{
z = 1.0/z;
c = f*z;
s = h*z;
}
f = c*g + s*y;
x = c*y - s*g;
for (long jj = 0; jj < m; ++jj)
{
y = a(jj,j);
z = a(jj,i);
a(jj,j) = y*c + z*s;
a(jj,i) = z*c - y*s;
}
}
rv1(l) = 0.0;
rv1(k) = f;
w(k) = x;
}
}
return true;
}
template <
typename T,
long N,
long NX,
typename MM1,
typename MM2,
typename MM3
>
bool ludcmp (
matrix<T,N,N,MM1>& a,
matrix<long,N,NX,MM2>& indx,
T& d,
matrix<T,N,NX,MM3> vv
)
/*!
( this function is derived from the one in numerical recipes in C chapter 2.3)
ensures
- #a == both the L and U matrices
- #indx == the permutation vector (see numerical recipes in C)
- #d == some other thing (see numerical recipes in C)
- #vv == some undefined value. this is just used for scratch space
- if (the matrix is singular and we can't do anything) then
- returns false
- else
- returns true
!*/
{
DLIB_ASSERT(indx.nc() == 1,"in dlib::nric::ludcmp() the indx matrix must be a column vector");
DLIB_ASSERT(vv.nc() == 1,"in dlib::nric::ludcmp() the vv matrix must be a column vector");
const long n = a.nr();
long imax = 0;
T big, dum, sum, temp;
d = 1.0;
for (long i = 0; i < n; ++i)
{
big = 0;
for (long j = 0; j < n; ++j)
{
if ((temp=std::abs(a(i,j))) > big)
big = temp;
}
if (big == 0.0)
{
return false;
}
vv(i) = 1/big;
}
for (long j = 0; j < n; ++j)
{
for (long i = 0; i < j; ++i)
{
sum = a(i,j);
for (long k = 0; k < i; ++k)
sum -= a(i,k)*a(k,j);
a(i,j) = sum;
}
big = 0;
for (long i = j; i < n; ++i)
{
sum = a(i,j);
for (long k = 0; k < j; ++k)
sum -= a(i,k)*a(k,j);
a(i,j) = sum;
if ( (dum=vv(i)*std::abs(sum)) >= big)
{
big = dum;
imax = i;
}
}
if (j != imax)
{
for (long k = 0; k < n; ++k)
{
dum = a(imax,k);
a(imax,k) = a(j,k);
a(j,k) = dum;
}
d = -d;
vv(imax) = vv(j);
}
indx(j) = imax;
if (j < n-1)
{
dum = 1/a(j,j);
for (long i = j+1; i < n; ++i)
a(i,j) *= dum;
}
}
return true;
}
// ----------------------------------------------------------------------------------------
template <
typename T,
long N,
long NX,
typename MM1,
typename MM2,
typename MM3
>
void lubksb (
const matrix<T,N,N,MM1>& a,
const matrix<long,N,NX,MM2>& indx,
matrix<T,N,NX,MM3>& b
)
/*!
( this function is derived from the one in numerical recipes in C chapter 2.3)
requires
- a == the LU decomposition you get from ludcmp()
- indx == the indx term you get out of ludcmp()
- b == the right hand side vector from the expression a*x = b
ensures
- #b == the solution vector x from the expression a*x = b
(basically, this function solves for x given b and a)
!*/
{
DLIB_ASSERT(indx.nc() == 1,"in dlib::nric::lubksb() the indx matrix must be a column vector");
DLIB_ASSERT(b.nc() == 1,"in dlib::nric::lubksb() the b matrix must be a column vector");
const long n = a.nr();
long i, ii = -1, ip, j;
T sum;
for (i = 0; i < n; ++i)
{
ip = indx(i);
sum=b(ip);
b(ip) = b(i);
if (ii != -1)
{
for (j = ii; j < i; ++j)
sum -= a(i,j)*b(j);
}
else if (sum)
{
ii = i;
}
b(i) = sum;
}
for (i = n-1; i >= 0; --i)
{
sum = b(i);
for (j = i+1; j < n; ++j)
sum -= a(i,j)*b(j);
b(i) = sum/a(i,i);
}
}
}
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
template <
typename OP
>
class matrix_zeroary_exp
{
public:
typedef typename OP::type type;
typedef matrix_zeroary_exp ref_type;
typedef typename OP::mem_manager_type mem_manager_type;
const static long NR = OP::NR;
const static long NC = OP::NC;
const typename OP::type operator() (
long r,
long c
) const { return OP::apply(r,c); }
template <typename U, long iNR, long iNC , typename MM>
bool aliases (
const matrix<U,iNR,iNC,MM>& item
) const { return false; }
template <typename U, long iNR, long iNC, typename MM >
bool destructively_aliases (
const matrix<U,iNR,iNC,MM>& item
) const { return false; }
long nr (
) const { return NR; }
long nc (
) const { return NC; }
const ref_type& ref(
) const { return *this; }
};
// ----------------------------------------------------------------------------------------
template <
typename S,
typename OP
>
class dynamic_matrix_scalar_unary_exp
{
/*!
REQUIREMENTS ON S
- must NOT be a matrix_exp or matrix_ref object (or
an object with a compatible interface).
!*/
public:
typedef typename OP::type type;
typedef dynamic_matrix_scalar_unary_exp ref_type;
typedef typename OP::mem_manager_type mem_manager_type;
const static long NR = OP::NR;
const static long NC = OP::NC;
dynamic_matrix_scalar_unary_exp (
long nr__,
long nc__,
const S& s_
) :
nr_(nr__),
nc_(nc__),
s(s_)
{
COMPILE_TIME_ASSERT(is_matrix<S>::value == false);
}
const typename OP::type operator() (
long r,
long c
) const { return OP::apply(s,r,c); }
template <typename U, long iNR, long iNC, typename MM >
bool aliases (
const matrix<U,iNR,iNC,MM>& item
) const { return false; }
template <typename U, long iNR, long iNC , typename MM>
bool destructively_aliases (
const matrix<U,iNR,iNC,MM>& item
) const { return false; }
const ref_type& ref(
) const { return *this; }
long nr (
) const { return nr_; }
long nc (
) const { return nc_; }
private:
const long nr_;
const long nc_;
const S s;
};
// ----------------------------------------------------------------------------------------
template <
typename S,
typename OP
>
class matrix_scalar_unary_exp
{
/*!
REQUIREMENTS ON S
- must NOT be a matrix_exp or matrix_ref object (or
an object with a compatible interface).
!*/
public:
typedef typename OP::type type;
typedef matrix_scalar_unary_exp ref_type;
typedef typename OP::mem_manager_type mem_manager_type;
const static long NR = OP::NR;
const static long NC = OP::NC;
matrix_scalar_unary_exp (
const S& s_
) :
s(s_)
{
COMPILE_TIME_ASSERT(is_matrix<S>::value == false);
}
const typename OP::type operator() (
long r,
long c
) const { return OP::apply(s,r,c); }
template <typename U, long iNR, long iNC, typename MM >
bool aliases (
const matrix<U,iNR,iNC,MM>& item
) const { return false; }
template <typename U, long iNR, long iNC, typename MM >
bool destructively_aliases (
const matrix<U,iNR,iNC,MM>& item
) const { return false; }
const ref_type& ref(
) const { return *this; }
long nr (
) const { return NR; }
long nc (
) const { return NC; }
private:
const S s;
};
// ----------------------------------------------------------------------------------------
template <
typename M,
typename OP_
>
class matrix_unary_exp
{
/*!
REQUIREMENTS ON M
- must be a matrix_exp or matrix_ref object (or
an object with a compatible interface).
!*/
typedef typename OP_::template op<M> OP;
public:
typedef typename OP::type type;
typedef matrix_unary_exp ref_type;
typedef typename OP::mem_manager_type mem_manager_type;
const static long NR = OP::NR;
const static long NC = OP::NC;
matrix_unary_exp (
const M& m_
) :
m(m_)
{}
const typename OP::type operator() (
long r,
long c
) const { return OP::apply(m,r,c); }
template <typename U, long iNR, long iNC, typename MM >
bool aliases (
const matrix<U,iNR,iNC,MM>& item
) const { return m.aliases(item); }
template <typename U, long iNR, long iNC, typename MM >
bool destructively_aliases (
const matrix<U,iNR,iNC,MM>& item
) const { return OP::destructively_aliases(m,item); }
const ref_type& ref(
) const { return *this; }
long nr (
) const { return OP::nr(m); }
long nc (
) const { return OP::nc(m); }
private:
const M m;
};
// ----------------------------------------------------------------------------------------
template <
typename M
>
class matrix_std_vector_exp
{
/*!
REQUIREMENTS ON M
- must be a std::vector object (or
an object with a compatible interface).
!*/
public:
typedef typename M::value_type type;
typedef matrix_std_vector_exp ref_type;
typedef typename memory_manager<char>::kernel_1a mem_manager_type;
const static long NR = 0;
const static long NC = 1;
matrix_std_vector_exp (
const M& m_
) :
m(m_)
{
}
const typename M::value_type operator() (
long r,
long
) const { return m[r]; }
template <typename U, long iNR, long iNC, typename MM>
bool aliases (
const matrix<U,iNR,iNC,MM>& item
) const { return false; }
template <typename U, long iNR, long iNC , typename MM>
bool destructively_aliases (
const matrix<U,iNR,iNC,MM>& item
) const { return false; }
const ref_type& ref(
) const { return *this; }
long nr (
) const { return m.size(); }
long nc (
) const { return 1; }
private:
const M& m;
};
// ----------------------------------------------------------------------------------------
template <
typename M
>
class matrix_vector_exp
{
/*!
REQUIREMENTS ON M
- must be a dlib::array object (or
an object with a compatible interface).
!*/
public:
typedef typename M::type type;
typedef matrix_vector_exp ref_type;
typedef typename M::mem_manager_type mem_manager_type;
const static long NR = 0;
const static long NC = 1;
matrix_vector_exp (
const M& m_
) :
m(m_)
{
}
const typename M::type operator() (
long r,
long
) const { return m[r]; }
template <ty