Mercurial > hg > segmenter-vamp-plugin

annotate armadillo-2.4.4/include/armadillo_bits/op_median_meat.hpp @ 0:8b6102e2a9b0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Armadillo Library
author maxzanoni76 <max.zanoni@eecs.qmul.ac.uk>
date Wed, 11 Apr 2012 09:27:06 +0100
parents
children
rev   line source
max@0 1 // Copyright (C) 2009-2011 NICTA (www.nicta.com.au)
max@0 2 // Copyright (C) 2009-2011 Conrad Sanderson
max@0 3 //
max@0 4 // This file is part of the Armadillo C++ library.
max@0 5 // It is provided without any warranty of fitness
max@0 6 // for any purpose. You can redistribute this file
max@0 7 // and/or modify it under the terms of the GNU
max@0 8 // Lesser General Public License (LGPL) as published
max@0 9 // by the Free Software Foundation, either version 3
max@0 10 // of the License or (at your option) any later version.
max@0 11 // (see http://www.opensource.org/licenses for more info)
max@0 12
max@0 13
max@0 14 //! \addtogroup op_median
max@0 15 //! @{
max@0 16
max@0 17
max@0 18
max@0 19 template<typename eT>
max@0 20 arma_inline
max@0 21 eT
max@0 22 op_median::robust_mean(const eT A, const eT B)
max@0 23 {
max@0 24 return A + (B - A)/eT(2);
max@0 25 }
max@0 26
max@0 27
max@0 28
max@0 29 //! find the median value of a std::vector (contents is modified)
max@0 30 template<typename eT>
max@0 31 inline
max@0 32 eT
max@0 33 op_median::direct_median(std::vector<eT>& X)
max@0 34 {
max@0 35 arma_extra_debug_sigprint();
max@0 36
max@0 37 const uword n_elem = X.size();
max@0 38 const uword half = n_elem/2;
max@0 39
max@0 40 std::sort(X.begin(), X.end());
max@0 41
max@0 42 if((n_elem % 2) == 0)
max@0 43 {
max@0 44 return op_median::robust_mean(X[half-1], X[half]);
max@0 45 }
max@0 46 else
max@0 47 {
max@0 48 return X[half];
max@0 49 }
max@0 50 }
max@0 51
max@0 52
max@0 53
max@0 54 template<typename eT>
max@0 55 inline
max@0 56 eT
max@0 57 op_median::direct_median(const eT* X, const uword n_elem)
max@0 58 {
max@0 59 arma_extra_debug_sigprint();
max@0 60
max@0 61 std::vector<eT> tmp(X, X+n_elem);
max@0 62
max@0 63 return op_median::direct_median(tmp);
max@0 64 }
max@0 65
max@0 66
max@0 67
max@0 68 template<typename eT>
max@0 69 inline
max@0 70 eT
max@0 71 op_median::direct_median(const subview<eT>& X)
max@0 72 {
max@0 73 arma_extra_debug_sigprint();
max@0 74
max@0 75 const uword X_n_elem = X.n_elem;
max@0 76
max@0 77 std::vector<eT> tmp(X_n_elem);
max@0 78
max@0 79 for(uword i=0; i<X_n_elem; ++i)
max@0 80 {
max@0 81 tmp[i] = X[i];
max@0 82 }
max@0 83
max@0 84 return op_median::direct_median(tmp);
max@0 85 }
max@0 86
max@0 87
max@0 88
max@0 89 template<typename eT>
max@0 90 inline
max@0 91 eT
max@0 92 op_median::direct_median(const diagview<eT>& X)
max@0 93 {
max@0 94 arma_extra_debug_sigprint();
max@0 95
max@0 96 const uword X_n_elem = X.n_elem;
max@0 97
max@0 98 std::vector<eT> tmp(X_n_elem);
max@0 99
max@0 100 for(uword i=0; i<X_n_elem; ++i)
max@0 101 {
max@0 102 tmp[i] = X[i];
max@0 103 }
max@0 104
max@0 105 return op_median::direct_median(tmp);
max@0 106 }
max@0 107
max@0 108
max@0 109
max@0 110 //! \brief
max@0 111 //! For each row or for each column, find the median value.
max@0 112 //! The result is stored in a dense matrix that has either one column or one row.
max@0 113 //! The dimension, for which the medians are found, is set via the median() function.
max@0 114 template<typename T1>
max@0 115 inline
max@0 116 void
max@0 117 op_median::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_median>& in)
max@0 118 {
max@0 119 arma_extra_debug_sigprint();
max@0 120
max@0 121 typedef typename T1::elem_type eT;
max@0 122
max@0 123 const unwrap_check<T1> tmp(in.m, out);
max@0 124 const Mat<eT>& X = tmp.M;
max@0 125
max@0 126 const uword X_n_rows = X.n_rows;
max@0 127 const uword X_n_cols = X.n_cols;
max@0 128
max@0 129 const uword dim = in.aux_uword_a;
max@0 130 arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or 1");
max@0 131
max@0 132 if(dim == 0) // in each column
max@0 133 {
max@0 134 arma_extra_debug_print("op_median::apply(), dim = 0");
max@0 135
max@0 136 arma_debug_check( (X_n_rows == 0), "median(): given object has zero rows" );
max@0 137
max@0 138 out.set_size(1, X_n_cols);
max@0 139
max@0 140 std::vector<eT> tmp_vec(X_n_rows);
max@0 141
max@0 142 for(uword col=0; col<X_n_cols; ++col)
max@0 143 {
max@0 144 const eT* colmem = X.colptr(col);
max@0 145
max@0 146 for(uword row=0; row<X_n_rows; ++row)
max@0 147 {
max@0 148 tmp_vec[row] = colmem[row];
max@0 149 }
max@0 150
max@0 151 out[col] = op_median::direct_median(tmp_vec);
max@0 152 }
max@0 153 }
max@0 154 else
max@0 155 if(dim == 1) // in each row
max@0 156 {
max@0 157 arma_extra_debug_print("op_median::apply(), dim = 1");
max@0 158
max@0 159 arma_debug_check( (X_n_cols == 0), "median(): given object has zero columns" );
max@0 160
max@0 161 out.set_size(X_n_rows, 1);
max@0 162
max@0 163 std::vector<eT> tmp_vec(X_n_cols);
max@0 164
max@0 165 for(uword row=0; row<X_n_rows; ++row)
max@0 166 {
max@0 167 for(uword col=0; col<X_n_cols; ++col)
max@0 168 {
max@0 169 tmp_vec[col] = X.at(row,col);
max@0 170 }
max@0 171
max@0 172 out[row] = op_median::direct_median(tmp_vec);
max@0 173 }
max@0 174 }
max@0 175 }
max@0 176
max@0 177
max@0 178
max@0 179 template<typename T>
max@0 180 arma_inline
max@0 181 std::complex<T>
max@0 182 op_median::robust_mean(const std::complex<T>& A, const std::complex<T>& B)
max@0 183 {
max@0 184 return A + (B - A)/T(2);
max@0 185 }
max@0 186
max@0 187
max@0 188
max@0 189 template<typename T>
max@0 190 inline
max@0 191 void
max@0 192 op_median::direct_cx_median_index
max@0 193 (
max@0 194 uword& out_index1,
max@0 195 uword& out_index2,
max@0 196 std::vector< arma_cx_median_packet<T> >& X
max@0 197 )
max@0 198 {
max@0 199 arma_extra_debug_sigprint();
max@0 200
max@0 201 const uword n_elem = X.size();
max@0 202 const uword half = n_elem/2;
max@0 203
max@0 204 std::sort(X.begin(), X.end());
max@0 205
max@0 206 if((n_elem % 2) == 0)
max@0 207 {
max@0 208 out_index1 = X[half-1].index;
max@0 209 out_index2 = X[half ].index;
max@0 210 }
max@0 211 else
max@0 212 {
max@0 213 out_index1 = X[half].index;
max@0 214 out_index2 = out_index1;
max@0 215 }
max@0 216 }
max@0 217
max@0 218
max@0 219
max@0 220 template<typename T>
max@0 221 inline
max@0 222 void
max@0 223 op_median::direct_cx_median_index
max@0 224 (
max@0 225 uword& out_index1,
max@0 226 uword& out_index2,
max@0 227 const std::complex<T>* X,
max@0 228 const uword n_elem
max@0 229 )
max@0 230 {
max@0 231 arma_extra_debug_sigprint();
max@0 232
max@0 233 std::vector< arma_cx_median_packet<T> > tmp(n_elem);
max@0 234
max@0 235 for(uword i=0; i<n_elem; ++i)
max@0 236 {
max@0 237 tmp[i].val = std::abs(X[i]);
max@0 238 tmp[i].index = i;
max@0 239 }
max@0 240
max@0 241 op_median::direct_cx_median_index(out_index1, out_index2, tmp);
max@0 242 }
max@0 243
max@0 244
max@0 245
max@0 246 template<typename T>
max@0 247 inline
max@0 248 void
max@0 249 op_median::direct_cx_median_index
max@0 250 (
max@0 251 uword& out_index1,
max@0 252 uword& out_index2,
max@0 253 const subview< std::complex<T> >&X
max@0 254 )
max@0 255 {
max@0 256 arma_extra_debug_sigprint();
max@0 257
max@0 258 const uword n_elem = X.n_elem;
max@0 259
max@0 260 std::vector< arma_cx_median_packet<T> > tmp(n_elem);
max@0 261
max@0 262 for(uword i=0; i<n_elem; ++i)
max@0 263 {
max@0 264 tmp[i].val = std::abs(X[i]);
max@0 265 tmp[i].index = i;
max@0 266 }
max@0 267
max@0 268 op_median::direct_cx_median_index(out_index1, out_index2, tmp);
max@0 269 }
max@0 270
max@0 271
max@0 272
max@0 273 template<typename T>
max@0 274 inline
max@0 275 void
max@0 276 op_median::direct_cx_median_index
max@0 277 (
max@0 278 uword& out_index1,
max@0 279 uword& out_index2,
max@0 280 const diagview< std::complex<T> >&X
max@0 281 )
max@0 282 {
max@0 283 arma_extra_debug_sigprint();
max@0 284
max@0 285 const uword n_elem = X.n_elem;
max@0 286
max@0 287 std::vector< arma_cx_median_packet<T> > tmp(n_elem);
max@0 288
max@0 289 for(uword i=0; i<n_elem; ++i)
max@0 290 {
max@0 291 tmp[i].val = std::abs(X[i]);
max@0 292 tmp[i].index = i;
max@0 293 }
max@0 294
max@0 295 op_median::direct_cx_median_index(out_index1, out_index2, tmp);
max@0 296 }
max@0 297
max@0 298
max@0 299
max@0 300 //! Implementation for complex numbers
max@0 301 template<typename T, typename T1>
max@0 302 inline
max@0 303 void
max@0 304 op_median::apply(Mat< std::complex<T> >& out, const Op<T1,op_median>& in)
max@0 305 {
max@0 306 arma_extra_debug_sigprint();
max@0 307
max@0 308 typedef typename std::complex<T> eT;
max@0 309
max@0 310 arma_type_check(( is_same_type<eT, typename T1::elem_type>::value == false ));
max@0 311
max@0 312 const unwrap_check<T1> tmp(in.m, out);
max@0 313 const Mat<eT>& X = tmp.M;
max@0 314
max@0 315 const uword X_n_rows = X.n_rows;
max@0 316 const uword X_n_cols = X.n_cols;
max@0 317
max@0 318 const uword dim = in.aux_uword_a;
max@0 319 arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or 1");
max@0 320
max@0 321 if(dim == 0) // in each column
max@0 322 {
max@0 323 arma_extra_debug_print("op_median::apply(), dim = 0");
max@0 324
max@0 325 arma_debug_check( (X_n_rows == 0), "median(): given object has zero rows" );
max@0 326
max@0 327 out.set_size(1, X_n_cols);
max@0 328
max@0 329 std::vector< arma_cx_median_packet<T> > tmp_vec(X_n_rows);
max@0 330
max@0 331 for(uword col=0; col<X_n_cols; ++col)
max@0 332 {
max@0 333 const eT* colmem = X.colptr(col);
max@0 334
max@0 335 for(uword row=0; row<X_n_rows; ++row)
max@0 336 {
max@0 337 tmp_vec[row].val = std::abs(colmem[row]);
max@0 338 tmp_vec[row].index = row;
max@0 339 }
max@0 340
max@0 341 uword index1;
max@0 342 uword index2;
max@0 343 op_median::direct_cx_median_index(index1, index2, tmp_vec);
max@0 344
max@0 345 out[col] = op_median::robust_mean(colmem[index1], colmem[index2]);
max@0 346 }
max@0 347 }
max@0 348 else
max@0 349 if(dim == 1) // in each row
max@0 350 {
max@0 351 arma_extra_debug_print("op_median::apply(), dim = 1");
max@0 352
max@0 353 arma_debug_check( (X_n_cols == 0), "median(): given object has zero columns" );
max@0 354
max@0 355 out.set_size(X_n_rows, 1);
max@0 356
max@0 357 std::vector< arma_cx_median_packet<T> > tmp_vec(X_n_cols);
max@0 358
max@0 359 for(uword row=0; row<X_n_rows; ++row)
max@0 360 {
max@0 361 for(uword col=0; col<X_n_cols; ++col)
max@0 362 {
max@0 363 tmp_vec[col].val = std::abs(X.at(row,col));
max@0 364 tmp_vec[row].index = col;
max@0 365 }
max@0 366
max@0 367 uword index1;
max@0 368 uword index2;
max@0 369 op_median::direct_cx_median_index(index1, index2, tmp_vec);
max@0 370
max@0 371 out[row] = op_median::robust_mean( X.at(row,index1), X.at(row,index2) );
max@0 372 }
max@0 373 }
max@0 374 }
max@0 375
max@0 376
max@0 377
max@0 378 //! @}
max@0 379