00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 template<typename eT>
00024 inline
00025 eT
00026 op_median::direct_median(std::vector<eT>& X)
00027 {
00028 arma_extra_debug_sigprint();
00029
00030 std::sort(X.begin(), X.end());
00031
00032 const u32 n_elem = X.size();
00033 const u32 half = n_elem/2;
00034
00035 if((n_elem % 2) == 0)
00036 {
00037 return (X[half-1] + X[half]) / eT(2);
00038 }
00039 else
00040 {
00041 return X[half];
00042 }
00043 }
00044
00045
00046
00047 template<typename eT>
00048 inline
00049 eT
00050 op_median::direct_median(const eT* X, const u32 n_elem)
00051 {
00052 arma_extra_debug_sigprint();
00053
00054 std::vector<eT> tmp(X, X+n_elem);
00055 return op_median::direct_median(tmp);
00056 }
00057
00058
00059
00060 template<typename eT>
00061 inline
00062 eT
00063 op_median::direct_median(const subview<eT>& X)
00064 {
00065 arma_extra_debug_sigprint();
00066
00067 std::vector<eT> tmp(X.n_elem);
00068
00069 for(u32 i=0; i<X.n_elem; ++i)
00070 {
00071 tmp[i] = X[i];
00072 }
00073
00074 return op_median::direct_median(tmp);
00075 }
00076
00077
00078
00079 template<typename eT>
00080 inline
00081 eT
00082 op_median::direct_median(const diagview<eT>& X)
00083 {
00084 arma_extra_debug_sigprint();
00085
00086 std::vector<eT> tmp(X.n_elem);
00087
00088 for(u32 i=0; i<X.n_elem; ++i)
00089 {
00090 tmp[i] = X[i];
00091 }
00092
00093 return op_median::direct_median(tmp);
00094 }
00095
00096
00097
00098
00099
00100
00101
00102 template<typename T1>
00103 inline
00104 void
00105 op_median::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_median>& in)
00106 {
00107 arma_extra_debug_sigprint();
00108
00109 typedef typename T1::elem_type eT;
00110
00111 const unwrap_check<T1> tmp(in.m, out);
00112 const Mat<eT>& X = tmp.M;
00113
00114 arma_debug_check( (X.n_elem == 0), "median(): given matrix has no elements" );
00115
00116 const u32 dim = in.aux_u32_a;
00117 arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or 1");
00118
00119
00120 if(dim == 0)
00121 {
00122 arma_extra_debug_print("op_median::apply(), dim = 0");
00123
00124 out.set_size(1, X.n_cols);
00125
00126 std::vector<eT> tmp_vec(X.n_rows);
00127
00128 for(u32 col=0; col<X.n_cols; ++col)
00129 {
00130 const eT* colmem = X.colptr(col);
00131
00132 for(u32 row=0; row<X.n_rows; ++row)
00133 {
00134 tmp_vec[row] = colmem[row];
00135 }
00136
00137 out[col] = op_median::direct_median(tmp_vec);
00138 }
00139 }
00140 else
00141 if(dim == 1)
00142 {
00143 arma_extra_debug_print("op_median::apply(), dim = 1");
00144
00145 out.set_size(X.n_rows, 1);
00146
00147 std::vector<eT> tmp_vec(X.n_cols);
00148
00149 for(u32 row=0; row<X.n_rows; ++row)
00150 {
00151 for(u32 col=0; col<X.n_cols; ++col)
00152 {
00153 tmp_vec[col] = X.at(row,col);
00154 }
00155
00156 out[row] = op_median::direct_median(tmp_vec);
00157 }
00158 }
00159
00160 }
00161
00162
00163
00164 template<typename T>
00165 inline
00166 void
00167 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, std::vector< arma_cx_median_packet<T> >& X)
00168 {
00169 arma_extra_debug_sigprint();
00170
00171 std::sort(X.begin(), X.end());
00172
00173 const u32 n_elem = X.size();
00174 const u32 half = n_elem/2;
00175
00176 if((n_elem % 2) == 0)
00177 {
00178 out_index1 = X[half-1].index;
00179 out_index2 = X[half].index;
00180 }
00181 else
00182 {
00183 out_index1 = X[half].index;
00184 out_index2 = X[half].index;
00185 }
00186
00187 }
00188
00189
00190
00191 template<typename T>
00192 inline
00193 void
00194 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, const std::complex<T>* X, const u32 n_elem)
00195 {
00196 arma_extra_debug_sigprint();
00197
00198 std::vector< arma_cx_median_packet<T> > tmp(n_elem);
00199
00200 for(u32 i=0; i<n_elem; ++i)
00201 {
00202 tmp[i].val = std::abs(X[i]);
00203 tmp[i].index = i;
00204 }
00205
00206 op_median::direct_cx_median_index(out_index1, out_index2, tmp);
00207 }
00208
00209
00210
00211 template<typename T>
00212 inline
00213 void
00214 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, const subview< std::complex<T> >&X)
00215 {
00216 arma_extra_debug_sigprint();
00217
00218 const u32 n_elem = X.n_elem;
00219
00220 std::vector< arma_cx_median_packet<T> > tmp(n_elem);
00221
00222 for(u32 i=0; i<n_elem; ++i)
00223 {
00224 tmp[i].val = std::abs(X[i]);
00225 tmp[i].index = i;
00226 }
00227
00228 op_median::direct_cx_median_index(out_index1, out_index2, tmp);
00229 }
00230
00231
00232
00233 template<typename T>
00234 inline
00235 void
00236 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, const diagview< std::complex<T> >&X)
00237 {
00238 arma_extra_debug_sigprint();
00239
00240 const u32 n_elem = X.n_elem;
00241
00242 std::vector< arma_cx_median_packet<T> > tmp(n_elem);
00243
00244 for(u32 i=0; i<n_elem; ++i)
00245 {
00246 tmp[i].val = std::abs(X[i]);
00247 tmp[i].index = i;
00248 }
00249
00250 op_median::direct_cx_median_index(out_index1, out_index2, tmp);
00251 }
00252
00253
00254
00255
00256 template<typename T, typename T1>
00257 inline
00258 void
00259 op_median::apply(Mat< std::complex<T> >& out, const Op<T1,op_median>& in)
00260 {
00261 arma_extra_debug_sigprint();
00262
00263 typedef typename std::complex<T> eT;
00264 isnt_same_type<eT, typename T1::elem_type>::check();
00265
00266 const unwrap_check<T1> tmp(in.m, out);
00267 const Mat<eT>& X = tmp.M;
00268
00269 arma_debug_check( (X.n_elem == 0), "median(): given matrix has no elements" );
00270
00271 const u32 dim = in.aux_u32_a;
00272 arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or 1");
00273
00274
00275 if(dim == 0)
00276 {
00277 arma_extra_debug_print("op_median::apply(), dim = 0");
00278
00279 out.set_size(1, X.n_cols);
00280
00281 std::vector< arma_cx_median_packet<T> > tmp_vec(X.n_rows);
00282
00283 for(u32 col=0; col<X.n_cols; ++col)
00284 {
00285 const eT* colmem = X.colptr(col);
00286
00287 for(u32 row=0; row<X.n_rows; ++row)
00288 {
00289 tmp_vec[row].val = std::abs(colmem[row]);
00290 tmp_vec[row].index = row;
00291 }
00292
00293 u32 index1;
00294 u32 index2;
00295 op_median::direct_cx_median_index(index1, index2, tmp_vec);
00296
00297 out[col] = (colmem[index1] + colmem[index2]) / T(2);
00298 }
00299 }
00300 else
00301 if(dim == 1)
00302 {
00303 arma_extra_debug_print("op_median::apply(), dim = 1");
00304
00305 out.set_size(X.n_rows, 1);
00306
00307 std::vector< arma_cx_median_packet<T> > tmp_vec(X.n_cols);
00308
00309 for(u32 row=0; row<X.n_rows; ++row)
00310 {
00311 for(u32 col=0; col<X.n_cols; ++col)
00312 {
00313 tmp_vec[col].val = std::abs(X.at(row,col));
00314 tmp_vec[row].index = col;
00315 }
00316
00317 u32 index1;
00318 u32 index2;
00319 op_median::direct_cx_median_index(index1, index2, tmp_vec);
00320
00321 out[row] = ( X.at(row,index1) + X.at(row,index2) ) / T(2);
00322 }
00323 }
00324
00325 }
00326
00327
00328
00329