manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Added method for sorting the singular values computed by the SingularValueDecomposition procedure

This commit is contained in:
Paolo Cignoni 2006-04-29 10:20:52 +00:00
parent 9df4f755ec
commit 26ce24dec3
1 changed files with 104 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

109
vcg/math/lin_algebra.h
View File

@ -150,7 +150,7 @@ namespace vcg
* \return * \return
*/ */
template <typename MATRIX_TYPE> template <typename MATRIX_TYPE>
static bool SingularValueDecomposition(MATRIX_TYPE &A, typename MATRIX_TYPE::ScalarType *W, MATRIX_TYPE &V, const int max_iters = 30) static bool SingularValueDecomposition(MATRIX_TYPE &A, typename MATRIX_TYPE::ScalarType *W, MATRIX_TYPE &V, const bool sort_w=false, const int max_iters=30)
{ {
typedef typename MATRIX_TYPE::ScalarType ScalarType; typedef typename MATRIX_TYPE::ScalarType ScalarType;
int m = (int) A.RowsNumber(); int m = (int) A.RowsNumber();
@ -394,10 +394,105 @@ namespace vcg
} }
} }
delete []rv1; delete []rv1;
if (sort_w)
Sort<MATRIX_TYPE>(A, W, V, false);
return convergence; return convergence;
}; };
/*!
* Sort the singular values computed by the <CODE>SingularValueDecomposition</CODE> procedure and
* modify the matrices <I>U</I> and <I>V</I> accordingly.
*/
template< typename MATRIX_TYPE >
void Sort(MATRIX_TYPE &U, typename MATRIX_TYPE::ScalarType W[], MATRIX_TYPE &V, bool ascending)
{
typedef typename MATRIX_TYPE::ScalarType ScalarType;
assert(U.ColumnsNumber()==V.ColumnsNumber());
int mu = U.RowsNumber();
int mv = V.RowsNumber();
int n = U.ColumnsNumber();
//ScalarType* u = &U[0][0];
//ScalarType* v = &V[0][0];
for (int i=0; i<n; i++)
{
int k = i;
ScalarType p = W[i];
if (ascending)
{
for (int j=i+1; j<n; j++)
{
if (W[j] < p)
{
k = j;
p = W[j];
}
}
}
else
{
for (int j=i+1; j<n; j++)
{
if (W[j] > p)
{
k = j;
p = W[j];
}
}
}
if (k != i)
{
W[k] = W[i]; // i.e.
W[i] = p; // swaps the i-th and the k-th elements
int j = mu;
//ScalarType* uji = u + i; // uji = &U[0][i]
//ScalarType* ujk = u + k; // ujk = &U[0][k]
//ScalarType* vji = v + i; // vji = &V[0][i]
//ScalarType* vjk = v + k; // vjk = &V[0][k]
//if (j)
//{
// for(;;) for( ; j!=0; --j, uji+=n, ujk+=n)
// { {
// p = *uji; p = *uji; // i.e.
// *uji = *ujk; *uji = *ujk; // swap( U[s][i], U[s][k] )
// *ujk = p; *ujk = p; //
// if (!(--j)) }
// break;
// uji += n;
// ujk += n;
// }
//}
for(int s=0; j!=0; ++s, --j)
std::swap(U[s][i], U[s][k]);
j = mv;
//if (j!=0)
//{
// for(;;) for ( ; j!=0; --j, vji+=n, ujk+=n)
// { {
// p = *vji; p = *vji; // i.e.
// *vji = *vjk; *vji = *vjk; // swap( V[s][i], V[s][k] )
// *vjk = p; *vjk = p; //
// if (!(--j)) }
// break;
// vji += n;
// vjk += n;
// }
//}
for (int s=0; j!=0; ++s, --j)
std::swap(V[s][i], V[s][k]);
}
}
}
/*! /*!
* Solves A·X = B for a vector X, where A is specified by the matrices <I>U<SUB>m×n</SUB></I>, * Solves A·X = B for a vector X, where A is specified by the matrices <I>U<SUB>m×n</SUB></I>,
* <I>W<SUB>n×1</SUB></I> and <I>V<SUB>n×n</SUB></I> as returned by <CODE>SingularValueDecomposition</CODE>. * <I>W<SUB>n×1</SUB></I> and <I>V<SUB>n×n</SUB></I> as returned by <CODE>SingularValueDecomposition</CODE>.
@ -414,23 +509,25 @@ namespace vcg
{ {
typedef typename MATRIX_TYPE::ScalarType ScalarType; typedef typename MATRIX_TYPE::ScalarType ScalarType;
unsigned int jj, j, i; unsigned int jj, j, i;
unsigned int columns_number = U.ColumnsNumber();
unsigned int rows_number = U.RowsNumber();
ScalarType s; ScalarType s;
ScalarType *tmp = new ScalarType[U.ColumnsNumber()]; ScalarType *tmp = new ScalarType[columns_number];
for (j=0; j<U.ColumnsNumber(); j++) //Calculate U^T * B. for (j=0; j<columns_number; j++) //Calculate U^T * B.
{ {
s = 0; s = 0;
if (W[j]!=0) //Nonzero result only if wj is nonzero. if (W[j]!=0) //Nonzero result only if wj is nonzero.
{ {
for (i=0; i<U.RowsNumber(); i++) for (i=0; i<rows_number; i++)
s += U[i][j]*b[i]; s += U[i][j]*b[i];
s /= W[j]; //This is the divide by wj . s /= W[j]; //This is the divide by wj .
} }
tmp[j]=s; tmp[j]=s;
} }
for (j=0;j<U.ColumnsNumber();j++) //Matrix multiply by V to get answer. for (j=0;j<columns_number;j++) //Matrix multiply by V to get answer.
{ {
s = 0; s = 0;
for (jj=0; jj<U.ColumnsNumber(); jj++) for (jj=0; jj<columns_number; jj++)
s += V[j][jj]*tmp[jj]; s += V[j][jj]*tmp[jj];
x[j]=s; x[j]=s;
} }