Go to the documentation of this file.00001
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;
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
00044
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;
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
00093
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;
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
00155
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;
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
00172
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;
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;
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_