template<class ArrayT>
| ArrayT* Impala::Core::Matrix::MatPInv | ( | ArrayT * | a, | |
| ArrayT * | u, | |||
| ArrayT * | s, | |||
| ArrayT * | v, | |||
| double | tolerance = -1 | |||
| ) | [inline] |
Definition at line 47 of file MatPInv.h.
References MatMul(), MatNrCol(), MatNrRow(), max, min, Impala::Core::Array::PixMax(), Impala::Core::Array::SetVal(), and Impala::Core::Array::Array2dTem< StorT, elemSize, ArithT >::SetValue().
00047 { 00048 00049 Mat* res; 00050 VecScalarReal64* Singulars; 00051 00052 int m = MatNrRow(a); 00053 int n = MatNrCol(a); 00054 00055 int ssize = std::min(s->H(),s->W()); 00056 00057 //Store the singular values in a vector 00058 Singulars = VecCreate<VecScalarReal64>(ssize); 00059 SetVal(Singulars,0); 00060 for(int i=0; i<ssize;i++) 00061 Singulars->SetValue(s->Value(i,i),i,0); 00062 00063 //If No tolerance is given, the default tolerance is: 00064 //eps*max(m,n)*max(Singulars) 00065 if(tolerance <= 0) 00066 tolerance = std::max(m,n) * PixMax(Singulars) * std::numeric_limits<double>::epsilon(); 00067 00068 //Find the number of singular values above tolerance 00069 int countAboveTolerance = 0; 00070 for(int i=0 ; i<ssize ; ++i) 00071 if(Singulars->Value(i,0) > tolerance) 00072 ++countAboveTolerance; 00073 00074 if (countAboveTolerance == 0) 00075 { 00076 res = MatCreate<Mat>(n,m); 00077 SetVal(res, 0); 00078 } 00079 else 00080 { 00081 //SPlus is the square matrix of reciprocals of singular values that are 00082 //above tolerance - it is assumed that the Singular value matrix is 00083 //sorted from largest to smallest singular values. 00084 Mat* SPlus = MatCreate<Mat>(countAboveTolerance,countAboveTolerance); 00085 SetVal(SPlus, 0); 00086 for(int i=0 ; i<countAboveTolerance ; ++i) 00087 SPlus->SetValue(1/Singulars->Value(i,0),i,i); 00088 00089 //We only need the first countAboveTolerance columns from V 00090 Mat* V =MatCreate<Mat>(n,countAboveTolerance); 00091 for(int i=0;i<countAboveTolerance;i++) 00092 for(int j=0;j<n;j++) 00093 V->SetValue(v->Value(i,j),i,j); 00094 00095 //We only need the first countAboveTolerance rows from U' 00096 //Note that U is transposed: 00097 Mat* UT =MatCreate<Mat>(countAboveTolerance,m); 00098 for(int i=0;i<countAboveTolerance;i++) 00099 for(int j=0;j<m;j++) 00100 UT->SetValue(u->Value(i,j),j,i); 00101 00102 Mat* SPlus_X_UT = MatMul(SPlus,UT); 00103 res = MatMul(V,SPlus_X_UT); 00104 //res = MatMul(V,MatMul(SPlus,UT)); 00105 delete V; 00106 delete SPlus; 00107 delete SPlus_X_UT; 00108 delete UT; 00109 } 00110 00111 delete Singulars; 00112 return res; 00113 }
Here is the call graph for this function:
