core/vnl/vnl_rank.txx
Go to the documentation of this file.
00001 // This is core/vnl/vnl_rank.txx
00002 #ifndef vnl_rank_txx_
00003 #define vnl_rank_txx_
00004 
00005 #include "vnl_rank.h"
00006 
00007 template <class T>
00008 vnl_matrix<T> vnl_rank_row_reduce(vnl_matrix<T> const& mat, vnl_rank_pivot_type t)
00009 {
00010   vnl_matrix<T> a = mat;
00011   bool changed = true;
00012   unsigned int m = a.rows(), n=a.columns();
00013   while (changed)
00014   {
00015     changed = false;
00016     for (unsigned int r=0; r<m; ++r)
00017     {
00018       unsigned int c=0; while (c<n && a[r][c] != 1 && a[r][c] != -1) ++c;
00019       if (c==n) continue;
00020       for (unsigned int s=0; s<m; ++s)
00021       {
00022         if (s==r || a[s][c] == 0) continue;
00023         for (unsigned int d=0; d<n; ++d)
00024           if (d!=c) a[s][d] -= a[r][d] * a[r][c] * a[s][c];
00025         a[s][c] = T(0);
00026         changed = true;
00027       }
00028     }
00029   }
00030   if (t == vnl_rank_pivot_one) return a;
00031   changed = true;
00032   while (changed)
00033   {
00034     changed = false;
00035     for (unsigned int r=0; r<m; ++r)
00036     {
00037       unsigned int c=0; while (c<n && a[r][c] == 0) ++c;
00038       if (c==n) continue; // zero row
00039       for (unsigned int s=0; s<m; ++s)
00040       {
00041         if (s==r) continue;
00042         T scale = a[s][c] / a[r][c];
00043         // Note that this can possibly be an integer division, so
00044         // it is *not* guaranteed that a[r][c] * scale == a[s][c] .
00045         if (scale == 0) continue;
00046         for (unsigned int d=0; d<n; ++d)
00047           if (d!=c) a[s][d] -= a[r][d] * scale;
00048         a[s][c] -= a[r][c] * scale;
00049         changed = true;
00050       }
00051     }
00052   }
00053   return a;
00054 }
00055 
00056 template <class T>
00057 vnl_matrix<T> vnl_rank_column_reduce(vnl_matrix<T> const& mat, vnl_rank_pivot_type t)
00058 {
00059   vnl_matrix<T> a = mat;
00060   bool changed = true;
00061   unsigned int m = a.rows(), n=a.columns();
00062   while (changed)
00063   {
00064     changed = false;
00065     for (unsigned int c=0; c<n; ++c)
00066     {
00067       unsigned int r=0; while (r<m && a[r][c] != 1 && a[r][c] != -1) ++r;
00068       if (r==m) continue;
00069       for (unsigned int d=0; d<n; ++d)
00070       {
00071         if (d==c || a[r][d] == 0) continue;
00072         for (unsigned int s=0; s<m; ++s)
00073           if (s!=r) a[s][d] -= a[s][c] * a[r][c] * a[r][d];
00074         a[r][d] = T(0);
00075         changed = true;
00076       }
00077     }
00078   }
00079   if (t == vnl_rank_pivot_one) return a;
00080   changed = true;
00081   while (changed)
00082   {
00083     changed = false;
00084     for (unsigned int c=0; c<n; ++c)
00085     {
00086       unsigned int r=0; while (r<m && a[r][c] == 0) ++r;
00087       if (r==m) continue; // zero row
00088       for (unsigned int d=0; d<n; ++d)
00089       {
00090         if (d==c) continue;
00091         T scale = a[r][d] / a[r][c];
00092         // Note that this can possibly be an integer division, so
00093         // it is *not* guaranteed that a[r][c] * scale == a[r][d] .
00094         if (scale == 0) continue;
00095         for (unsigned int s=0; s<m; ++s)
00096           if (s!=r) a[s][d] -= a[s][c] * scale;
00097         a[r][d] -= a[r][c] * scale;
00098         changed = true;
00099       }
00100     }
00101   }
00102   return a;
00103 }
00104 
00105 template <class T>
00106 vnl_matrix<T> vnl_rank_row_column_reduce(vnl_matrix<T> const& mat, vnl_rank_pivot_type t)
00107 {
00108   vnl_matrix<T> a = mat;
00109   bool changed = true;
00110   unsigned int m = a.rows(), n=a.columns();
00111   while (changed)
00112   {
00113     changed = false;
00114     for (unsigned int r=0; r<m; ++r)
00115     {
00116       unsigned int c=0; while (c<n && a[r][c] != 1 && a[r][c] != -1) ++c;
00117       if (c==n) continue;
00118       for (unsigned int s=0; s<m; ++s)
00119       {
00120         if (s==r || a[s][c] == 0) continue;
00121         for (unsigned int d=0; d<n; ++d)
00122           if (d!=c) a[s][d] -= a[r][d] * a[r][c] * a[s][c];
00123         a[s][c] = T(0);
00124         changed = true;
00125       }
00126     }
00127     for (unsigned int c=0; c<n; ++c)
00128     {
00129       unsigned int r=0; while (r<m && a[r][c] != 1 && a[r][c] != -1) ++r;
00130       if (r==m) continue;
00131       for (unsigned int d=0; d<n; ++d)
00132       {
00133         if (d==c || a[r][d] == 0) continue;
00134         for (unsigned int s=0; s<m; ++s)
00135           if (s!=r) a[s][d] -= a[s][c] * a[r][c] * a[r][d];
00136         a[r][d] = T(0);
00137         changed = true;
00138       }
00139     }
00140   }
00141   if (t == vnl_rank_pivot_one) return a;
00142   changed = true;
00143   while (changed)
00144   {
00145     changed = false;
00146     for (unsigned int r=0; r<m; ++r)
00147     {
00148       unsigned int c=0; while (c<n && a[r][c] == 0) ++c;
00149       if (c==n) continue; // zero row
00150       for (unsigned int s=0; s<m; ++s)
00151       {
00152         if (s==r) continue;
00153         T scale = a[s][c] / a[r][c];
00154         // Note that this can possibly be an integer division, so
00155         // it is *not* guaranteed that a[r][c] * scale == a[s][c] .
00156         if (scale == 0) continue;
00157         for (unsigned int d=0; d<n; ++d)
00158           if (d!=c) a[s][d] -= a[r][d] * scale;
00159         a[s][c] -= a[r][c] * scale;
00160         changed = true;
00161       }
00162     }
00163     for (unsigned int c=0; c<n; ++c)
00164     {
00165       unsigned int r=0; while (r<m && a[r][c] == 0) ++r;
00166       if (r==m) continue; // zero row
00167       for (unsigned int d=0; d<n; ++d)
00168       {
00169         if (d==c) continue;
00170         T scale = a[r][d] / a[r][c];
00171         // Note that this can possibly be an integer division, so
00172         // it is *not* guaranteed that a[r][c] * scale == a[r][d] .
00173         if (scale == 0) continue;
00174         for (unsigned int s=0; s<m; ++s)
00175           if (s!=r) a[s][d] -= a[s][c] * scale;
00176         a[r][d] -= a[r][c] * scale;
00177         changed = true;
00178       }
00179     }
00180   }
00181   return a;
00182 }
00183 
00184 template <class T>
00185 unsigned int vnl_rank(vnl_matrix<T> const& mat, vnl_rank_type t)
00186 {
00187   unsigned int rank = 0;
00188   if (t == vnl_rank_row)
00189   {
00190     vnl_matrix<T> a = vnl_rank_row_reduce(mat, vnl_rank_pivot_all);
00191     for (unsigned int r=0; r<a.rows(); ++r)
00192     {
00193       unsigned int c=0;
00194       while (c<a.columns() && a[r][c] == 0) ++c;
00195       if (c!=a.columns()) ++rank; // not all elements in row r are 0
00196     }
00197   }
00198   else
00199   {
00200     vnl_matrix<T> a = (t == vnl_rank_column) ? vnl_rank_column_reduce(mat,vnl_rank_pivot_all) :
00201                                                vnl_rank_row_column_reduce(mat,vnl_rank_pivot_all);
00202     for (unsigned int c=0; c<a.columns(); ++c)
00203     {
00204       unsigned int r=0;
00205       while (r<a.rows() && a[r][c] == 0) ++r;
00206       if (r!=a.rows()) ++rank; // not all elements in column c are 0
00207     }
00208   }
00209   return rank;
00210 }
00211 
00212 #undef VNL_RANK_INSTANTIATE
00213 #define VNL_RANK_INSTANTIATE(T) \
00214 template vnl_matrix<T > vnl_rank_row_reduce(vnl_matrix<T > const&, vnl_rank_pivot_type);\
00215 template vnl_matrix<T > vnl_rank_column_reduce(vnl_matrix<T > const&, vnl_rank_pivot_type);\
00216 template vnl_matrix<T > vnl_rank_row_column_reduce(vnl_matrix<T > const&, vnl_rank_pivot_type);\
00217 template unsigned int vnl_rank(vnl_matrix<T > const&, vnl_rank_type)
00218 
00219 #endif // vnl_rank_txx_