[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/accumulator.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2011-2012 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 #ifndef VIGRA_ACCUMULATOR_HXX 00037 #define VIGRA_ACCUMULATOR_HXX 00038 00039 #ifdef _MSC_VER 00040 #pragma warning (disable: 4503) 00041 #endif 00042 00043 #include "accumulator-grammar.hxx" 00044 #include "config.hxx" 00045 #include "metaprogramming.hxx" 00046 #include "bit_array.hxx" 00047 #include "static_assert.hxx" 00048 #include "mathutil.hxx" 00049 #include "utilities.hxx" 00050 #include "multi_iterator_coupled.hxx" 00051 #include "matrix.hxx" 00052 #include "multi_math.hxx" 00053 #include "eigensystem.hxx" 00054 #include "histogram.hxx" 00055 #include <algorithm> 00056 #include <iostream> 00057 00058 namespace vigra { 00059 00060 /** \defgroup FeatureAccumulators Feature Accumulators 00061 00062 The namespace <tt>vigra::acc</tt> provides the function \ref vigra::acc::extractFeatures() along with associated statistics functors and accumulator classes. Together, they provide a framework for efficient compution of a wide variety of statistical features, both globally for an entire image, and locally for each region defined by a label array. Many different statistics can be composed out of a small number of fundamental statistics and suitable modifiers. The user simply selects the desired statistics by means of their <i>tags</i> (see below), and a template meta-program automatically generates an efficient functor that computes exactly those statistics. 00063 00064 The function \ref acc::extractFeatures() "extractFeatures()" scans the data in as few passes as the selected statstics permit (usually one or two passes are sufficient). Statistics are computed by accurate incremental algorithms, whose internal state is maintained by accumulator objects. The state is updated by passing data to the accumulator one sample at a time. Accumulators are grouped within an accumulator chain. Dependencies between accumulators in the accumulator chain are automatically resolved and missing dependencies are inserted. For example, to compute the mean, you also need to count the number of samples. This allows accumulators to offload some of their computations on other accumulators, making the algorithms more efficient. Each accumulator only sees data in the appropriate pass through the data, called its "working pass". 00065 00066 <b>\#include</b> <vigra/accumulator.hxx> 00067 00068 <b>Basic statistics:</b> 00069 - PowerSum<N> (@f$ \sum_i x_i^N @f$) 00070 - AbsPowerSum<N> (@f$ \sum_i |x_i|^N @f$) 00071 - Skewness, UnbiasedSkewness 00072 - Kurtosis, UnbiasedKurtosis 00073 - Minimum, Maximum 00074 - FlatScatterMatrix (flattened upper-triangular part of scatter matrix) 00075 - 4 histogram classes (see \ref histogram "below") 00076 - StandardQuantiles (0%, 10%, 25%, 50%, 75%, 90%, 100%) 00077 - ArgMinWeight, ArgMaxWeight (store data or coordinate where weight assumes its minimal or maximal value) 00078 - CoordinateSystem (identity matrix of appropriate size) 00079 00080 <b>Modifiers:</b> (S is the statistc to be modified) 00081 - Normalization 00082 <table border="0"> 00083 <tr><td> DivideByCount<S> </td><td> S/Count </td></tr> 00084 <tr><td> RootDivideByCount<S> </td><td> sqrt( S/Count ) </td></tr> 00085 <tr><td> DivideUnbiased<S> </td><td> S/(Count-1) </td></tr> 00086 <tr><td> RootDivideUnbiased<S> </td><td> sqrt( S/(Count-1) ) </td></tr> 00087 </table> 00088 - Data preparation: 00089 <table border="0"> 00090 <tr><td> Central<S> </td><td> substract mean before computing S </td></tr> 00091 <tr><td> Principal<S> </td><td> project onto PCA eigenvectors </td></tr> 00092 <tr><td> Whitened<S> </td><td> scale to unit variance after PCA </td></tr> 00093 <tr><td> Coord<S> </td><td> compute S from pixel coordinates rather than from pixel values </td></tr> 00094 <tr><td> Weighted<S> </td><td> compute weighted version of S </td></tr> 00095 <tr><td> Global<S> </td><td> compute S globally rather than per region (per region is default if labels are given) </td></tr> 00096 </table> 00097 00098 Aliases for a couple of important features are implemented (mainly as <tt>typedef FullName Alias</tt>). The alias names are equivalent to full names. 00099 Here are some examples for supported alias names (these examples also show how to compose statistics from the fundamental statistics and modifiers): 00100 00101 <table border="0"> 00102 <tr><th> Alias </th><th> Full Name </th></tr> 00103 <tr><td> Count </td><td> PowerSum<0> </td></tr> 00104 <tr><td> Sum </td><td> PowerSum<1> </td></tr> 00105 <tr><td> SumOfSquares </td><td> PowerSum<2> </td></tr> 00106 <tr><td> Mean </td><td> DivideByCount<PowerSum<1>> </td></tr> 00107 <tr><td> RootMeanSquares </td><td> RootDivideByCount<PowerSum<2>> </td></tr> 00108 <tr><td> Moment<N> </td><td> DivideByCount<PowerSum<N>> </td></tr> 00109 <tr><td> Variance </td><td> DivideByCount<Central<PowerSum<2>>> </td></tr> 00110 <tr><td> StdDev </td><td> RootDivideByCount<Central<PowerSum<2>>> </td></tr> 00111 <tr><td> Covariance </td><td> DivideByCount<FlatScatterMatrix> </td></tr> 00112 <tr><td> RegionCenter </td><td> Coord<Mean> </td></tr> 00113 <tr><td> CenterOfMass </td><td> Weighted<Coord<Mean>> </td></tr> 00114 </table> 00115 00116 There are a few <b>rules for composing statistics</b>: 00117 - modifiers can be specified in any order, but are internally transformed to standard order: Global<Weighted<Coord<normalization<data preparation<basic statistic 00118 - only one normalization modifier and one data preparation modifier (Central or Principal or Whitened) is permitted 00119 - Count ignores all modifiers except Global and Weighted 00120 - Sum ignores Central and Principal, because sum would be zero 00121 - ArgMinWeight and ArgMaxWeight are automatically Weighted 00122 00123 00124 Here is an example how to use \ref acc::AccumulatorChain to compute statistics. (To use Weighted<> or Coord<> modifiers, see below): 00125 00126 \code 00127 #include <vigra/multi_array.hxx> 00128 #include <vigra/impex.hxx> 00129 #include <vigra/accumulator.hxx> 00130 using namespace vigra::acc; 00131 typedef double DataType; 00132 int size = 1000; 00133 vigra::MultiArray<2, DataType> data(vigra::Shape2(size, size)); 00134 00135 AccumulatorChain<DataType, 00136 Select<Variance, Mean, StdDev, Minimum, Maximum, RootMeanSquares, Skewness, Covariance> > 00137 a; 00138 00139 std::cout << "passes required: " << a.passesRequired() << std::endl; 00140 extractFeatures(data.begin(), data.end(), a); 00141 00142 std::cout << "Mean: " << get<Mean>(a) << std::endl; 00143 std::cout << "Variance: " << get<Variance>(a) << std::endl; 00144 \endcode 00145 00146 The \ref acc::AccumulatorChain object contains the selected statistics and their dependencies. Statistics have to be wrapped with \ref acc::Select. The statistics are computed with the acc::extractFeatures function and the statistics can be accessed with acc::get . 00147 00148 Rules and notes: 00149 - order of statistics in Select<> is arbitrary 00150 - up to 20 statistics in Select<>, but Select<> can be nested 00151 - dependencies are automatically inserted 00152 - duplicates are automatically removed 00153 - extractFeatures() does as many passes through the data as necessary 00154 - each accumulator only sees data in the appropriate pass (its "working pass") 00155 00156 The Accumulators can also be used with vector-valued data (vigra::RGBValue, vigra::TinyVector, vigra::MultiArray or vigra::MultiArrayView): 00157 00158 \code 00159 typedef vigra::RGBValue<double> DataType; 00160 AccumulatorChain<DataType, Select<...> > a; 00161 ... 00162 \endcode 00163 00164 To compute <b>weighted statistics</b> (Weighted<>) or <b>statistics over coordinates</b> (Coord<>), the accumulator chain can be used with \ref CoupledScanOrderIterator. The coupled iterator provides simultaneous access to several images (e.g. weight and data) and pixel coordinates. The first parameter in the accumulator chain is the type of the CoupledHandle. The indeces at which the CoupledHandle holds the data, weights etc. can be specified inside the Select wrapper. 00165 00166 These <b>index specifiers</b> are: (INDEX is of type int) 00167 00168 - DataArg<INDEX>: CoupledHandle holds data at index 'INDEX' (default INDEX=1) 00169 - LabelArg<INDEX>: CoupledHandle holds labels at index 'INDEX' (default INDEX=2) 00170 - WeightArg<INDEX>: CoupledHandle holds weights at index 'INDEX' (default INDEX=outermost index) 00171 00172 Pixel coordinates are always at index 0. 00173 00174 \code 00175 using namespace vigra::acc; 00176 vigra::MultiArray<3, double> data(...), weights(...); 00177 typedef vigra::CoupledIteratorType<3, double, double>::type Iterator; //type of the CoupledScanOrderIterator 00178 typedef Iterator::value_type Handle; //type of the corresponding CoupledHandle 00179 00180 AccumulatorChain<Handle, 00181 Select<DataArg<1>, WeightArg<2>, //where to look in the Handle (coordinates are always arg 0) 00182 Mean, Variance, //statistics over values 00183 Coord<Mean>, Coord<Variance>, //statistics over coordinates, 00184 Weighted<Mean>, Weighted<Variance>, //weighted values, 00185 Weighted<Coord<Mean> > > > //weighted coordinates. 00186 a; 00187 00188 Iterator start = createCoupledIterator(data, weights); //coord->index 0, data->index 1, weights->index 2 00189 Iterator end = start.getEndIterator(); 00190 00191 extractFeatures(start,end,a); 00192 \endcode 00193 00194 To compute <b>region statistics</b>, use \ref acc::AccumulatorChainArray : 00195 00196 \code 00197 using namespace vigra::acc; 00198 vigra::MultiArray<3, double> data(...); 00199 vigra::MultiArray<3, int> labels(...); 00200 typedef vigra::CoupledIteratorType<3, double, int>::type Iterator; 00201 typedef Iterator::value_type Handle; 00202 00203 AccumulatorChainArray<Handle, 00204 Select<DataArg<1>, LabelArg<2>, //where to look in the Handle (coordinates are always arg 0) 00205 Mean, Variance, //per-region statistics over values 00206 Coord<Mean>, Coord<Variance>, //per-region statistics over coordinates 00207 Global<Mean>, Global<Variance> > > //global statistics 00208 a; 00209 00210 Iterator start = createCoupledIterator(data, labels); 00211 Iterator end = start.getEndIterator(); 00212 00213 a.ignoreLabel(0); //statistics will not be computed for region 0 (e.g. background) 00214 00215 extractFeatures(start,end,a); 00216 00217 int regionlabel = ...; 00218 std::cout << get<Mean>(a, regionlabel) << std::endl; //get Mean of region with label 'regionlabel' 00219 \endcode 00220 00221 00222 In some application it will be known only at run-time which statistics have to be computed. An Accumulator with <b>run-time activation</b> is provided by the \ref acc::DynamicAccumulatorChain class. One specifies a set of statistics at compile-time and from this set one can activate the needed statistics at run-time: 00223 00224 \code 00225 using namespace vigra::acc; 00226 vigra::MultiArray<2, double> data(...); 00227 DynamicAccumulatorChain<double, 00228 Select<Mean, Minimum, Maximum, Variance, StdDev> > a; // at compile-time 00229 activate<Mean>(a); //at run-time 00230 a.activate("Minimum"); //same as activate<Minimum>(a) (alias names are not recognized) 00231 00232 extractFeatures(data.begin(), data.end(), a); 00233 std::cout << "Mean: " << get<Mean>(a) << std::endl; //ok 00234 //std::cout << "Maximum: " << get<Maximum>(a) << std::endl; // run-time error because Maximum not activated 00235 \endcode 00236 00237 Likewise, for run-time activation of region statistics, use \ref acc::DynamicAccumulatorChainArray. 00238 00239 <b>Accumulator merging</b> (e.g. for parallelization or hierarchical segmentation) is possible for many accumulators: 00240 00241 \code 00242 using namespace vigra::acc; 00243 vigra::MultiArray<2, double> data(...); 00244 AccumulatorChain<double, Select<Mean, Variance, Skewness> > a, a1, a2; 00245 00246 extractFeatures(data.begin(), data.end(), a); //process entire data set at once 00247 extractFeatures(data.begin(), data.begin()+data.size()/2, a1); //process first half 00248 extractFeatures(data.begin()+data.size()/2, data.end(), a2); //process second half 00249 a1 += a2; // merge: a1 now equals a0 (with numerical tolerances) 00250 \endcode 00251 00252 Not all statistics can be merged (e.g. Principal<A> usually cannot, except for some important specializations). A statistic can be merged if the "+=" operator is supported (see the documentation of that particular statistic). If the accumulator chain only requires one pass to collect the data, it is also possible to just apply the extractFeatures() function repeatedly: 00253 00254 \code 00255 using namespace vigra::acc; 00256 vigra::MultiArray<2, double> data(...); 00257 AccumulatorChain<double, Select<Mean, Variance> > a; 00258 00259 extractFeatures(data.begin(), data.begin()+data.size()/2, a); // this works because 00260 extractFeatures(data.begin()+data.size()/2, data.end(), a); // all statistics only work in pass 1 00261 00262 \endcode 00263 00264 00265 \anchor histogram 00266 Four kinds of <b>histograms</b> are currently implemented: 00267 00268 <table border="0"> 00269 <tr><td> IntegerHistogram </td><td> Data values are equal to bin indices </td></tr> 00270 <tr><td> UserRangeHistogram </td><td> User provides lower and upper bounds for linear range mapping from values to indices. </td></tr> 00271 <tr><td> AutoRangeHistogram </td><td> Range mapping bounds are defiend by minimum and maximum of the data (2 passes needed!) </td></tr> 00272 <tr><td> GlobalRangeHistogram </td><td> Likewise, but use global min/max rather than region min/max as AutoRangeHistogram will </td></tr> 00273 </table> 00274 00275 00276 00277 - The number of bins is specified at compile time (as template parameter int BinCount) or at run-time (if BinCount is zero at compile time). In the first case the return type of the accumulator is TinyVector<double, BinCount> (number of bins cannot be changed). In the second case, the return type is MultiArray<1, double> and the number of bins must be set before seeing data (see example below). 00278 - If UserRangeHistogram is used, the bounds for the linear range mapping from values to indices must be set before seeing data (see below). 00279 - Options can be set by passing an instance of HistogramOptions to the accumulator chain (same options for all histograms in the chain) or by directly calling the appropriate member functions of the accumulators. 00280 - Merging is supported if the range mapping of the histograms is the same. 00281 - Histogram accumulators have two members for outliers (left_outliers, right_outliers). 00282 00283 With the StandardQuantiles class, <b>histogram quantiles</b> (0%, 10%, 25%, 50%, 75%, 90%, 100%) are computed from a given histgram using linear interpolation. The return type is TinyVector<double, 7> . 00284 00285 \anchor acc_hist_options Usage: 00286 \code 00287 using namespace vigra::acc; 00288 typedef double DataType; 00289 vigra::MultiArray<2, DataType> data(...); 00290 00291 typedef UserRangeHistogram<40> SomeHistogram; //binCount set at compile time 00292 typedef UserRangeHistogram<0> SomeHistogram2; // binCount must be set at run-time 00293 typedef AutoRangeHistogram<0> SomeHistogram3; 00294 typedef StandardQuantiles<SomeHistogram3> Quantiles3; 00295 00296 AccumulatorChain<DataType, Select<SomeHistogram, SomeHistogram2, SomeHistogram3, Quantiles3> > a; 00297 00298 //set options for all histograms in the accumulator chain: 00299 vigra::HistogramOptions histogram_opt; 00300 histogram_opt = histogram_opt.setBinCount(50); 00301 //histogram_opt = histogram_opt.setMinMax(0.1, 0.9); // this would set min/max for all three histograms, but range bounds 00302 // shall be set automatically by min/max of data for SomeHistogram3 00303 a.setHistogramOptions(histogram_opt); 00304 00305 // set options for a specific histogram in the accumulator chain: 00306 getAccumulator<SomeHistogram>(a).setMinMax(0.1, 0.9); // number of bins must be set before setting min/max 00307 getAccumulator<SomeHistogram2>(a).setMinMax(0.0, 1.0); 00308 00309 extractFeatures(data.begin(), data.end(), a); 00310 00311 vigra::TinyVector<double, 40> hist = get<SomeHistogram>(a); 00312 vigra::MultiArray<1, double> hist2 = get<SomeHistogram2>(a); 00313 vigra::TinyVector<double, 7> quant = get<Quantiles3>(a); 00314 double right_outliers = getAccumulator<SomeHistogram>(a).right_outliers; 00315 \endcode 00316 00317 00318 00319 */ 00320 00321 00322 /** This namespace contains the accumulator classes, fundamental statistics and modifiers. See \ref FeatureAccumulators for examples of usage. 00323 */ 00324 namespace acc { 00325 00326 /****************************************************************************/ 00327 /* */ 00328 /* infrastructure */ 00329 /* */ 00330 /****************************************************************************/ 00331 00332 /// \brief Wrapper for MakeTypeList that additionally performs tag standardization. 00333 00334 template <class T01=void, class T02=void, class T03=void, class T04=void, class T05=void, 00335 class T06=void, class T07=void, class T08=void, class T09=void, class T10=void, 00336 class T11=void, class T12=void, class T13=void, class T14=void, class T15=void, 00337 class T16=void, class T17=void, class T18=void, class T19=void, class T20=void> 00338 struct Select 00339 : public MakeTypeList< 00340 typename StandardizeTag<T01>::type, typename StandardizeTag<T02>::type, typename StandardizeTag<T03>::type, 00341 typename StandardizeTag<T04>::type, typename StandardizeTag<T05>::type, typename StandardizeTag<T06>::type, 00342 typename StandardizeTag<T07>::type, typename StandardizeTag<T08>::type, typename StandardizeTag<T09>::type, 00343 typename StandardizeTag<T10>::type, typename StandardizeTag<T11>::type, typename StandardizeTag<T12>::type, 00344 typename StandardizeTag<T13>::type, typename StandardizeTag<T14>::type, typename StandardizeTag<T15>::type, 00345 typename StandardizeTag<T16>::type, typename StandardizeTag<T17>::type, typename StandardizeTag<T18>::type, 00346 typename StandardizeTag<T19>::type, typename StandardizeTag<T20>::type 00347 > 00348 {}; 00349 00350 // enable nesting of Select<> expressions 00351 template <class T01, class T02, class T03, class T04, class T05, 00352 class T06, class T07, class T08, class T09, class T10, 00353 class T11, class T12, class T13, class T14, class T15, 00354 class T16, class T17, class T18, class T19, class T20> 00355 struct StandardizeTag<Select<T01, T02, T03, T04, T05, 00356 T06, T07, T08, T09, T10, 00357 T11, T12, T13, T14, T15, 00358 T16, T17, T18, T19, T20>, 00359 Select<T01, T02, T03, T04, T05, 00360 T06, T07, T08, T09, T10, 00361 T11, T12, T13, T14, T15, 00362 T16, T17, T18, T19, T20> > 00363 { 00364 typedef typename Select<T01, T02, T03, T04, T05, 00365 T06, T07, T08, T09, T10, 00366 T11, T12, T13, T14, T15, 00367 T16, T17, T18, T19, T20>::type type; 00368 }; 00369 00370 struct AccumulatorBegin 00371 { 00372 typedef Select<> Dependencies; 00373 00374 static std::string const & name() 00375 { 00376 static const std::string n("AccumulatorBegin (internal)"); 00377 return n; 00378 } 00379 00380 template <class T, class BASE> 00381 struct Impl 00382 : public BASE 00383 {}; 00384 }; 00385 00386 00387 struct AccumulatorEnd; 00388 struct DataArgTag; 00389 struct WeightArgTag; 00390 struct LabelArgTag; 00391 struct CoordArgTag; 00392 struct LabelDispatchTag; 00393 00394 struct Error__Global_statistics_are_only_defined_for_AccumulatorChainArray; 00395 00396 /** \brief Specifies index of labels in CoupledHandle. 00397 00398 LabelArg<INDEX> tells the acc::AccumulatorChainArray which index of the Handle contains the labels. (Note that coordinates are always index 0) 00399 */ 00400 template <int INDEX> 00401 class LabelArg 00402 { 00403 public: 00404 typedef Select<> Dependencies; 00405 00406 static std::string const & name() 00407 { 00408 static const std::string n = std::string("LabelArg<") + asString(INDEX) + "> (internal)"; 00409 return n; 00410 } 00411 00412 template <class T, class BASE> 00413 struct Impl 00414 : public BASE 00415 { 00416 typedef LabelArgTag Tag; 00417 typedef void value_type; 00418 typedef void result_type; 00419 00420 static const int value = INDEX; 00421 static const unsigned int workInPass = 0; 00422 }; 00423 }; 00424 00425 template <int INDEX> 00426 class CoordArg 00427 { 00428 public: 00429 typedef Select<> Dependencies; 00430 00431 static std::string const & name() 00432 { 00433 static const std::string n = std::string("CoordArg<") + asString(INDEX) + "> (internal)"; 00434 return n; 00435 } 00436 00437 template <class T, class BASE> 00438 struct Impl 00439 : public BASE 00440 { 00441 typedef CoordArgTag Tag; 00442 typedef void value_type; 00443 typedef void result_type; 00444 00445 static const int value = INDEX; 00446 static const unsigned int workInPass = 0; 00447 }; 00448 }; 00449 00450 template <class T, class TAG, class NEXT=AccumulatorEnd> 00451 struct AccumulatorBase; 00452 00453 template <class Tag, class A> 00454 struct LookupTag; 00455 00456 template <class Tag, class A, class TargetTag=typename A::Tag> 00457 struct LookupDependency; 00458 00459 #ifndef _MSC_VER // compiler bug? (causes 'ambiguous overload error') 00460 00461 template <class TAG, class A> 00462 typename LookupTag<TAG, A>::reference 00463 getAccumulator(A & a); 00464 00465 template <class TAG, class A> 00466 typename LookupDependency<TAG, A>::result_type 00467 getDependency(A const & a); 00468 00469 #endif 00470 00471 namespace detail { 00472 00473 /****************************************************************************/ 00474 /* */ 00475 /* internal tag handling meta-functions */ 00476 /* */ 00477 /****************************************************************************/ 00478 00479 // we must make sure that Arg<INDEX> tags are at the end of the chain because 00480 // all other tags potentially depend on them 00481 template <class T> 00482 struct PushArgTagToTail 00483 { 00484 typedef T type; 00485 }; 00486 00487 #define VIGRA_PUSHARGTAG(TAG) \ 00488 template <int INDEX, class TAIL> \ 00489 struct PushArgTagToTail<TypeList<TAG<INDEX>, TAIL> > \ 00490 { \ 00491 typedef typename Push<TAIL, TypeList<TAG<INDEX> > >::type type; \ 00492 }; 00493 00494 VIGRA_PUSHARGTAG(DataArg) 00495 VIGRA_PUSHARGTAG(WeightArg) 00496 VIGRA_PUSHARGTAG(CoordArg) 00497 VIGRA_PUSHARGTAG(LabelArg) 00498 00499 #undef VIGRA_PUSHARGTAG 00500 00501 // Insert the dependencies of the selected functors into the TypeList and sort 00502 // the list such that dependencies come after the functors using them. Make sure 00503 // that each functor is contained only once. 00504 template <class T> 00505 struct AddDependencies; 00506 00507 template <class HEAD, class TAIL> 00508 struct AddDependencies<TypeList<HEAD, TAIL> > 00509 { 00510 typedef typename AddDependencies<TAIL>::type TailWithDependencies; 00511 typedef typename StandardizeDependencies<HEAD>::type HeadDependencies; 00512 typedef typename AddDependencies<HeadDependencies>::type TransitiveHeadDependencies; 00513 typedef TypeList<HEAD, TransitiveHeadDependencies> HeadWithDependencies; 00514 typedef typename PushUnique<HeadWithDependencies, TailWithDependencies>::type UnsortedDependencies; 00515 typedef typename PushArgTagToTail<UnsortedDependencies>::type type; 00516 }; 00517 00518 template <> 00519 struct AddDependencies<void> 00520 { 00521 typedef void type; 00522 }; 00523 00524 // Helper class to activate dependencies at runtime (i.e. when activate<Tag>(accu) is called, 00525 // activate() must also be called for Tag's dependencies). 00526 template <class Dependencies> 00527 struct ActivateDependencies; 00528 00529 template <class HEAD, class TAIL> 00530 struct ActivateDependencies<TypeList<HEAD, TAIL> > 00531 { 00532 template <class Chain, class ActiveFlags> 00533 static void exec(ActiveFlags & flags) 00534 { 00535 LookupTag<HEAD, Chain>::type::activateImpl(flags); 00536 ActivateDependencies<TAIL>::template exec<Chain>(flags); 00537 } 00538 00539 template <class Chain, class ActiveFlags, class GlobalFlags> 00540 static void exec(ActiveFlags & flags, GlobalFlags & gflags) 00541 { 00542 LookupTag<HEAD, Chain>::type::template activateImpl<Chain>(flags, gflags); 00543 ActivateDependencies<TAIL>::template exec<Chain>(flags, gflags); 00544 } 00545 }; 00546 00547 template <class HEAD, class TAIL> 00548 struct ActivateDependencies<TypeList<Global<HEAD>, TAIL> > 00549 { 00550 template <class Chain, class ActiveFlags, class GlobalFlags> 00551 static void exec(ActiveFlags & flags, GlobalFlags & gflags) 00552 { 00553 LookupTag<Global<HEAD>, Chain>::type::activateImpl(gflags); 00554 ActivateDependencies<TAIL>::template exec<Chain>(flags, gflags); 00555 } 00556 }; 00557 00558 template <> 00559 struct ActivateDependencies<void> 00560 { 00561 template <class Chain, class ActiveFlags> 00562 static void exec(ActiveFlags &) 00563 {} 00564 00565 template <class Chain, class ActiveFlags, class GlobalFlags> 00566 static void exec(ActiveFlags &, GlobalFlags &) 00567 {} 00568 }; 00569 00570 template <class List> 00571 struct SeparateGlobalAndRegionTags; 00572 00573 template <class HEAD, class TAIL> 00574 struct SeparateGlobalAndRegionTags<TypeList<HEAD, TAIL> > 00575 { 00576 typedef SeparateGlobalAndRegionTags<TAIL> Inner; 00577 typedef TypeList<HEAD, typename Inner::RegionTags> RegionTags; 00578 typedef typename Inner::GlobalTags GlobalTags; 00579 }; 00580 00581 template <class HEAD, class TAIL> 00582 struct SeparateGlobalAndRegionTags<TypeList<Global<HEAD>, TAIL> > 00583 { 00584 typedef SeparateGlobalAndRegionTags<TAIL> Inner; 00585 typedef typename Inner::RegionTags RegionTags; 00586 typedef TypeList<HEAD, typename Inner::GlobalTags> GlobalTags; 00587 }; 00588 00589 template <int INDEX, class TAIL> 00590 struct SeparateGlobalAndRegionTags<TypeList<DataArg<INDEX>, TAIL> > 00591 { 00592 typedef SeparateGlobalAndRegionTags<TAIL> Inner; 00593 typedef TypeList<DataArg<INDEX>, typename Inner::RegionTags> RegionTags; 00594 typedef TypeList<DataArg<INDEX>, typename Inner::GlobalTags> GlobalTags; 00595 }; 00596 00597 template <int INDEX, class TAIL> 00598 struct SeparateGlobalAndRegionTags<TypeList<LabelArg<INDEX>, TAIL> > 00599 { 00600 typedef SeparateGlobalAndRegionTags<TAIL> Inner; 00601 typedef typename Inner::RegionTags RegionTags; 00602 typedef TypeList<LabelArg<INDEX>, typename Inner::GlobalTags> GlobalTags; 00603 }; 00604 00605 template <int INDEX, class TAIL> 00606 struct SeparateGlobalAndRegionTags<TypeList<WeightArg<INDEX>, TAIL> > 00607 { 00608 typedef SeparateGlobalAndRegionTags<TAIL> Inner; 00609 typedef TypeList<WeightArg<INDEX>, typename Inner::RegionTags> RegionTags; 00610 typedef TypeList<WeightArg<INDEX>, typename Inner::GlobalTags> GlobalTags; 00611 }; 00612 00613 template <int INDEX, class TAIL> 00614 struct SeparateGlobalAndRegionTags<TypeList<CoordArg<INDEX>, TAIL> > 00615 { 00616 typedef SeparateGlobalAndRegionTags<TAIL> Inner; 00617 typedef TypeList<CoordArg<INDEX>, typename Inner::RegionTags> RegionTags; 00618 typedef TypeList<CoordArg<INDEX>, typename Inner::GlobalTags> GlobalTags; 00619 }; 00620 00621 template <> 00622 struct SeparateGlobalAndRegionTags<void> 00623 { 00624 typedef void RegionTags; 00625 typedef void GlobalTags; 00626 }; 00627 00628 /****************************************************************************/ 00629 /* */ 00630 /* helper classes to handle tags at runtime via strings */ 00631 /* */ 00632 /****************************************************************************/ 00633 00634 template <class Accumulators> 00635 struct CollectAccumulatorNames; 00636 00637 template <class HEAD, class TAIL> 00638 struct CollectAccumulatorNames<TypeList<HEAD, TAIL> > 00639 { 00640 template <class BackInsertable> 00641 static void exec(BackInsertable & a, bool skipInternals=true) 00642 { 00643 if(!skipInternals || HEAD::name().find("internal") == std::string::npos) 00644 a.push_back(HEAD::name()); 00645 CollectAccumulatorNames<TAIL>::exec(a, skipInternals); 00646 } 00647 }; 00648 00649 template <> 00650 struct CollectAccumulatorNames<void> 00651 { 00652 template <class BackInsertable> 00653 static void exec(BackInsertable & a, bool skipInternals=true) 00654 {} 00655 }; 00656 00657 template <class T> 00658 struct ApplyVisitorToTag; 00659 00660 template <class HEAD, class TAIL> 00661 struct ApplyVisitorToTag<TypeList<HEAD, TAIL> > 00662 { 00663 template <class Accu, class Visitor> 00664 static bool exec(Accu & a, std::string const & tag, Visitor const & v) 00665 { 00666 static const std::string name = normalizeString(HEAD::name()); 00667 if(name == tag) 00668 { 00669 v.template exec<HEAD>(a); 00670 return true; 00671 } 00672 else 00673 { 00674 return ApplyVisitorToTag<TAIL>::exec(a, tag, v); 00675 } 00676 } 00677 }; 00678 00679 template <> 00680 struct ApplyVisitorToTag<void> 00681 { 00682 template <class Accu, class Visitor> 00683 static bool exec(Accu & a, std::string const & tag, Visitor const & v) 00684 { 00685 return false; 00686 } 00687 }; 00688 00689 struct ActivateTag_Visitor 00690 { 00691 template <class TAG, class Accu> 00692 void exec(Accu & a) const 00693 { 00694 a.template activate<TAG>(); 00695 } 00696 }; 00697 00698 struct TagIsActive_Visitor 00699 { 00700 mutable bool result; 00701 00702 template <class TAG, class Accu> 00703 void exec(Accu & a) const 00704 { 00705 result = a.template isActive<TAG>(); 00706 } 00707 }; 00708 00709 /****************************************************************************/ 00710 /* */ 00711 /* histogram initialization functors */ 00712 /* */ 00713 /****************************************************************************/ 00714 00715 template <class TAG> 00716 struct SetHistogramBincount 00717 { 00718 template <class Accu> 00719 static void exec(Accu & a, HistogramOptions const & options) 00720 {} 00721 }; 00722 00723 template <template <int> class Histogram> 00724 struct SetHistogramBincount<Histogram<0> > 00725 { 00726 template <class Accu> 00727 static void exec(Accu & a, HistogramOptions const & options) 00728 { 00729 a.setBinCount(options.binCount); 00730 } 00731 }; 00732 00733 template <class TAG> 00734 struct ApplyHistogramOptions 00735 { 00736 template <class Accu> 00737 static void exec(Accu & a, HistogramOptions const & options) 00738 {} 00739 }; 00740 00741 template <class TAG> 00742 struct ApplyHistogramOptions<StandardQuantiles<TAG> > 00743 { 00744 template <class Accu> 00745 static void exec(Accu & a, HistogramOptions const & options) 00746 {} 00747 }; 00748 00749 template <class TAG, template <class> class MODIFIER> 00750 struct ApplyHistogramOptions<MODIFIER<TAG> > 00751 : public ApplyHistogramOptions<TAG> 00752 {}; 00753 00754 template <> 00755 struct ApplyHistogramOptions<IntegerHistogram<0> > 00756 { 00757 template <class Accu> 00758 static void exec(Accu & a, HistogramOptions const & options) 00759 { 00760 SetHistogramBincount<IntegerHistogram<0> >::exec(a, options); 00761 } 00762 }; 00763 00764 template <int BinCount> 00765 struct ApplyHistogramOptions<UserRangeHistogram<BinCount> > 00766 { 00767 template <class Accu> 00768 static void exec(Accu & a, HistogramOptions const & options) 00769 { 00770 SetHistogramBincount<UserRangeHistogram<BinCount> >::exec(a, options); 00771 if(a.scale_ == 0.0 && options.validMinMax()) 00772 a.setMinMax(options.minimum, options.maximum); 00773 } 00774 }; 00775 00776 template <int BinCount> 00777 struct ApplyHistogramOptions<AutoRangeHistogram<BinCount> > 00778 { 00779 template <class Accu> 00780 static void exec(Accu & a, HistogramOptions const & options) 00781 { 00782 SetHistogramBincount<AutoRangeHistogram<BinCount> >::exec(a, options); 00783 if(a.scale_ == 0.0 && options.validMinMax()) 00784 a.setMinMax(options.minimum, options.maximum); 00785 } 00786 }; 00787 00788 template <int BinCount> 00789 struct ApplyHistogramOptions<GlobalRangeHistogram<BinCount> > 00790 { 00791 template <class Accu> 00792 static void exec(Accu & a, HistogramOptions const & options) 00793 { 00794 SetHistogramBincount<GlobalRangeHistogram<BinCount> >::exec(a, options); 00795 if(a.scale_ == 0.0) 00796 { 00797 if(options.validMinMax()) 00798 a.setMinMax(options.minimum, options.maximum); 00799 else 00800 a.setRegionAutoInit(options.local_auto_init); 00801 } 00802 } 00803 }; 00804 00805 /****************************************************************************/ 00806 /* */ 00807 /* internal accumulator chain classes */ 00808 /* */ 00809 /****************************************************************************/ 00810 00811 // AccumulatorEndImpl has the following functionalities: 00812 // * marks end of accumulator chain by the AccumulatorEnd tag 00813 // * provides empty implementation of standard accumulator functions 00814 // * provides active_accumulators_ flags for run-time activation of dynamic accumulators 00815 // * provides is_dirty_ flags for caching accumulators 00816 // * hold the GlobalAccumulatorHandle for global accumulator lookup from region accumulators 00817 template <unsigned LEVEL, class GlobalAccumulatorHandle> 00818 struct AccumulatorEndImpl 00819 { 00820 typedef typename GlobalAccumulatorHandle::type GlobalAccumulatorType; 00821 00822 typedef AccumulatorEnd Tag; 00823 typedef void value_type; 00824 typedef bool result_type; 00825 typedef BitArray<LEVEL> AccumulatorFlags; 00826 00827 static const unsigned int workInPass = 0; 00828 static const int index = -1; 00829 static const unsigned level = LEVEL; 00830 00831 AccumulatorFlags active_accumulators_; 00832 mutable AccumulatorFlags is_dirty_; 00833 GlobalAccumulatorHandle globalAccumulator_; 00834 00835 template <class GlobalAccumulator> 00836 void setGlobalAccumulator(GlobalAccumulator const * a) 00837 { 00838 globalAccumulator_.pointer_ = a; 00839 } 00840 00841 static std::string name() 00842 { 00843 return "AccumulatorEnd (internal)"; 00844 } 00845 00846 bool operator()() const { return false; } 00847 bool get() const { return false; } 00848 00849 template <unsigned, class U> 00850 void pass(U const &) 00851 {} 00852 00853 template <unsigned, class U> 00854 void pass(U const &, double) 00855 {} 00856 00857 template <class U> 00858 void merge(U const &) 00859 {} 00860 00861 template <class U> 00862 void resize(U const &) 00863 {} 00864 00865 void activate() 00866 {} 00867 00868 bool isActive() const 00869 { 00870 return false; 00871 } 00872 00873 template <class Flags> 00874 static void activateImpl(Flags &) 00875 {} 00876 00877 template <class Accu, class Flags1, class Flags2> 00878 static void activateImpl(Flags1 &, Flags2 &) 00879 {} 00880 00881 template <class Flags> 00882 static bool isActiveImpl(Flags const &) 00883 { 00884 return true; 00885 } 00886 00887 void applyHistogramOptions(HistogramOptions const &) 00888 {} 00889 00890 static unsigned int passesRequired() 00891 { 00892 return 0; 00893 } 00894 00895 static unsigned int passesRequired(AccumulatorFlags const &) 00896 { 00897 return 0; 00898 } 00899 00900 void reset() 00901 { 00902 active_accumulators_.clear(); 00903 is_dirty_.clear(); 00904 } 00905 00906 template <int which> 00907 void setDirtyImpl() const 00908 { 00909 is_dirty_.template set<which>(); 00910 } 00911 00912 template <int which> 00913 void setCleanImpl() const 00914 { 00915 is_dirty_.template reset<which>(); 00916 } 00917 00918 template <int which> 00919 bool isDirtyImpl() const 00920 { 00921 return is_dirty_.template test<which>(); 00922 } 00923 }; 00924 00925 // DecoratorImpl implement the functionality of Decorator below 00926 template <class A, unsigned CurrentPass, bool allowRuntimeActivation, unsigned WorkPass=A::workInPass> 00927 struct DecoratorImpl 00928 { 00929 template <class T> 00930 static void exec(A & a, T const & t) 00931 {} 00932 00933 template <class T> 00934 static void exec(A & a, T const & t, double weight) 00935 {} 00936 }; 00937 00938 template <class A, unsigned CurrentPass> 00939 struct DecoratorImpl<A, CurrentPass, false, CurrentPass> 00940 { 00941 template <class T> 00942 static void exec(A & a, T const & t) 00943 { 00944 a.update(t); 00945 } 00946 00947 template <class T> 00948 static void exec(A & a, T const & t, double weight) 00949 { 00950 a.update(t, weight); 00951 } 00952 00953 static typename A::result_type get(A const & a) 00954 { 00955 return a(); 00956 } 00957 00958 static void merge(A & a, A const & o) 00959 { 00960 a += o; 00961 } 00962 00963 template <class T> 00964 static void resize(A & a, T const & t) 00965 { 00966 a.reshape(t); 00967 } 00968 00969 static void applyHistogramOptions(A & a, HistogramOptions const & options) 00970 { 00971 ApplyHistogramOptions<typename A::Tag>::exec(a, options); 00972 } 00973 00974 static unsigned int passesRequired() 00975 { 00976 static const unsigned int A_workInPass = A::workInPass; 00977 return std::max(A_workInPass, A::InternalBaseType::passesRequired()); 00978 } 00979 }; 00980 00981 template <class A, unsigned CurrentPass> 00982 struct DecoratorImpl<A, CurrentPass, true, CurrentPass> 00983 { 00984 static bool isActive(A const & a) 00985 { 00986 return A::isActiveImpl(getAccumulator<AccumulatorEnd>(a).active_accumulators_); 00987 } 00988 00989 template <class T> 00990 static void exec(A & a, T const & t) 00991 { 00992 if(isActive(a)) 00993 a.update(t); 00994 } 00995 00996 template <class T> 00997 static void exec(A & a, T const & t, double weight) 00998 { 00999 if(isActive(a)) 01000 a.update(t, weight); 01001 } 01002 01003 static typename A::result_type get(A const & a) 01004 { 01005 static const std::string message = std::string("get(accumulator): attempt to access inactive statistic '") + 01006 typeid(typename A::Tag).name() + "'."; 01007 vigra_precondition(isActive(a), message); 01008 return a(); 01009 } 01010 01011 static void merge(A & a, A const & o) 01012 { 01013 if(isActive(a)) 01014 a += o; 01015 } 01016 01017 template <class T> 01018 static void resize(A & a, T const & t) 01019 { 01020 if(isActive(a)) 01021 a.reshape(t); 01022 } 01023 01024 static void applyHistogramOptions(A & a, HistogramOptions const & options) 01025 { 01026 if(isActive(a)) 01027 ApplyHistogramOptions<typename A::Tag>::exec(a, options); 01028 } 01029 01030 template <class ActiveFlags> 01031 static unsigned int passesRequired(ActiveFlags const & flags) 01032 { 01033 static const unsigned int A_workInPass = A::workInPass; 01034 return A::isActiveImpl(flags) 01035 ? std::max(A_workInPass, A::InternalBaseType::passesRequired(flags)) 01036 : A::InternalBaseType::passesRequired(flags); 01037 } 01038 }; 01039 01040 // Generic reshape function (expands to a no-op when T has fixed shape, and to 01041 // the appropriate specialized call otherwise). Shape is an instance of MultiArrayShape<N>::type. 01042 template <class T, class Shape> 01043 void reshapeImpl(T &, Shape const &) 01044 {} 01045 01046 template <class T, class Shape, class Initial> 01047 void reshapeImpl(T &, Shape const &, Initial const & = T()) 01048 {} 01049 01050 template <unsigned int N, class T, class Alloc, class Shape> 01051 void reshapeImpl(MultiArray<N, T, Alloc> & a, Shape const & s, T const & initial = T()) 01052 { 01053 MultiArray<N, T, Alloc>(s, initial).swap(a); 01054 } 01055 01056 template <class T, class Alloc, class Shape> 01057 void reshapeImpl(Matrix<T, Alloc> & a, Shape const & s, T const & initial = T()) 01058 { 01059 Matrix<T, Alloc>(s, initial).swap(a); 01060 } 01061 01062 // generic functions to create suitable shape objects from various input data types 01063 template <unsigned int N, class T, class Stride> 01064 inline typename MultiArrayShape<N>::type 01065 shapeOf(MultiArrayView<N, T, Stride> const & a) 01066 { 01067 return a.shape(); 01068 } 01069 01070 template <class T, int N> 01071 inline Shape1 01072 shapeOf(TinyVector<T, N> const &) 01073 { 01074 return Shape1(N); 01075 } 01076 01077 template <class T, class NEXT> 01078 inline CoupledHandle<T, NEXT> const & 01079 shapeOf(CoupledHandle<T, NEXT> const & t) 01080 { 01081 return t; 01082 } 01083 01084 #define VIGRA_SHAPE_OF(type) \ 01085 inline Shape1 \ 01086 shapeOf(type) \ 01087 { \ 01088 return Shape1(1); \ 01089 } 01090 01091 VIGRA_SHAPE_OF(unsigned char) 01092 VIGRA_SHAPE_OF(signed char) 01093 VIGRA_SHAPE_OF(unsigned short) 01094 VIGRA_SHAPE_OF(short) 01095 VIGRA_SHAPE_OF(unsigned int) 01096 VIGRA_SHAPE_OF(int) 01097 VIGRA_SHAPE_OF(unsigned long) 01098 VIGRA_SHAPE_OF(long) 01099 VIGRA_SHAPE_OF(unsigned long long) 01100 VIGRA_SHAPE_OF(long long) 01101 VIGRA_SHAPE_OF(float) 01102 VIGRA_SHAPE_OF(double) 01103 VIGRA_SHAPE_OF(long double) 01104 01105 #undef VIGRA_SHAPE_OF 01106 01107 // LabelDispatch is only used in AccumulatorChainArrays and has the following functionalities: 01108 // * hold an accumulator chain for global statistics 01109 // * hold an array of accumulator chains (one per region) for region statistics 01110 // * forward data to the appropriate chains 01111 // * allocate the region array with appropriate size 01112 // * store and forward activation requests 01113 // * compute required number of passes as maximum from global and region accumulators 01114 template <class T, class GlobalAccumulators, class RegionAccumulators> 01115 struct LabelDispatch 01116 { 01117 typedef LabelDispatchTag Tag; 01118 typedef GlobalAccumulators GlobalAccumulatorChain; 01119 typedef RegionAccumulators RegionAccumulatorChain; 01120 typedef typename LookupTag<AccumulatorEnd, RegionAccumulatorChain>::type::AccumulatorFlags ActiveFlagsType; 01121 typedef ArrayVector<RegionAccumulatorChain> RegionAccumulatorArray; 01122 01123 typedef LabelDispatch type; 01124 typedef LabelDispatch & reference; 01125 typedef LabelDispatch const & const_reference; 01126 typedef GlobalAccumulatorChain InternalBaseType; 01127 01128 typedef T const & argument_type; 01129 typedef argument_type first_argument_type; 01130 typedef double second_argument_type; 01131 typedef RegionAccumulatorChain & result_type; 01132 01133 static const int index = GlobalAccumulatorChain::index + 1; 01134 01135 GlobalAccumulatorChain next_; 01136 RegionAccumulatorArray regions_; 01137 HistogramOptions region_histogram_options_; 01138 MultiArrayIndex ignore_label_; 01139 ActiveFlagsType active_region_accumulators_; 01140 01141 template <class IndexDefinition, class TagFound=typename IndexDefinition::Tag> 01142 struct LabelIndexSelector 01143 { 01144 static const int value = 2; // default: CoupledHandle holds labels at index 2 01145 01146 template <class U, class NEXT> 01147 static MultiArrayIndex exec(CoupledHandle<U, NEXT> const & t) 01148 { 01149 return (MultiArrayIndex)get<value>(t); 01150 } 01151 }; 01152 01153 template <class IndexDefinition> 01154 struct LabelIndexSelector<IndexDefinition, LabelArgTag> 01155 { 01156 static const int value = IndexDefinition::value; 01157 01158 template <class U, class NEXT> 01159 static MultiArrayIndex exec(CoupledHandle<U, NEXT> const & t) 01160 { 01161 return (MultiArrayIndex)get<value>(t); 01162 } 01163 }; 01164 01165 template <class TAG> 01166 struct ActivateImpl 01167 { 01168 typedef typename LookupTag<TAG, type>::type TargetAccumulator; 01169 01170 static void activate(GlobalAccumulatorChain & globals, RegionAccumulatorArray & regions, 01171 ActiveFlagsType & flags) 01172 { 01173 TargetAccumulator::template activateImpl<LabelDispatch>( 01174 flags, getAccumulator<AccumulatorEnd>(globals).active_accumulators_); 01175 for(unsigned int k=0; k<regions.size(); ++k) 01176 getAccumulator<AccumulatorEnd>(regions[k]).active_accumulators_ = flags; 01177 } 01178 01179 static bool isActive(GlobalAccumulatorChain const &, ActiveFlagsType const & flags) 01180 { 01181 return TargetAccumulator::isActiveImpl(flags); 01182 } 01183 }; 01184 01185 template <class TAG> 01186 struct ActivateImpl<Global<TAG> > 01187 { 01188 static void activate(GlobalAccumulatorChain & globals, RegionAccumulatorArray &, ActiveFlagsType &) 01189 { 01190 LookupTag<TAG, GlobalAccumulatorChain>::type::activateImpl(getAccumulator<AccumulatorEnd>(globals).active_accumulators_); 01191 } 01192 01193 static bool isActive(GlobalAccumulatorChain const & globals, ActiveFlagsType const &) 01194 { 01195 return LookupTag<TAG, GlobalAccumulatorChain>::type::isActiveImpl(getAccumulator<AccumulatorEnd>(globals).active_accumulators_); 01196 } 01197 }; 01198 01199 template <int INDEX> 01200 struct ActivateImpl<LabelArg<INDEX> > 01201 { 01202 static void activate(GlobalAccumulatorChain &, RegionAccumulatorArray &, ActiveFlagsType &) 01203 {} 01204 01205 static bool isActive(GlobalAccumulatorChain const & globals, ActiveFlagsType const &) 01206 { 01207 return getAccumulator<LabelArg<INDEX> >(globals).isActive(); 01208 } 01209 }; 01210 01211 typedef typename LookupTag<LabelArgTag, GlobalAccumulatorChain>::type FindLabelIndex; 01212 01213 LabelDispatch() 01214 : next_(), 01215 regions_(), 01216 region_histogram_options_(), 01217 ignore_label_(-1), 01218 active_region_accumulators_() 01219 {} 01220 01221 LabelDispatch(LabelDispatch const & o) 01222 : next_(o.next_), 01223 regions_(o.regions_), 01224 region_histogram_options_(o.region_histogram_options_), 01225 ignore_label_(o.ignore_label_), 01226 active_region_accumulators_(o.active_region_accumulators_) 01227 { 01228 for(unsigned int k=0; k<regions_.size(); ++k) 01229 { 01230 getAccumulator<AccumulatorEnd>(regions_[k]).setGlobalAccumulator(&next_); 01231 } 01232 } 01233 01234 MultiArrayIndex maxRegionLabel() const 01235 { 01236 return (MultiArrayIndex)regions_.size() - 1; 01237 } 01238 01239 void setMaxRegionLabel(unsigned maxlabel) 01240 { 01241 if(maxRegionLabel() == (MultiArrayIndex)maxlabel) 01242 return; 01243 unsigned int oldSize = regions_.size(); 01244 regions_.resize(maxlabel + 1); 01245 for(unsigned int k=oldSize; k<regions_.size(); ++k) 01246 { 01247 getAccumulator<AccumulatorEnd>(regions_[k]).setGlobalAccumulator(&next_); 01248 getAccumulator<AccumulatorEnd>(regions_[k]).active_accumulators_ = active_region_accumulators_; 01249 regions_[k].applyHistogramOptions(region_histogram_options_); 01250 } 01251 } 01252 01253 void ignoreLabel(MultiArrayIndex l) 01254 { 01255 ignore_label_ = l; 01256 } 01257 01258 void applyHistogramOptions(HistogramOptions const & options) 01259 { 01260 applyHistogramOptions(options, options); 01261 } 01262 01263 void applyHistogramOptions(HistogramOptions const & regionoptions, HistogramOptions const & globaloptions) 01264 { 01265 region_histogram_options_ = regionoptions; 01266 next_.applyHistogramOptions(globaloptions); 01267 } 01268 01269 template <class U> 01270 void resize(U const & t) 01271 { 01272 if(regions_.size() == 0) 01273 { 01274 static const int labelIndex = LabelIndexSelector<FindLabelIndex>::value; 01275 typedef typename CoupledHandleCast<labelIndex, T>::type LabelHandle; 01276 typedef typename LabelHandle::value_type LabelType; 01277 typedef MultiArrayView<LabelHandle::dimensions, LabelType, StridedArrayTag> LabelArray; 01278 LabelArray labelArray(t.shape(), cast<labelIndex>(t).strides(), const_cast<LabelType *>(cast<labelIndex>(t).ptr())); 01279 01280 LabelType minimum, maximum; 01281 labelArray.minmax(&minimum, &maximum); 01282 setMaxRegionLabel(maximum); 01283 } 01284 next_.resize(t); 01285 // FIXME: only call resize when label k actually exists? 01286 for(unsigned int k=0; k<regions_.size(); ++k) 01287 regions_[k].resize(t); 01288 } 01289 01290 template <unsigned N> 01291 void pass(T const & t) 01292 { 01293 if(LabelIndexSelector<FindLabelIndex>::exec(t) != ignore_label_) 01294 { 01295 next_.template pass<N>(t); 01296 regions_[LabelIndexSelector<FindLabelIndex>::exec(t)].template pass<N>(t); 01297 } 01298 } 01299 01300 template <unsigned N> 01301 void pass(T const & t, double weight) 01302 { 01303 if(LabelIndexSelector<FindLabelIndex>::exec(t) != ignore_label_) 01304 { 01305 next_.template pass<N>(t, weight); 01306 regions_[LabelIndexSelector<FindLabelIndex>::exec(t)].template pass<N>(t, weight); 01307 } 01308 } 01309 01310 static unsigned int passesRequired() 01311 { 01312 return std::max(GlobalAccumulatorChain::passesRequired(), RegionAccumulatorChain::passesRequired()); 01313 } 01314 01315 unsigned int passesRequiredDynamic() const 01316 { 01317 return std::max(GlobalAccumulatorChain::passesRequired(getAccumulator<AccumulatorEnd>(next_).active_accumulators_), 01318 RegionAccumulatorChain::passesRequired(active_region_accumulators_)); 01319 } 01320 01321 void reset() 01322 { 01323 next_.reset(); 01324 01325 active_region_accumulators_.clear(); 01326 RegionAccumulatorArray().swap(regions_); 01327 // FIXME: or is it better to just reset the region accumulators? 01328 // for(unsigned int k=0; k<regions_.size(); ++k) 01329 // regions_[k].reset(); 01330 } 01331 01332 template <class TAG> 01333 void activate() 01334 { 01335 ActivateImpl<TAG>::activate(next_, regions_, active_region_accumulators_); 01336 } 01337 01338 void activateAll() 01339 { 01340 getAccumulator<AccumulatorEnd>(next_).active_accumulators_.set(); 01341 active_region_accumulators_.set(); 01342 for(unsigned int k=0; k<regions_.size(); ++k) 01343 getAccumulator<AccumulatorEnd>(regions_[k]).active_accumulators_.set(); 01344 } 01345 01346 template <class TAG> 01347 bool isActive() const 01348 { 01349 return ActivateImpl<TAG>::isActive(next_, active_region_accumulators_); 01350 } 01351 01352 void merge(LabelDispatch const & o) 01353 { 01354 for(unsigned int k=0; k<regions_.size(); ++k) 01355 regions_[k].merge(o.regions_[k]); 01356 next_.merge(o.next_); 01357 } 01358 01359 void merge(unsigned i, unsigned j) 01360 { 01361 regions_[i].merge(regions_[j]); 01362 regions_[j].reset(); 01363 getAccumulator<AccumulatorEnd>(regions_[j]).active_accumulators_ = active_region_accumulators_; 01364 } 01365 01366 template <class ArrayLike> 01367 void merge(LabelDispatch const & o, ArrayLike const & labelMapping) 01368 { 01369 MultiArrayIndex newMaxLabel = std::max<MultiArrayIndex>(maxRegionLabel(), *argMax(labelMapping.begin(), labelMapping.end())); 01370 setMaxRegionLabel(newMaxLabel); 01371 for(unsigned int k=0; k<labelMapping.size(); ++k) 01372 regions_[labelMapping[k]].merge(o.regions_[k]); 01373 next_.merge(o.next_); 01374 } 01375 }; 01376 01377 template <class TargetTag, class TagList> 01378 struct FindNextTag; 01379 01380 template <class TargetTag, class HEAD, class TAIL> 01381 struct FindNextTag<TargetTag, TypeList<HEAD, TAIL> > 01382 { 01383 typedef typename FindNextTag<TargetTag, TAIL>::type type; 01384 }; 01385 01386 template <class TargetTag, class TAIL> 01387 struct FindNextTag<TargetTag, TypeList<TargetTag, TAIL> > 01388 { 01389 typedef typename TAIL::Head type; 01390 }; 01391 01392 template <class TargetTag> 01393 struct FindNextTag<TargetTag, TypeList<TargetTag, void> > 01394 { 01395 typedef void type; 01396 }; 01397 01398 template <class TargetTag> 01399 struct FindNextTag<TargetTag, void> 01400 { 01401 typedef void type; 01402 }; 01403 01404 // AccumulatorFactory creates the decorator hierarchy for the given TAG and configuration CONFIG 01405 template <class TAG, class CONFIG, unsigned LEVEL=0> 01406 struct AccumulatorFactory 01407 { 01408 typedef typename FindNextTag<TAG, typename CONFIG::TagList>::type NextTag; 01409 typedef typename AccumulatorFactory<NextTag, CONFIG, LEVEL+1>::type NextType; 01410 typedef typename CONFIG::InputType InputType; 01411 01412 template <class T> 01413 struct ConfigureTag 01414 { 01415 typedef TAG type; 01416 }; 01417 01418 // When InputType is a CoupledHandle, some tags need to be wrapped into 01419 // DataFromHandle<> and/or Weighted<> modifiers. The following code does 01420 // this when appropriate. 01421 template <class T, class NEXT> 01422 struct ConfigureTag<CoupledHandle<T, NEXT> > 01423 { 01424 typedef typename StandardizeTag<DataFromHandle<TAG> >::type WrappedTag; 01425 typedef typename IfBool<(!HasModifierPriority<WrappedTag, WeightingPriority>::value && ShouldBeWeighted<WrappedTag>::value), 01426 Weighted<WrappedTag>, WrappedTag>::type type; 01427 }; 01428 01429 typedef typename ConfigureTag<InputType>::type UseTag; 01430 01431 // base class of the decorator hierarchy: default (possibly empty) 01432 // implementations of all members 01433 struct AccumulatorBase 01434 { 01435 typedef AccumulatorBase ThisType; 01436 typedef TAG Tag; 01437 typedef NextType InternalBaseType; 01438 typedef InputType input_type; 01439 typedef input_type const & argument_type; 01440 typedef argument_type first_argument_type; 01441 typedef double second_argument_type; 01442 typedef void result_type; 01443 01444 static const unsigned int workInPass = 1; 01445 static const int index = InternalBaseType::index + 1; 01446 01447 InternalBaseType next_; 01448 01449 static std::string name() 01450 { 01451 return TAG::name(); 01452 } 01453 01454 template <class ActiveFlags> 01455 static void activateImpl(ActiveFlags & flags) 01456 { 01457 flags.template set<index>(); 01458 typedef typename StandardizeDependencies<Tag>::type StdDeps; 01459 detail::ActivateDependencies<StdDeps>::template exec<ThisType>(flags); 01460 } 01461 01462 template <class Accu, class ActiveFlags, class GlobalFlags> 01463 static void activateImpl(ActiveFlags & flags, GlobalFlags & gflags) 01464 { 01465 flags.template set<index>(); 01466 typedef typename StandardizeDependencies<Tag>::type StdDeps; 01467 detail::ActivateDependencies<StdDeps>::template exec<Accu>(flags, gflags); 01468 } 01469 01470 template <class ActiveFlags> 01471 static bool isActiveImpl(ActiveFlags & flags) 01472 { 01473 return flags.template test<index>(); 01474 } 01475 01476 void setDirty() const 01477 { 01478 next_.template setDirtyImpl<index>(); 01479 } 01480 01481 template <int INDEX> 01482 void setDirtyImpl() const 01483 { 01484 next_.template setDirtyImpl<INDEX>(); 01485 } 01486 01487 void setClean() const 01488 { 01489 next_.template setCleanImpl<index>(); 01490 } 01491 01492 template <int INDEX> 01493 void setCleanImpl() const 01494 { 01495 next_.template setCleanImpl<INDEX>(); 01496 } 01497 01498 bool isDirty() const 01499 { 01500 return next_.template isDirtyImpl<index>(); 01501 } 01502 01503 template <int INDEX> 01504 bool isDirtyImpl() const 01505 { 01506 return next_.template isDirtyImpl<INDEX>(); 01507 } 01508 01509 void reset() 01510 {} 01511 01512 template <class Shape> 01513 void reshape(Shape const &) 01514 {} 01515 01516 void operator+=(AccumulatorBase const &) 01517 {} 01518 01519 template <class U> 01520 void update(U const &) 01521 {} 01522 01523 template <class U> 01524 void update(U const &, double) 01525 {} 01526 01527 template <class TargetTag> 01528 typename LookupDependency<TargetTag, ThisType>::result_type 01529 call_getDependency() const 01530 { 01531 return getDependency<TargetTag>(*this); 01532 } 01533 }; 01534 01535 // The middle class(es) of the decorator hierarchy implement the actual feature computation. 01536 typedef typename UseTag::template Impl<InputType, AccumulatorBase> AccumulatorImpl; 01537 01538 // outer class of the decorator hierarchy. It has the following functionalities 01539 // * ensure that only active accumulators are called in a dynamic accumulator chain 01540 // * ensure that each accumulator is only called in its desired pass as defined in A::workInPass 01541 // * determine how many passes through the data are required 01542 struct Accumulator 01543 : public AccumulatorImpl 01544 { 01545 typedef Accumulator type; 01546 typedef Accumulator & reference; 01547 typedef Accumulator const & const_reference; 01548 typedef AccumulatorImpl A; 01549 01550 static const unsigned int workInPass = A::workInPass; 01551 static const bool allowRuntimeActivation = CONFIG::allowRuntimeActivation; 01552 01553 template <class T> 01554 void resize(T const & t) 01555 { 01556 this->next_.resize(t); 01557 DecoratorImpl<Accumulator, workInPass, allowRuntimeActivation>::resize(*this, t); 01558 } 01559 01560 void reset() 01561 { 01562 this->next_.reset(); 01563 A::reset(); 01564 } 01565 01566 typename A::result_type get() const 01567 { 01568 return DecoratorImpl<A, workInPass, allowRuntimeActivation>::get(*this); 01569 } 01570 01571 template <unsigned N, class T> 01572 void pass(T const & t) 01573 { 01574 this->next_.template pass<N>(t); 01575 DecoratorImpl<Accumulator, N, allowRuntimeActivation>::exec(*this, t); 01576 } 01577 01578 template <unsigned N, class T> 01579 void pass(T const & t, double weight) 01580 { 01581 this->next_.template pass<N>(t, weight); 01582 DecoratorImpl<Accumulator, N, allowRuntimeActivation>::exec(*this, t, weight); 01583 } 01584 01585 void merge(Accumulator const & o) 01586 { 01587 DecoratorImpl<Accumulator, Accumulator::workInPass, allowRuntimeActivation>::merge(*this, o); 01588 this->next_.merge(o.next_); 01589 } 01590 01591 void applyHistogramOptions(HistogramOptions const & options) 01592 { 01593 DecoratorImpl<Accumulator, workInPass, allowRuntimeActivation>::applyHistogramOptions(*this, options); 01594 this->next_.applyHistogramOptions(options); 01595 } 01596 01597 static unsigned int passesRequired() 01598 { 01599 return DecoratorImpl<Accumulator, workInPass, allowRuntimeActivation>::passesRequired(); 01600 } 01601 01602 template <class ActiveFlags> 01603 static unsigned int passesRequired(ActiveFlags const & flags) 01604 { 01605 return DecoratorImpl<Accumulator, workInPass, allowRuntimeActivation>::passesRequired(flags); 01606 } 01607 }; 01608 01609 typedef Accumulator type; 01610 }; 01611 01612 template <class CONFIG, unsigned LEVEL> 01613 struct AccumulatorFactory<void, CONFIG, LEVEL> 01614 { 01615 typedef AccumulatorEndImpl<LEVEL, typename CONFIG::GlobalAccumulatorHandle> type; 01616 }; 01617 01618 struct InvalidGlobalAccumulatorHandle 01619 { 01620 typedef Error__Global_statistics_are_only_defined_for_AccumulatorChainArray type; 01621 01622 InvalidGlobalAccumulatorHandle() 01623 : pointer_(0) 01624 {} 01625 01626 type const * pointer_; 01627 }; 01628 01629 // helper classes to create an accumulator chain from a TypeList 01630 // if dynamic=true, a dynamic accumulator will be created 01631 // if dynamic=false, a plain accumulator will be created 01632 template <class T, class Selected, bool dynamic=false, class GlobalHandle=InvalidGlobalAccumulatorHandle> 01633 struct ConfigureAccumulatorChain 01634 #ifndef DOXYGEN 01635 : public ConfigureAccumulatorChain<T, typename AddDependencies<typename Selected::type>::type, dynamic> 01636 #endif 01637 {}; 01638 01639 template <class T, class HEAD, class TAIL, bool dynamic, class GlobalHandle> 01640 struct ConfigureAccumulatorChain<T, TypeList<HEAD, TAIL>, dynamic, GlobalHandle> 01641 { 01642 typedef TypeList<HEAD, TAIL> TagList; 01643 typedef T InputType; 01644 static const bool allowRuntimeActivation = dynamic; 01645 typedef GlobalHandle GlobalAccumulatorHandle; 01646 01647 typedef typename AccumulatorFactory<HEAD, ConfigureAccumulatorChain>::type type; 01648 }; 01649 01650 template <class T, class Selected, bool dynamic=false> 01651 struct ConfigureAccumulatorChainArray 01652 #ifndef DOXYGEN 01653 : public ConfigureAccumulatorChainArray<T, typename AddDependencies<typename Selected::type>::type, dynamic> 01654 #endif 01655 {}; 01656 01657 template <class T, class HEAD, class TAIL, bool dynamic> 01658 struct ConfigureAccumulatorChainArray<T, TypeList<HEAD, TAIL>, dynamic> 01659 { 01660 typedef TypeList<HEAD, TAIL> TagList; 01661 typedef SeparateGlobalAndRegionTags<TagList> TagSeparator; 01662 typedef typename TagSeparator::GlobalTags GlobalTags; 01663 typedef typename TagSeparator::RegionTags RegionTags; 01664 typedef typename ConfigureAccumulatorChain<T, GlobalTags, dynamic>::type GlobalAccumulatorChain; 01665 01666 struct GlobalAccumulatorHandle 01667 { 01668 typedef GlobalAccumulatorChain type; 01669 01670 GlobalAccumulatorHandle() 01671 : pointer_(0) 01672 {} 01673 01674 type const * pointer_; 01675 }; 01676 01677 typedef typename ConfigureAccumulatorChain<T, RegionTags, dynamic, GlobalAccumulatorHandle>::type RegionAccumulatorChain; 01678 01679 typedef LabelDispatch<T, GlobalAccumulatorChain, RegionAccumulatorChain> type; 01680 }; 01681 01682 } // namespace detail 01683 01684 /****************************************************************************/ 01685 /* */ 01686 /* accumulator chain */ 01687 /* */ 01688 /****************************************************************************/ 01689 01690 // Implement the high-level interface of an accumulator chain 01691 template <class T, class NEXT> 01692 class AccumulatorChainImpl 01693 { 01694 public: 01695 typedef NEXT InternalBaseType; 01696 typedef AccumulatorBegin Tag; 01697 typedef typename InternalBaseType::argument_type argument_type; 01698 typedef typename InternalBaseType::first_argument_type first_argument_type; 01699 typedef typename InternalBaseType::second_argument_type second_argument_type; 01700 typedef void value_type; 01701 typedef typename InternalBaseType::result_type result_type; 01702 01703 static const int staticSize = InternalBaseType::index; 01704 01705 InternalBaseType next_; 01706 01707 /** \brief Current pass of the accumulator chain. 01708 */ 01709 unsigned int current_pass_; 01710 01711 AccumulatorChainImpl() 01712 : current_pass_(0) 01713 {} 01714 01715 /** Set options for all histograms in the accumulator chain. See histogram accumulators for possible options. The function is ignored if there is no histogram in the accumulator chain. 01716 */ 01717 void setHistogramOptions(HistogramOptions const & options) 01718 { 01719 next_.applyHistogramOptions(options); 01720 } 01721 01722 01723 /** Set regional and global options for all histograms in the accumulator chain. 01724 */ 01725 void setHistogramOptions(HistogramOptions const & regionoptions, HistogramOptions const & globaloptions) 01726 { 01727 next_.applyHistogramOptions(regionoptions, globaloptions); 01728 } 01729 01730 /** Reset current_pass_ of the accumulator chain to 'reset_to_pass'. 01731 */ 01732 void reset(unsigned int reset_to_pass = 0) 01733 { 01734 current_pass_ = reset_to_pass; 01735 if(reset_to_pass == 0) 01736 next_.reset(); 01737 } 01738 01739 template <unsigned N> 01740 void update(T const & t) 01741 { 01742 if(current_pass_ == N) 01743 { 01744 next_.template pass<N>(t); 01745 } 01746 else if(current_pass_ < N) 01747 { 01748 current_pass_ = N; 01749 if(N == 1) 01750 next_.resize(detail::shapeOf(t)); 01751 next_.template pass<N>(t); 01752 } 01753 else 01754 { 01755 std::string message("AccumulatorChain::update(): cannot return to pass "); 01756 message << N << " after working on pass " << current_pass_ << "."; 01757 vigra_precondition(false, message); 01758 } 01759 } 01760 01761 template <unsigned N> 01762 void update(T const & t, double weight) 01763 { 01764 if(current_pass_ == N) 01765 { 01766 next_.template pass<N>(t, weight); 01767 } 01768 else if(current_pass_ < N) 01769 { 01770 current_pass_ = N; 01771 if(N == 1) 01772 next_.resize(detail::shapeOf(t)); 01773 next_.template pass<N>(t, weight); 01774 } 01775 else 01776 { 01777 std::string message("AccumulatorChain::update(): cannot return to pass "); 01778 message << N << " after working on pass " << current_pass_ << "."; 01779 vigra_precondition(false, message); 01780 } 01781 } 01782 01783 /** Equivalent to merge(o) . 01784 */ 01785 void operator+=(AccumulatorChainImpl const & o) 01786 { 01787 merge(o); 01788 } 01789 01790 /** Merge the accumulator chain with accumulator chain 'o'. This only works if all selected statistics in the accumulator chain support the '+=' operator. See the documentations of the particular statistics for support information. 01791 */ 01792 void merge(AccumulatorChainImpl const & o) 01793 { 01794 next_.merge(o.next_); 01795 } 01796 01797 result_type operator()() const 01798 { 01799 return next_.get(); 01800 } 01801 01802 void operator()(T const & t) 01803 { 01804 update<1>(t); 01805 } 01806 01807 void operator()(T const & t, double weight) 01808 { 01809 update<1>(t, weight); 01810 } 01811 01812 void updatePass2(T const & t) 01813 { 01814 update<2>(t); 01815 } 01816 01817 void updatePass2(T const & t, double weight) 01818 { 01819 update<2>(t, weight); 01820 } 01821 01822 /** Upate all accumulators in the accumulator chain that work in pass N with data t. Requirement: 0 < N < 6 and N >= current_pass_ . If N < current_pass_ call reset() first. 01823 */ 01824 void updatePassN(T const & t, unsigned int N) 01825 { 01826 switch (N) 01827 { 01828 case 1: update<1>(t); break; 01829 case 2: update<2>(t); break; 01830 case 3: update<3>(t); break; 01831 case 4: update<4>(t); break; 01832 case 5: update<5>(t); break; 01833 default: 01834 vigra_precondition(false, 01835 "AccumulatorChain::updatePassN(): 0 < N < 6 required."); 01836 } 01837 } 01838 01839 /** Upate all accumulators in the accumulator chain that work in pass N with data t and weight. Requirement: 0 < N < 6 and N >= current_pass_ . If N < current_pass_ call reset() first. 01840 */ 01841 void updatePassN(T const & t, double weight, unsigned int N) 01842 { 01843 switch (N) 01844 { 01845 case 1: update<1>(t, weight); break; 01846 case 2: update<2>(t, weight); break; 01847 case 3: update<3>(t, weight); break; 01848 case 4: update<4>(t, weight); break; 01849 case 5: update<5>(t, weight); break; 01850 default: 01851 vigra_precondition(false, 01852 "AccumulatorChain::updatePassN(): 0 < N < 6 required."); 01853 } 01854 } 01855 01856 /** Return the number of passes required to compute all statistics in the accumulator chain. 01857 */ 01858 unsigned int passesRequired() const 01859 { 01860 return InternalBaseType::passesRequired(); 01861 } 01862 }; 01863 01864 01865 01866 // Create an accumulator chain containing the Selected statistics and their dependencies. 01867 01868 /** \brief Create an accumulator chain containing the selected statistics and their dependencies. 01869 01870 AccumulatorChain is used to compute global statistics which have to be selected at compile time. 01871 01872 The template parameters are as follows: 01873 - T: The input type 01874 - either element type of the data(e.g. double, int, RGBValue, ...) 01875 - or type of CoupledHandle (for simultaneous access to coordinates and/or weights) 01876 - Selected: statistics to be computed and index specifier for the CoupledHandle, wrapped with Select 01877 01878 Usage: 01879 \code 01880 typedef double DataType; 01881 AccumulatorChain<DataType, Select<Variance, Mean, Minimum, ...> > accumulator; 01882 \endcode 01883 01884 Usage, using CoupledHandle: 01885 \code 01886 const int dim = 3; //dimension of MultiArray 01887 typedef double DataType; 01888 typedef double WeightType; 01889 typedef vigra::CoupledIteratorType<dim, DataType, WeightType>::HandleType Handle; 01890 AccumulatorChain<Handle, Select<DataArg<1>, WeightArg<2>, Mean,...> > a; 01891 \endcode 01892 01893 See \ref FeatureAccumulators for more information and examples of use. 01894 */ 01895 template <class T, class Selected, bool dynamic=false> 01896 class AccumulatorChain 01897 #ifndef DOXYGEN // hide AccumulatorChainImpl from documentation 01898 : public AccumulatorChainImpl<T, typename detail::ConfigureAccumulatorChain<T, Selected, dynamic>::type> 01899 #endif 01900 { 01901 public: 01902 // \brief TypeList of Tags in the accumulator chain (?). 01903 typedef typename detail::ConfigureAccumulatorChain<T, Selected, dynamic>::TagList AccumulatorTags; 01904 01905 /** Before having seen data (current_pass_==0), the shape of the data can be changed... (?) 01906 */ 01907 template <class U, int N> 01908 void reshape(TinyVector<U, N> const & s) 01909 { 01910 vigra_precondition(this->current_pass_ == 0, 01911 "AccumulatorChain::reshape(): cannot reshape after seeing data. Call AccumulatorChain::reset() first."); 01912 this->next_.resize(s); 01913 this->current_pass_ = 1; 01914 } 01915 01916 /** Return the names of all tags in the accumulator chain (selected statistics and their dependencies). 01917 */ 01918 static ArrayVector<std::string> const & tagNames() 01919 { 01920 static const ArrayVector<std::string> n = collectTagNames(); 01921 return n; 01922 } 01923 01924 01925 #ifdef DOXYGEN // hide AccumulatorChainImpl from documentation 01926 01927 /** Set options for all histograms in the accumulator chain. See histogram accumulators for possible options. The function is ignored if there is no histogram in the accumulator chain. 01928 */ 01929 void setHistogramOptions(HistogramOptions const & options); 01930 01931 /** Reset current_pass_ of the accumulator chain to 'reset_to_pass'. */ 01932 void reset(unsigned int reset_to_pass = 0); 01933 01934 /** Equivalent to merge(o) . */ 01935 void operator+=(AccumulatorChainImpl const & o); 01936 01937 /** Merge the accumulator chain with accumulator chain 'o'. This only works if all selected statistics in the accumulator chain support the '+=' operator. See the documentations of the particular statistics for support information. 01938 */ 01939 void merge(AccumulatorChainImpl const & o); 01940 01941 /** Upate all accumulators in the accumulator chain that work in pass N with data t. Requirement: 0 < N < 6 and N >= current_pass_ . If N < current_pass_ call reset first. 01942 */ 01943 void updatePassN(T const & t, unsigned int N); 01944 01945 /** Upate all accumulators in the accumulator chain that work in pass N with data t and weight. Requirement: 0 < N < 6 and N >= current_pass_ . If N < current_pass_ call reset first. 01946 */ 01947 void updatePassN(T const & t, double weight, unsigned int N); 01948 01949 /** Return the number of passes required to compute all statistics in the accumulator chain. 01950 */ 01951 unsigned int passesRequired() const; 01952 01953 #endif 01954 01955 private: 01956 static ArrayVector<std::string> collectTagNames() 01957 { 01958 ArrayVector<std::string> n; 01959 detail::CollectAccumulatorNames<AccumulatorTags>::exec(n); 01960 std::sort(n.begin(), n.end()); 01961 return n; 01962 } 01963 }; 01964 01965 01966 // Create a dynamic accumulator chain containing the Selected statistics and their dependencies. 01967 // Statistics will only be computed if activate<Tag>() is called at runtime. 01968 /** \brief Create a dynamic accumulator chain containing the selected statistics and their dependencies. 01969 01970 DynamicAccumulatorChain is used to compute global statistics with run-time activation. A set of statistics is selected at run-time and from this set statistics can be activated at run-time by calling activate<stat>() or activate(std::string stat). 01971 01972 The template parameters are as follows: 01973 - T: The input type 01974 - either element type of the data(e.g. double, int, RGBValue, ...) 01975 - or type of CoupledHandle (for access to coordinates and/or weights) 01976 - Selected: statistics to be computed and index specifier for the CoupledHandle, wrapped with Select 01977 01978 Usage: 01979 \code 01980 typedef double DataType; 01981 DynamicAccumulatorChain<DataType, Select<Variance, Mean, Minimum, ...> > accumulator; 01982 \endcode 01983 01984 Usage, using CoupledHandle: 01985 \code 01986 const int dim = 3; //dimension of MultiArray 01987 typedef double DataType; 01988 typedef double WeightType; 01989 typedef vigra::CoupledIteratorType<dim, DataType, WeightType>::HandleType Handle; 01990 DynamicAccumulatorChain<Handle, Select<DataArg<1>, WeightArg<2>, Mean,...> > a; 01991 \endcode 01992 01993 See \ref FeatureAccumulators for more information and examples of use. 01994 */ 01995 template <class T, class Selected> 01996 class DynamicAccumulatorChain 01997 : public AccumulatorChain<T, Selected, true> 01998 { 01999 public: 02000 typedef typename AccumulatorChain<T, Selected, true>::InternalBaseType InternalBaseType; 02001 typedef typename DynamicAccumulatorChain::AccumulatorTags AccumulatorTags; 02002 02003 /** Activate statistic 'tag'. Alias names are not recognized. If the statistic is not in the accumulator chain a PreconditionViolation is thrown. 02004 */ 02005 void activate(std::string tag) 02006 { 02007 vigra_precondition(activateImpl(tag), 02008 std::string("DynamicAccumulatorChain::activate(): Tag '") + tag + "' not found."); 02009 } 02010 02011 /** %activate<TAG>() activates statistic 'TAG'. If the statistic is not in the accumulator chain it is ignored. (?) 02012 */ 02013 template <class TAG> 02014 void activate() 02015 { 02016 LookupTag<TAG, DynamicAccumulatorChain>::type::activateImpl(getAccumulator<AccumulatorEnd>(*this).active_accumulators_); 02017 } 02018 02019 /** Activate all statistics in the accumulator chain. 02020 */ 02021 void activateAll() 02022 { 02023 getAccumulator<AccumulatorEnd>(*this).active_accumulators_.set(); 02024 } 02025 /** Return true if the statistic 'tag' is active, i.e. activate(std::string tag) or activate<TAG>() has been called. If the statistic is not in the accumulator chain a PreconditionViolation is thrown. (Note that alias names are not recognized.) 02026 */ 02027 bool isActive(std::string tag) const 02028 { 02029 detail::TagIsActive_Visitor v; 02030 vigra_precondition(isActiveImpl(tag, v), 02031 std::string("DynamicAccumulatorChain::isActive(): Tag '") + tag + "' not found."); 02032 return v.result; 02033 } 02034 02035 /** %isActive<TAG>() returns true if statistic 'TAG' is active, i.e. activate(std::string tag) or activate<TAG>() has been called. If the statistic is not in the accumulator chain, true is returned. (?) 02036 */ 02037 template <class TAG> 02038 bool isActive() const 02039 { 02040 return LookupTag<TAG, DynamicAccumulatorChain>::type::isActiveImpl(getAccumulator<AccumulatorEnd>(*this).active_accumulators_); 02041 } 02042 02043 /** Return names of all statistics in the accumulator chain that are active. 02044 */ 02045 ArrayVector<std::string> activeNames() const 02046 { 02047 ArrayVector<std::string> res; 02048 for(unsigned k=0; k<DynamicAccumulatorChain::tagNames().size(); ++k) 02049 if(isActive(DynamicAccumulatorChain::tagNames()[k])) 02050 res.push_back(DynamicAccumulatorChain::tagNames()[k]); 02051 return res; 02052 } 02053 02054 /** Return number of passes required to compute the active statistics in the accumulator chain. 02055 */ 02056 unsigned int passesRequired() const 02057 { 02058 return InternalBaseType::passesRequired(getAccumulator<AccumulatorEnd>(*this).active_accumulators_); 02059 } 02060 02061 protected: 02062 02063 bool activateImpl(std::string tag) 02064 { 02065 return detail::ApplyVisitorToTag<AccumulatorTags>::exec(*this, 02066 normalizeString(tag), detail::ActivateTag_Visitor()); 02067 } 02068 02069 bool isActiveImpl(std::string tag, detail::TagIsActive_Visitor & v) const 02070 { 02071 return detail::ApplyVisitorToTag<AccumulatorTags>::exec(*this, normalizeString(tag), v); 02072 } 02073 }; 02074 02075 /** \brief Create an array of accumulator chains containing the selected per-region and global statistics and their dependencies. 02076 02077 AccumulatorChainArray is used to compute per-region statistics (as well as global statistics). The statistics are selected at compile-time. An array of accumulator chains (one per region) for region statistics is created and one accumulator chain for global statistics. The region labels always start at 0. Use the Global modifier to compute global statistics (by default per-region statistics are computed). 02078 02079 The template parameters are as follows: 02080 - T: The input type, type of CoupledHandle (for access to coordinates, labels and weights) 02081 - Selected: statistics to be computed and index specifier for the CoupledHandle, wrapped with Select 02082 02083 Usage: 02084 \code 02085 const int dim = 3; //dimension of MultiArray 02086 typedef double DataType; 02087 typedef double WeightType; 02088 typedef unsigned int LabelType; 02089 typedef vigra::CoupledIteratorType<dim, DataType, WeightType, LabelType>::HandleType Handle; 02090 AccumulatorChainArray<Handle, Select<DataArg<1>, WeightArg<2>, LabelArg<3>, Mean, Variance, ...> > a; 02091 \endcode 02092 02093 See \ref FeatureAccumulators for more information and examples of use. 02094 */ 02095 template <class T, class Selected, bool dynamic=false> 02096 class AccumulatorChainArray 02097 #ifndef DOXYGEN //hide AccumulatorChainImpl vom documentation 02098 : public AccumulatorChainImpl<T, typename detail::ConfigureAccumulatorChainArray<T, Selected, dynamic>::type> 02099 #endif 02100 { 02101 public: 02102 typedef typename detail::ConfigureAccumulatorChainArray<T, Selected, dynamic> Creator; 02103 typedef typename Creator::TagList AccumulatorTags; 02104 typedef typename Creator::GlobalTags GlobalTags; 02105 typedef typename Creator::RegionTags RegionTags; 02106 02107 /** Statistics will not be computed for label l. Note that only one label can be ignored. 02108 */ 02109 void ignoreLabel(MultiArrayIndex l) 02110 { 02111 this->next_.ignoreLabel(l); 02112 } 02113 02114 /** Set the maximum region label (e.g. for merging two accumulator chains). 02115 */ 02116 void setMaxRegionLabel(unsigned label) 02117 { 02118 this->next_.setMaxRegionLabel(label); 02119 } 02120 02121 /** %Maximum region label. (equal to regionCount() - 1) 02122 */ 02123 MultiArrayIndex maxRegionLabel() const 02124 { 02125 return this->next_.maxRegionLabel(); 02126 } 02127 02128 /** Number of Regions. (equal to maxRegionLabel() + 1) 02129 */ 02130 unsigned int regionCount() const 02131 { 02132 return this->next_.regions_.size(); 02133 } 02134 02135 /** Merge region i with region j. 02136 */ 02137 void merge(unsigned i, unsigned j) 02138 { 02139 vigra_precondition(i <= maxRegionLabel() && j <= maxRegionLabel(), 02140 "AccumulatorChainArray::merge(): region labels out of range."); 02141 this->next_.merge(i, j); 02142 } 02143 02144 /** Merge with accumulator chain o. maxRegionLabel() of the two accumulators must be equal. 02145 */ 02146 void merge(AccumulatorChainArray const & o) 02147 { 02148 if(maxRegionLabel() == -1) 02149 setMaxRegionLabel(o.maxRegionLabel()); 02150 vigra_precondition(maxRegionLabel() == o.maxRegionLabel(), 02151 "AccumulatorChainArray::merge(): maxRegionLabel must be equal."); 02152 this->next_.merge(o.next_); 02153 } 02154 02155 /** Merge with accumulator chain o using a mapping between labels of the two accumulators. Label l of accumulator chain o is mapped to labelMapping[l]. Hence, all elements of labelMapping must be <= maxRegionLabel() and size of labelMapping must match o.regionCount(). 02156 */ 02157 template <class ArrayLike> 02158 void merge(AccumulatorChainArray const & o, ArrayLike const & labelMapping) 02159 { 02160 vigra_precondition(labelMapping.size() == o.regionCount(), 02161 "AccumulatorChainArray::merge(): labelMapping.size() must match regionCount() of RHS."); 02162 this->next_.merge(o.next_, labelMapping); 02163 } 02164 02165 /** Return names of all tags in the accumulator chain (selected statistics and their dependencies). 02166 */ 02167 static ArrayVector<std::string> const & tagNames() 02168 { 02169 static const ArrayVector<std::string> n = collectTagNames(); 02170 return n; 02171 } 02172 02173 02174 #ifdef DOXYGEN // hide AccumulatorChainImpl from documentation 02175 02176 /** \copydoc AccumulatorChain::setHistogramOptions(HistogramOptions const &) */ 02177 void setHistogramOptions(HistogramOptions const & options); 02178 02179 /** Set regional and global options for all histograms in the accumulator chain. 02180 */ 02181 void setHistogramOptions(HistogramOptions const & regionoptions, HistogramOptions const & globaloptions); 02182 02183 /** \copydoc AccumulatorChain::reset() */ 02184 void reset(unsigned int reset_to_pass = 0); 02185 02186 /** \copydoc AccumulatorChain::operator+=() */ 02187 void operator+=(AccumulatorChainImpl const & o); 02188 02189 /** \copydoc AccumulatorChain::updatePassN(T const &,unsigned int) */ 02190 void updatePassN(T const & t, unsigned int N); 02191 02192 /** \copydoc AccumulatorChain::updatePassN(T const &,double,unsigned int) */ 02193 void updatePassN(T const & t, double weight, unsigned int N); 02194 02195 #endif 02196 02197 02198 02199 private: 02200 static ArrayVector<std::string> collectTagNames() 02201 { 02202 ArrayVector<std::string> n; 02203 detail::CollectAccumulatorNames<AccumulatorTags>::exec(n); 02204 std::sort(n.begin(), n.end()); 02205 return n; 02206 } 02207 }; 02208 02209 /** \brief Create an array of dynamic accumulator chains containing the selected per-region and global statistics and their dependencies. 02210 02211 02212 DynamicAccumulatorChainArray is used to compute per-region statistics (as well as global statistics) with run-time activation. A set of statistics is selected at run-time and from this set statistics can be activated at run-time by calling activate<stat>() or activate(std::string stat). 02213 02214 The template parameters are as follows: 02215 - T: The input type, type of CoupledHandle (for access to coordinates, labels and weights) 02216 - Selected: statistics to be computed and index specifier for the CoupledHandle, wrapped with Select 02217 02218 Usage: 02219 \code 02220 const int dim = 3; //dimension of MultiArray 02221 typedef double DataType; 02222 typedef double WeightType; 02223 typedef unsigned int LabelType; 02224 typedef vigra::CoupledIteratorType<dim, DataType, WeightType, LabelType>::HandleType Handle; 02225 DynamicAccumulatorChainArray<Handle, Select<DataArg<1>, WeightArg<2>, LabelArg<3>, Mean, Variance, ...> > a; 02226 \endcode 02227 02228 See \ref FeatureAccumulators for more information and examples of use. 02229 */ 02230 template <class T, class Selected> 02231 class DynamicAccumulatorChainArray 02232 : public AccumulatorChainArray<T, Selected, true> 02233 { 02234 public: 02235 typedef typename DynamicAccumulatorChainArray::AccumulatorTags AccumulatorTags; 02236 02237 /** \copydoc DynamicAccumulatorChain::activate(std::string tag) */ 02238 void activate(std::string tag) 02239 { 02240 vigra_precondition(activateImpl(tag), 02241 std::string("DynamicAccumulatorChainArray::activate(): Tag '") + tag + "' not found."); 02242 } 02243 02244 /** \copydoc DynamicAccumulatorChain::activate() */ 02245 template <class TAG> 02246 void activate() 02247 { 02248 this->next_.template activate<TAG>(); 02249 } 02250 02251 /** \copydoc DynamicAccumulatorChain::activateAll() */ 02252 void activateAll() 02253 { 02254 this->next_.activateAll(); 02255 } 02256 02257 /** Return true if the statistic 'tag' is active, i.e. activate(std::string tag) or activate<TAG>() has been called. If the statistic is not in the accumulator chain a PreconditionViolation is thrown. (Note that alias names are not recognized.) 02258 */ 02259 bool isActive(std::string tag) const 02260 { 02261 detail::TagIsActive_Visitor v; 02262 vigra_precondition(isActiveImpl(tag, v), 02263 std::string("DynamicAccumulatorChainArray::isActive(): Tag '") + tag + "' not found."); 02264 return v.result; 02265 } 02266 02267 /** %isActive<TAG>() returns true if statistic 'TAG' is active, i.e. activate(std::string tag) or activate<TAG>() has been called. If the statistic is not in the accumulator chain, true is returned. (?) 02268 */ 02269 template <class TAG> 02270 bool isActive() const 02271 { 02272 return this->next_.template isActive<TAG>(); 02273 } 02274 02275 /** \copydoc DynamicAccumulatorChain::activeNames() */ 02276 ArrayVector<std::string> activeNames() const 02277 { 02278 ArrayVector<std::string> res; 02279 for(unsigned k=0; k<DynamicAccumulatorChainArray::tagNames().size(); ++k) 02280 if(isActive(DynamicAccumulatorChainArray::tagNames()[k])) 02281 res.push_back(DynamicAccumulatorChainArray::tagNames()[k]); 02282 return res; 02283 } 02284 02285 /** \copydoc DynamicAccumulatorChain::passesRequired() */ 02286 unsigned int passesRequired() const 02287 { 02288 return this->next_.passesRequiredDynamic(); 02289 } 02290 02291 protected: 02292 02293 bool activateImpl(std::string tag) 02294 { 02295 return detail::ApplyVisitorToTag<AccumulatorTags>::exec(this->next_, 02296 normalizeString(tag), detail::ActivateTag_Visitor()); 02297 } 02298 02299 bool isActiveImpl(std::string tag, detail::TagIsActive_Visitor & v) const 02300 { 02301 return detail::ApplyVisitorToTag<AccumulatorTags>::exec(this->next_, normalizeString(tag), v); 02302 } 02303 }; 02304 02305 /****************************************************************************/ 02306 /* */ 02307 /* generic access functions */ 02308 /* */ 02309 /****************************************************************************/ 02310 02311 template <class TAG> 02312 struct Error__Attempt_to_access_inactive_statistic; 02313 02314 namespace detail { 02315 02316 // accumulator lookup rules: find the accumulator that implements TAG 02317 02318 // When A does not implement TAG, continue search in A::InternalBaseType. 02319 template <class TAG, class A, class FromTag=typename A::Tag> 02320 struct LookupTagImpl 02321 #ifndef DOXYGEN 02322 : public LookupTagImpl<TAG, typename A::InternalBaseType> 02323 #endif 02324 {}; 02325 02326 // 'const A' is treated like A, except that the reference member is now const. 02327 template <class TAG, class A, class FromTag> 02328 struct LookupTagImpl<TAG, A const, FromTag> 02329 : public LookupTagImpl<TAG, A> 02330 { 02331 typedef typename LookupTagImpl<TAG, A>::type const & reference; 02332 typedef typename LookupTagImpl<TAG, A>::type const * pointer; 02333 }; 02334 02335 // When A implements TAG, report its type and associated information. 02336 template <class TAG, class A> 02337 struct LookupTagImpl<TAG, A, TAG> 02338 { 02339 typedef TAG Tag; 02340 typedef A type; 02341 typedef A & reference; 02342 typedef A * pointer; 02343 typedef typename A::value_type value_type; 02344 typedef typename A::result_type result_type; 02345 }; 02346 02347 // Again, 'const A' is treated like A, except that the reference member is now const. 02348 template <class TAG, class A> 02349 struct LookupTagImpl<TAG, A const, TAG> 02350 : public LookupTagImpl<TAG, A, TAG> 02351 { 02352 typedef typename LookupTagImpl<TAG, A, TAG>::type const & reference; 02353 typedef typename LookupTagImpl<TAG, A, TAG>::type const * pointer; 02354 }; 02355 02356 // Recursion termination: when we end up in AccumulatorEnd without finding a 02357 // suitable A, we stop and report an error 02358 template <class TAG, class A> 02359 struct LookupTagImpl<TAG, A, AccumulatorEnd> 02360 { 02361 typedef TAG Tag; 02362 typedef A type; 02363 typedef A & reference; 02364 typedef A * pointer; 02365 typedef Error__Attempt_to_access_inactive_statistic<TAG> value_type; 02366 typedef Error__Attempt_to_access_inactive_statistic<TAG> result_type; 02367 }; 02368 02369 // ... except when we are actually looking for AccumulatorEnd 02370 template <class A> 02371 struct LookupTagImpl<AccumulatorEnd, A, AccumulatorEnd> 02372 { 02373 typedef AccumulatorEnd Tag; 02374 typedef A type; 02375 typedef A & reference; 02376 typedef A * pointer; 02377 typedef void value_type; 02378 typedef void result_type; 02379 }; 02380 02381 // ... or we are looking for a global statistic, in which case 02382 // we continue the serach via A::GlobalAccumulatorType, but remember that 02383 // we are actually looking for a global tag. 02384 template <class TAG, class A> 02385 struct LookupTagImpl<Global<TAG>, A, AccumulatorEnd> 02386 : public LookupTagImpl<TAG, typename A::GlobalAccumulatorType> 02387 { 02388 typedef Global<TAG> Tag; 02389 }; 02390 02391 // When we encounter the LabelDispatch accumulator, we continue the 02392 // search via LabelDispatch::RegionAccumulatorChain by default 02393 template <class TAG, class A> 02394 struct LookupTagImpl<TAG, A, LabelDispatchTag> 02395 : public LookupTagImpl<TAG, typename A::RegionAccumulatorChain> 02396 {}; 02397 02398 // ... except when we are looking for a global statistic, in which case 02399 // we continue via LabelDispatch::GlobalAccumulatorChain, but remember that 02400 // we are actually looking for a global tag. 02401 template <class TAG, class A> 02402 struct LookupTagImpl<Global<TAG>, A, LabelDispatchTag> 02403 : public LookupTagImpl<TAG, typename A::GlobalAccumulatorChain> 02404 { 02405 typedef Global<TAG> Tag; 02406 }; 02407 02408 // ... or we are looking for the LabelDispatch accumulator itself 02409 template <class A> 02410 struct LookupTagImpl<LabelDispatchTag, A, LabelDispatchTag> 02411 { 02412 typedef LabelDispatchTag Tag; 02413 typedef A type; 02414 typedef A & reference; 02415 typedef A * pointer; 02416 typedef void value_type; 02417 typedef void result_type; 02418 }; 02419 02420 } // namespace detail 02421 02422 // Lookup the accumulator in the chain A that implements the given TAG. 02423 template <class Tag, class A> 02424 struct LookupTag 02425 : public detail::LookupTagImpl<typename StandardizeTag<Tag>::type, A> 02426 {}; 02427 02428 // Lookup the dependency TAG of the accumulator A. 02429 // This template ensures that dependencies are used with matching modifiers. 02430 // Specifically, if you search for Count as a dependency of Weighted<Mean>, the search 02431 // actually returns Weighted<Count>, wheras Count will be returned for plain Mean. 02432 template <class Tag, class A, class TargetTag> 02433 struct LookupDependency 02434 : public detail::LookupTagImpl< 02435 typename TransferModifiers<TargetTag, typename StandardizeTag<Tag>::type>::type, A> 02436 {}; 02437 02438 02439 namespace detail { 02440 02441 // CastImpl applies the same rules as LookupTagImpl, but returns a reference to an 02442 // accumulator instance rather than an accumulator type 02443 template <class Tag, class FromTag, class reference> 02444 struct CastImpl 02445 { 02446 template <class A> 02447 static reference exec(A & a) 02448 { 02449 return CastImpl<Tag, typename A::InternalBaseType::Tag, reference>::exec(a.next_); 02450 } 02451 02452 template <class A> 02453 static reference exec(A & a, MultiArrayIndex label) 02454 { 02455 return CastImpl<Tag, typename A::InternalBaseType::Tag, reference>::exec(a.next_, label); 02456 } 02457 }; 02458 02459 template <class Tag, class reference> 02460 struct CastImpl<Tag, Tag, reference> 02461 { 02462 template <class A> 02463 static reference exec(A & a) 02464 { 02465 return const_cast<reference>(a); 02466 } 02467 02468 template <class A> 02469 static reference exec(A & a, MultiArrayIndex) 02470 { 02471 vigra_precondition(false, 02472 "getAccumulator(): region accumulators can only be queried for AccumulatorChainArray."); 02473 return a; 02474 } 02475 }; 02476 02477 template <class Tag, class reference> 02478 struct CastImpl<Tag, AccumulatorEnd, reference> 02479 { 02480 template <class A> 02481 static reference exec(A & a) 02482 { 02483 return a; 02484 } 02485 02486 template <class A> 02487 static reference exec(A & a, MultiArrayIndex) 02488 { 02489 return a; 02490 } 02491 }; 02492 02493 template <class Tag, class reference> 02494 struct CastImpl<Global<Tag>, AccumulatorEnd, reference> 02495 { 02496 template <class A> 02497 static reference exec(A & a) 02498 { 02499 return CastImpl<Tag, typename A::GlobalAccumulatorType::Tag, reference>::exec(*a.globalAccumulator_.pointer_); 02500 } 02501 }; 02502 02503 template <class reference> 02504 struct CastImpl<AccumulatorEnd, AccumulatorEnd, reference> 02505 { 02506 template <class A> 02507 static reference exec(A & a) 02508 { 02509 return a; 02510 } 02511 02512 template <class A> 02513 static reference exec(A & a, MultiArrayIndex) 02514 { 02515 return a; 02516 } 02517 }; 02518 02519 template <class Tag, class reference> 02520 struct CastImpl<Tag, LabelDispatchTag, reference> 02521 { 02522 template <class A> 02523 static reference exec(A & a) 02524 { 02525 vigra_precondition(false, 02526 "getAccumulator(): a region label is required when a region accumulator is queried."); 02527 return CastImpl<Tag, typename A::RegionAccumulatorChain::Tag, reference>::exec(a.regions_[0]); 02528 } 02529 02530 template <class A> 02531 static reference exec(A & a, MultiArrayIndex label) 02532 { 02533 return CastImpl<Tag, typename A::RegionAccumulatorChain::Tag, reference>::exec(a.regions_[label]); 02534 } 02535 }; 02536 02537 template <class Tag, class reference> 02538 struct CastImpl<Global<Tag>, LabelDispatchTag, reference> 02539 { 02540 template <class A> 02541 static reference exec(A & a) 02542 { 02543 return CastImpl<Tag, typename A::GlobalAccumulatorChain::Tag, reference>::exec(a.next_); 02544 } 02545 }; 02546 02547 template <class reference> 02548 struct CastImpl<LabelDispatchTag, LabelDispatchTag, reference> 02549 { 02550 template <class A> 02551 static reference exec(A & a) 02552 { 02553 return a; 02554 } 02555 }; 02556 02557 } // namespace detail 02558 02559 // Get a reference to the accumulator TAG in the accumulator chain A 02560 /** Get a reference to the accumulator 'TAG' in the accumulator chain 'a'. This can be useful for example to update a certain accumulator with data, set individual options or get information about a certain accumulator.\n 02561 Example of use (set options): 02562 \code 02563 vigra::MultiArray<2, double> data(...); 02564 typedef UserRangeHistogram<40> SomeHistogram; //binCount set at compile time 02565 typedef UserRangeHistogram<0> SomeHistogram2; // binCount must be set at run-time 02566 AccumulatorChain<DataType, Select<SomeHistogram, SomeHistogram2> > a; 02567 02568 getAccumulator<SomeHistogram>(a).setMinMax(0.1, 0.9); 02569 getAccumulator<SomeHistogram2>(a).setMinMax(0.0, 1.0); 02570 02571 extractFeatures(data.begin(), data.end(), a); 02572 \endcode 02573 02574 Example of use (get information): 02575 \code 02576 vigra::MultiArray<2, double> data(...)); 02577 AccumulatorChain<double, Select<Mean, Skewness> > a; 02578 02579 std::cout << "passes required for all statistics: " << a.passesRequired() << std::endl; //skewness needs two passes 02580 std::cout << "passes required by Mean: " << getAccumulator<Mean>(a).passesRequired() << std::endl; 02581 \endcode 02582 See \ref FeatureAccumulators for more information about feature computation via accumulators. 02583 */ 02584 template <class TAG, class A> 02585 inline typename LookupTag<TAG, A>::reference 02586 getAccumulator(A & a) 02587 { 02588 typedef typename LookupTag<TAG, A>::Tag StandardizedTag; 02589 typedef typename LookupTag<TAG, A>::reference reference; 02590 return detail::CastImpl<StandardizedTag, typename A::Tag, reference>::exec(a); 02591 } 02592 02593 // Get a reference to the accumulator TAG for region 'label' in the accumulator chain A 02594 /** Get a reference to the accumulator 'TAG' for region 'label' in the accumulator chain 'a'. 02595 */ 02596 template <class TAG, class A> 02597 inline typename LookupTag<TAG, A>::reference 02598 getAccumulator(A & a, MultiArrayIndex label) 02599 { 02600 typedef typename LookupTag<TAG, A>::Tag StandardizedTag; 02601 typedef typename LookupTag<TAG, A>::reference reference; 02602 return detail::CastImpl<StandardizedTag, typename A::Tag, reference>::exec(a, label); 02603 } 02604 02605 // get the result of the accumulator specified by TAG 02606 /** Get the result of the accumulator 'TAG' in the accumulator chain 'a'.\n 02607 Example of use: 02608 \code 02609 vigra::MultiArray<2, double> data(...); 02610 AccumulatorChain<DataType, Select<Variance, Mean, StdDev> > a; 02611 extractFeatures(data.begin(), data.end(), a); 02612 double mean = get<Mean>(a); 02613 \endcode 02614 See \ref FeatureAccumulators for more information about feature computation via accumulators. 02615 */ 02616 template <class TAG, class A> 02617 inline typename LookupTag<TAG, A>::result_type 02618 get(A const & a) 02619 { 02620 return getAccumulator<TAG>(a).get(); 02621 } 02622 02623 // get the result of the accumulator TAG for region 'label' 02624 /** Get the result of the accumulator 'TAG' for region 'label' in the accumulator chain 'a'.\n 02625 Example of use: 02626 \code 02627 vigra::MultiArray<2, double> data(...); 02628 vigra::MultiArray<2, int> labels(...); 02629 typedef vigra::CoupledIteratorType<2, double, int>::type Iterator; 02630 typedef Iterator::value_type Handle; 02631 02632 AccumulatorChainArray<Handle, 02633 Select<DataArg<1>, LabelArg<2>, Mean, Variance> > a; 02634 02635 Iterator start = createCoupledIterator(data, labels); 02636 Iterator end = start.getEndIterator(); 02637 extractFeatures(start,end,a); 02638 02639 double mean_of_region_1 = get<Mean>(a,1); 02640 double mean_of_background = get<Mean>(a,0); 02641 \endcode 02642 See \ref FeatureAccumulators for more information about feature computation via accumulators. 02643 */ 02644 template <class TAG, class A> 02645 inline typename LookupTag<TAG, A>::result_type 02646 get(A const & a, MultiArrayIndex label) 02647 { 02648 return getAccumulator<TAG>(a, label).get(); 02649 } 02650 02651 // Get the result of the accumulator specified by TAG without checking if the accumulator is active. 02652 // This must be used within an accumulator implementation to access dependencies because 02653 // it applies the approprate modifiers to the given TAG. It must not be used in other situations. 02654 // FIXME: is there a shorter name? 02655 template <class TAG, class A> 02656 inline typename LookupDependency<TAG, A>::result_type 02657 getDependency(A const & a) 02658 { 02659 typedef typename LookupDependency<TAG, A>::Tag StandardizedTag; 02660 typedef typename LookupDependency<TAG, A>::reference reference; 02661 return detail::CastImpl<StandardizedTag, typename A::Tag, reference>::exec(a)(); 02662 } 02663 02664 // activate the dynamic accumulator specified by Tag 02665 /** Activate the dynamic accumulator 'Tag' in the dynamic accumulator chain 'a'. Same as a.activate<Tag>() (see DynamicAccumulatorChain::activate<Tag>() or DynamicAccumulatorChainArray::activate<Tag>()). For run-time activation use DynamicAccumulatorChain::activate(std::string tag) or DynamicAccumulatorChainArray::activate(std::string tag) instead.\n 02666 See \ref FeatureAccumulators for more information about feature computation via accumulators. 02667 */ 02668 template <class Tag, class A> 02669 inline void 02670 activate(A & a) 02671 { 02672 a.template activate<Tag>(); 02673 } 02674 02675 // check if the dynamic accumulator specified by Tag is active 02676 /** Check if the dynamic accumulator 'Tag' in the accumulator chain 'a' is active. Same as a.isActive<Tag>() (see DynamicAccumulatorChain::isActive<Tag>() or DynamicAccumulatorChainArray::isActive<Tag>()). At run-time, use DynamicAccumulatorChain::isActive(std::string tag) const or DynamicAccumulatorChainArray::isActive(std::string tag) const instead.\n 02677 See \ref FeatureAccumulators for more information about feature computation via accumulators. 02678 */ 02679 template <class Tag, class A> 02680 inline bool 02681 isActive(A const & a) 02682 { 02683 return a.template isActive<Tag>(); 02684 } 02685 02686 /****************************************************************************/ 02687 /* */ 02688 /* generic loops */ 02689 /* */ 02690 /****************************************************************************/ 02691 02692 /** Generic loop to collect the statistics of the accumulator chain 'a' in as many passes over the data as necessary.\n 02693 02694 Example of use: 02695 \code 02696 vigra::MultiArray<3, double> data(...); 02697 vigra::MultiArray<3, int> labels(...); 02698 typedef vigra::CoupledIteratorType<3, double, int>::type Iterator; 02699 typedef Iterator::value_type Handle; 02700 02701 AccumulatorChainArray<Handle, 02702 Select<DataArg<1>, LabelArg<2>, Mean, Variance> > a; 02703 02704 Iterator start = createCoupledIterator(data, labels); 02705 Iterator end = start.getEndIterator(); 02706 02707 extractFeatures(start,end,a); 02708 \endcode 02709 See \ref FeatureAccumulators for more information about feature computation via accumulators. 02710 */ 02711 template <class ITERATOR, class ACCUMULATOR> 02712 void extractFeatures(ITERATOR start, ITERATOR end, ACCUMULATOR & a) 02713 { 02714 for(unsigned int k=1; k <= a.passesRequired(); ++k) 02715 for(ITERATOR i=start; i < end; ++i) 02716 a.updatePassN(*i, k); 02717 } 02718 02719 /****************************************************************************/ 02720 /* */ 02721 /* AccumulatorResultTraits */ 02722 /* */ 02723 /****************************************************************************/ 02724 02725 template <class T> 02726 struct AccumulatorResultTraits 02727 { 02728 typedef T type; 02729 typedef T element_type; 02730 typedef double element_promote_type; 02731 typedef T MinmaxType; 02732 typedef element_promote_type SumType; 02733 typedef element_promote_type FlatCovarianceType; 02734 typedef element_promote_type CovarianceType; 02735 }; 02736 02737 template <class T, int N> 02738 struct AccumulatorResultTraits<TinyVector<T, N> > 02739 { 02740 typedef TinyVector<T, N> type; 02741 typedef T element_type; 02742 typedef double element_promote_type; 02743 typedef TinyVector<T, N> MinmaxType; 02744 typedef TinyVector<element_promote_type, N> SumType; 02745 typedef TinyVector<element_promote_type, N*(N+1)/2> FlatCovarianceType; 02746 typedef Matrix<element_promote_type> CovarianceType; 02747 }; 02748 02749 // (?) beign change 02750 template <class T, unsigned int RED_IDX, unsigned int GREEN_IDX, unsigned int BLUE_IDX> 02751 struct AccumulatorResultTraits<RGBValue<T, RED_IDX, GREEN_IDX, BLUE_IDX> > 02752 { 02753 typedef RGBValue<T> type; 02754 typedef T element_type; 02755 typedef double element_promote_type; 02756 typedef RGBValue<T> MinmaxType; 02757 typedef RGBValue<element_promote_type> SumType; 02758 typedef TinyVector<element_promote_type, 3*(3+1)/2> FlatCovarianceType; 02759 typedef Matrix<element_promote_type> CovarianceType; 02760 }; 02761 // end change 02762 02763 02764 template <unsigned int N, class T, class Stride> 02765 struct AccumulatorResultTraits<MultiArrayView<N, T, Stride> > 02766 { 02767 typedef MultiArrayView<N, T, Stride> type; 02768 typedef T element_type; 02769 typedef double element_promote_type; 02770 typedef MultiArray<N, T> MinmaxType; 02771 typedef MultiArray<N, element_promote_type> SumType; 02772 typedef MultiArray<1, element_promote_type> FlatCovarianceType; 02773 typedef Matrix<element_promote_type> CovarianceType; 02774 }; 02775 02776 template <unsigned int N, class T, class Alloc> 02777 struct AccumulatorResultTraits<MultiArray<N, T, Alloc> > 02778 { 02779 typedef MultiArrayView<N, T, Alloc> type; 02780 typedef T element_type; 02781 typedef double element_promote_type; 02782 typedef MultiArray<N, T> MinmaxType; 02783 typedef MultiArray<N, element_promote_type> SumType; 02784 typedef MultiArray<1, element_promote_type> FlatCovarianceType; 02785 typedef Matrix<element_promote_type> CovarianceType; 02786 }; 02787 02788 /****************************************************************************/ 02789 /* */ 02790 /* modifier implementations */ 02791 /* */ 02792 /****************************************************************************/ 02793 02794 /** \brief Modifier. Compute statistic globally rather than per region. 02795 02796 This modifier only works when labels are given (with (Dynamic)AccumulatorChainArray), in which case statistics are computed per-region by default. 02797 */ 02798 template <class TAG> 02799 class Global 02800 { 02801 public: 02802 typedef typename StandardizeTag<TAG>::type TargetTag; 02803 typedef typename TargetTag::Dependencies Dependencies; 02804 02805 static std::string const & name() 02806 { 02807 static const std::string n = std::string("Global<") + TargetTag::name() + " >"; 02808 return n; 02809 } 02810 }; 02811 02812 /** \brief Specifies index of data in CoupledHandle. 02813 02814 If AccumulatorChain is used with CoupledIterator, DataArg<INDEX> tells the accumulator chain which index of the Handle contains the data. (Coordinates are always index 0) 02815 */ 02816 template <int INDEX> 02817 class DataArg 02818 { 02819 public: 02820 typedef Select<> Dependencies; 02821 02822 static std::string const & name() 02823 { 02824 static const std::string n = std::string("DataArg<") + asString(INDEX) + "> (internal)"; 02825 return n; 02826 } 02827 02828 template <class T, class BASE> 02829 struct Impl 02830 : public BASE 02831 { 02832 typedef DataArgTag Tag; 02833 typedef void value_type; 02834 typedef void result_type; 02835 02836 static const int value = INDEX; 02837 static const unsigned int workInPass = 0; 02838 }; 02839 }; 02840 02841 // Tags are automatically wrapped with DataFromHandle if CoupledHandle used 02842 template <class TAG> 02843 class DataFromHandle 02844 { 02845 public: 02846 typedef typename StandardizeTag<TAG>::type TargetTag; 02847 typedef typename TargetTag::Dependencies Dependencies; 02848 02849 static std::string const & name() 02850 { 02851 static const std::string n = std::string("DataFromHandle<") + TargetTag::name() + " > (internal)"; 02852 return n; 02853 } 02854 02855 template <class IndexDefinition, class TagFound=typename IndexDefinition::Tag> 02856 struct DataIndexSelector 02857 { 02858 static const int value = 1; // default: CoupledHandle holds data at index 1 02859 02860 template <class U, class NEXT> 02861 static typename CoupledHandleCast<value, CoupledHandle<U, NEXT> >::type::const_reference 02862 exec(CoupledHandle<U, NEXT> const & t) 02863 { 02864 return get<value>(t); 02865 } 02866 }; 02867 02868 template <class IndexDefinition> 02869 struct DataIndexSelector<IndexDefinition, DataArgTag> 02870 { 02871 static const int value = IndexDefinition::value; 02872 02873 template <class U, class NEXT> 02874 static typename CoupledHandleCast<value, CoupledHandle<U, NEXT> >::type::const_reference 02875 exec(CoupledHandle<U, NEXT> const & t) 02876 { 02877 return get<value>(t); 02878 } 02879 }; 02880 02881 template <class T, class BASE> 02882 struct SelectInputType 02883 { 02884 typedef typename LookupTag<DataArgTag, BASE>::type FindDataIndex; 02885 typedef DataIndexSelector<FindDataIndex> DataIndex; 02886 typedef typename CoupledHandleCast<DataIndex::value, T>::type::value_type type; 02887 }; 02888 02889 template <class T, class BASE> 02890 struct Impl 02891 : public TargetTag::template Impl<typename SelectInputType<T, BASE>::type, BASE> 02892 { 02893 typedef SelectInputType<T, BASE> InputTypeSelector; 02894 typedef typename InputTypeSelector::DataIndex DataIndex; 02895 typedef typename InputTypeSelector::type input_type; 02896 typedef input_type const & argument_type; 02897 typedef argument_type first_argument_type; 02898 02899 typedef typename TargetTag::template Impl<input_type, BASE> ImplType; 02900 02901 using ImplType::reshape; 02902 02903 template <class U, class NEXT> 02904 void reshape(CoupledHandle<U, NEXT> const & t) 02905 { 02906 ImplType::reshape(detail::shapeOf(DataIndex::exec(t))); 02907 } 02908 02909 template <class U, class NEXT> 02910 void update(CoupledHandle<U, NEXT> const & t) 02911 { 02912 ImplType::update(DataIndex::exec(t)); 02913 } 02914 02915 template <class U, class NEXT> 02916 void update(CoupledHandle<U, NEXT> const & t, double weight) 02917 { 02918 ImplType::update(DataIndex::exec(t), weight); 02919 } 02920 }; 02921 }; 02922 02923 /** \brief Modifier. Compute statistic from pixel coordinates rather than from pixel values. 02924 02925 AccumulatorChain must be used with CoupledIterator in order to have access to pixel coordinates. 02926 */ 02927 template <class TAG> 02928 class Coord 02929 { 02930 public: 02931 typedef typename StandardizeTag<TAG>::type TargetTag; 02932 typedef typename TargetTag::Dependencies Dependencies; 02933 02934 static std::string const & name() 02935 { 02936 static const std::string n = std::string("Coord<") + TargetTag::name() + " >"; 02937 return n; 02938 } 02939 02940 template <class IndexDefinition, class TagFound=typename IndexDefinition::Tag> 02941 struct CoordIndexSelector 02942 { 02943 static const int value = 0; // default: CoupledHandle holds coordinates at index 0 02944 02945 template <class U, class NEXT> 02946 static typename CoupledHandleCast<value, CoupledHandle<U, NEXT> >::type::const_reference 02947 exec(CoupledHandle<U, NEXT> const & t) 02948 { 02949 return get<value>(t); 02950 } 02951 }; 02952 02953 template <class IndexDefinition> 02954 struct CoordIndexSelector<IndexDefinition, CoordArgTag> 02955 { 02956 static const int value = IndexDefinition::value; 02957 02958 template <class U, class NEXT> 02959 static typename CoupledHandleCast<value, CoupledHandle<U, NEXT> >::type::const_reference 02960 exec(CoupledHandle<U, NEXT> const & t) 02961 { 02962 return get<value>(t); 02963 } 02964 }; 02965 02966 template <class T, class BASE> 02967 struct SelectInputType 02968 { 02969 typedef typename LookupTag<CoordArgTag, BASE>::type FindDataIndex; 02970 typedef CoordIndexSelector<FindDataIndex> CoordIndex; 02971 typedef typename CoupledHandleCast<CoordIndex::value, T>::type::value_type type; 02972 }; 02973 02974 template <class T, class BASE> 02975 struct Impl 02976 : public TargetTag::template Impl<typename SelectInputType<T, BASE>::type, BASE> 02977 { 02978 typedef SelectInputType<T, BASE> InputTypeSelector; 02979 typedef typename InputTypeSelector::CoordIndex CoordIndex; 02980 typedef typename InputTypeSelector::type input_type; 02981 typedef input_type const & argument_type; 02982 typedef argument_type first_argument_type; 02983 02984 typedef typename TargetTag::template Impl<input_type, BASE> ImplType; 02985 02986 using ImplType::reshape; 02987 02988 template <class U, class NEXT> 02989 void reshape(CoupledHandle<U, NEXT> const & t) 02990 { 02991 ImplType::reshape(detail::shapeOf(CoordIndex::exec(t))); 02992 } 02993 02994 template <class U, class NEXT> 02995 void update(CoupledHandle<U, NEXT> const & t) 02996 { 02997 ImplType::update(CoordIndex::exec(t)); 02998 } 02999 03000 template <class U, class NEXT> 03001 void update(CoupledHandle<U, NEXT> const & t, double weight) 03002 { 03003 ImplType::update(CoordIndex::exec(t), weight); 03004 } 03005 }; 03006 }; 03007 03008 /** \brief Specifies index of data in CoupledHandle. 03009 03010 If AccumulatorChain is used with CoupledIterator, WeightArg<INDEX> tells the accumulator chain which index of the Handle contains the weights. (Note that coordinates are always index 0.) 03011 */ 03012 template <int INDEX> 03013 class WeightArg 03014 { 03015 public: 03016 typedef Select<> Dependencies; 03017 03018 static std::string const & name() 03019 { 03020 static const std::string n = std::string("WeightArg<") + asString(INDEX) + "> (internal)"; 03021 return n; 03022 } 03023 03024 template <class T, class BASE> 03025 struct Impl 03026 : public BASE 03027 { 03028 typedef WeightArgTag Tag; 03029 typedef void value_type; 03030 typedef void result_type; 03031 03032 static const int value = INDEX; 03033 static const unsigned int workInPass = 0; 03034 }; 03035 }; 03036 03037 /** \brief Compute weighted version of the statistic. 03038 */ 03039 template <class TAG> 03040 class Weighted 03041 { 03042 public: 03043 typedef typename StandardizeTag<TAG>::type TargetTag; 03044 typedef typename TargetTag::Dependencies Dependencies; 03045 03046 static std::string const & name() 03047 { 03048 static const std::string n = std::string("Weighted<") + TargetTag::name() + " >"; 03049 return n; 03050 } 03051 03052 template <class IndexDefinition, class TagFound=typename IndexDefinition::Tag> 03053 struct WeightIndexSelector 03054 { 03055 template <class U, class NEXT> 03056 static double exec(CoupledHandle<U, NEXT> const & t) 03057 { 03058 return (double)*t; // default: CoupledHandle holds weights at the last (outermost) index 03059 } 03060 }; 03061 03062 template <class IndexDefinition> 03063 struct WeightIndexSelector<IndexDefinition, WeightArgTag> 03064 { 03065 template <class U, class NEXT> 03066 static double exec(CoupledHandle<U, NEXT> const & t) 03067 { 03068 return (double)get<IndexDefinition::value>(t); 03069 } 03070 }; 03071 03072 template <class T, class BASE> 03073 struct Impl 03074 : public TargetTag::template Impl<T, BASE> 03075 { 03076 typedef typename TargetTag::template Impl<T, BASE> ImplType; 03077 03078 typedef typename LookupTag<WeightArgTag, BASE>::type FindWeightIndex; 03079 03080 template <class U, class NEXT> 03081 void update(CoupledHandle<U, NEXT> const & t) 03082 { 03083 ImplType::update(t, WeightIndexSelector<FindWeightIndex>::exec(t)); 03084 } 03085 }; 03086 }; 03087 03088 // Centralize by subtracting the mean and cache the result 03089 class Centralize 03090 { 03091 public: 03092 typedef Select<Mean> Dependencies; 03093 03094 static std::string const & name() 03095 { 03096 static const std::string n("Centralize (internal)"); 03097 return n; 03098 } 03099 03100 template <class U, class BASE> 03101 struct Impl 03102 : public BASE 03103 { 03104 static const unsigned int workInPass = 2; 03105 03106 typedef typename AccumulatorResultTraits<U>::element_promote_type element_type; 03107 typedef typename AccumulatorResultTraits<U>::SumType value_type; 03108 typedef value_type const & result_type; 03109 03110 mutable value_type value_; 03111 03112 Impl() 03113 : value_() // call default constructor explicitly to ensure zero initialization 03114 {} 03115 03116 void reset() 03117 { 03118 value_ = element_type(); 03119 } 03120 03121 template <class Shape> 03122 void reshape(Shape const & s) 03123 { 03124 detail::reshapeImpl(value_, s); 03125 } 03126 03127 void update(U const & t) const 03128 { 03129 using namespace vigra::multi_math; 03130 value_ = t - getDependency<Mean>(*this); 03131 } 03132 03133 void update(U const & t, double) const 03134 { 03135 update(t); 03136 } 03137 03138 result_type operator()(U const & t) const 03139 { 03140 update(t); 03141 return value_; 03142 } 03143 03144 result_type operator()() const 03145 { 03146 return value_; 03147 } 03148 }; 03149 }; 03150 03151 /** \brief Modifier. Substract mean before computing statistic. 03152 03153 Works in pass 2, %operator+=() not supported (merging not supported). 03154 */ 03155 template <class TAG> 03156 class Central 03157 { 03158 public: 03159 typedef typename StandardizeTag<TAG>::type TargetTag; 03160 typedef Select<Centralize, typename TargetTag::Dependencies> Dependencies; 03161 03162 static std::string const & name() 03163 { 03164 static const std::string n = std::string("Central<") + TargetTag::name() + " >"; 03165 return n; 03166 } 03167 03168 template <class U, class BASE> 03169 struct Impl 03170 : public TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE> 03171 { 03172 typedef typename TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE> ImplType; 03173 03174 static const unsigned int workInPass = 2; 03175 03176 void operator+=(Impl const & o) 03177 { 03178 vigra_precondition(false, 03179 "Central<...>::operator+=(): not supported."); 03180 } 03181 03182 template <class T> 03183 void update(T const & t) 03184 { 03185 ImplType::update(getDependency<Centralize>(*this)); 03186 } 03187 03188 template <class T> 03189 void update(T const & t, double weight) 03190 { 03191 ImplType::update(getDependency<Centralize>(*this), weight); 03192 } 03193 }; 03194 }; 03195 03196 // alternative implementation without caching 03197 // 03198 // template <class TAG> 03199 // class Central 03200 // { 03201 // public: 03202 // typedef typename StandardizeTag<TAG>::type TargetTag; 03203 // typedef TypeList<Mean, typename TransferModifiers<Central<TargetTag>, typename TargetTag::Dependencies::type>::type> Dependencies; 03204 03205 // template <class U, class BASE> 03206 // struct Impl 03207 // : public TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE> 03208 // { 03209 // typedef typename TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE> ImplType; 03210 03211 // static const unsigned int workInPass = 2; 03212 03213 // void operator+=(Impl const & o) 03214 // { 03215 // vigra_precondition(false, 03216 // "Central<...>::operator+=(): not supported."); 03217 // } 03218 03219 // template <class T> 03220 // void update(T const & t) 03221 // { 03222 // ImplType::update(t - getDependency<Mean>(*this)); 03223 // } 03224 03225 // template <class T> 03226 // void update(T const & t, double weight) 03227 // { 03228 // ImplType::update(t - getDependency<Mean>(*this), weight); 03229 // } 03230 // }; 03231 // }; 03232 03233 03234 class PrincipalProjection 03235 { 03236 public: 03237 typedef Select<Centralize, Principal<CoordinateSystem> > Dependencies; 03238 03239 static std::string const & name() 03240 { 03241 static const std::string n("PrincipalProjection (internal)"); 03242 return n; 03243 } 03244 03245 template <class U, class BASE> 03246 struct Impl 03247 : public BASE 03248 { 03249 static const unsigned int workInPass = 2; 03250 03251 typedef typename AccumulatorResultTraits<U>::element_promote_type element_type; 03252 typedef typename AccumulatorResultTraits<U>::SumType value_type; 03253 typedef value_type const & result_type; 03254 03255 mutable value_type value_; 03256 03257 Impl() 03258 : value_() // call default constructor explicitly to ensure zero initialization 03259 {} 03260 03261 void reset() 03262 { 03263 value_ = element_type(); 03264 } 03265 03266 template <class Shape> 03267 void reshape(Shape const & s) 03268 { 03269 detail::reshapeImpl(value_, s); 03270 } 03271 03272 void update(U const & t) const 03273 { 03274 for(unsigned int k=0; k<t.size(); ++k) 03275 { 03276 value_[k] = getDependency<Principal<CoordinateSystem> >(*this)(0, k)*getDependency<Centralize>(*this)[0]; 03277 for(unsigned int d=1; d<t.size(); ++d) 03278 value_[k] += getDependency<Principal<CoordinateSystem> >(*this)(d, k)*getDependency<Centralize>(*this)[d]; 03279 } 03280 } 03281 03282 void update(U const & t, double) const 03283 { 03284 update(t); 03285 } 03286 03287 result_type operator()(U const & t) const 03288 { 03289 getAccumulator<Centralize>(*this).update(t); 03290 update(t); 03291 return value_; 03292 } 03293 03294 result_type operator()() const 03295 { 03296 return value_; 03297 } 03298 }; 03299 }; 03300 03301 /** \brief Modifier. Project onto PCA eigenvectors. 03302 03303 Works in pass 2, %operator+=() not supported (merging not supported). 03304 */ 03305 template <class TAG> 03306 class Principal 03307 { 03308 public: 03309 typedef typename StandardizeTag<TAG>::type TargetTag; 03310 typedef Select<PrincipalProjection, typename TargetTag::Dependencies> Dependencies; 03311 03312 static std::string const & name() 03313 { 03314 static const std::string n = std::string("Principal<") + TargetTag::name() + " >"; 03315 return n; 03316 } 03317 03318 template <class U, class BASE> 03319 struct Impl 03320 : public TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE> 03321 { 03322 typedef typename TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE> ImplType; 03323 03324 static const unsigned int workInPass = 2; 03325 03326 void operator+=(Impl const & o) 03327 { 03328 vigra_precondition(false, 03329 "Principal<...>::operator+=(): not supported."); 03330 } 03331 03332 template <class T> 03333 void update(T const & t) 03334 { 03335 ImplType::update(getDependency<PrincipalProjection>(*this)); 03336 } 03337 03338 template <class T> 03339 void update(T const & t, double weight) 03340 { 03341 ImplType::update(getDependency<PrincipalProjection>(*this), weight); 03342 } 03343 }; 03344 }; 03345 03346 /* 03347 important notes on modifiers: 03348 * upon accumulator creation, modifiers are reordered so that data preparation is innermost, 03349 and data access is outermost, e.g.: 03350 Coord<DivideByCount<Principal<PowerSum<2> > > > 03351 * modifiers are automatically transfered to dependencies as appropriate 03352 * modifiers for lookup (getAccumulator and get) of dependent accumulators are automatically adjusted 03353 * modifiers must adjust workInPass for the contained accumulator as appropriate 03354 * we may implement convenience versions of Select that apply a modifier to all 03355 contained tags at once 03356 * weighted accumulators have their own Count object when used together 03357 with unweighted ones (this is as yet untested - FIXME) 03358 * certain accumulators must remain unchanged when wrapped in certain modifiers: 03359 * Count: always except for Weighted<Count> and CoordWeighted<Count> 03360 * Sum: data preparation modifiers 03361 * FlatScatterMatrixImpl, CovarianceEigensystemImpl: Principal and Whitened 03362 * will it be useful to implement initPass<N>() or finalizePass<N>() ? 03363 */ 03364 03365 /****************************************************************************/ 03366 /* */ 03367 /* the actual accumulators */ 03368 /* */ 03369 /****************************************************************************/ 03370 03371 /** \brief Basic statistic. Identity matrix of appropriate size. 03372 */ 03373 class CoordinateSystem 03374 { 03375 public: 03376 typedef Select<> Dependencies; 03377 03378 static std::string const & name() 03379 { 03380 static const std::string n("CoordinateSystem"); 03381 return n; 03382 } 03383 03384 template <class U, class BASE> 03385 struct Impl 03386 : public BASE 03387 { 03388 typedef double element_type; 03389 typedef Matrix<double> value_type; 03390 typedef value_type const & result_type; 03391 03392 value_type value_; 03393 03394 Impl() 03395 : value_() // call default constructor explicitly to ensure zero initialization 03396 {} 03397 03398 void reset() 03399 { 03400 value_ = element_type(); 03401 } 03402 03403 template <class Shape> 03404 void reshape(Shape const & s) 03405 { 03406 detail::reshapeImpl(value_, s); 03407 } 03408 03409 result_type operator()() const 03410 { 03411 return value_; 03412 } 03413 }; 03414 }; 03415 03416 template <class BASE, class T, 03417 class ElementType=typename AccumulatorResultTraits<T>::element_promote_type, 03418 class SumType=typename AccumulatorResultTraits<T>::SumType> 03419 struct SumBaseImpl 03420 : public BASE 03421 { 03422 typedef ElementType element_type; 03423 typedef SumType value_type; 03424 typedef value_type const & result_type; 03425 03426 value_type value_; 03427 03428 SumBaseImpl() 03429 : value_() // call default constructor explicitly to ensure zero initialization 03430 {} 03431 03432 void reset() 03433 { 03434 value_ = element_type(); 03435 } 03436 03437 template <class Shape> 03438 void reshape(Shape const & s) 03439 { 03440 detail::reshapeImpl(value_, s); 03441 } 03442 03443 void operator+=(SumBaseImpl const & o) 03444 { 03445 value_ += o.value_; 03446 } 03447 03448 result_type operator()() const 03449 { 03450 return value_; 03451 } 03452 }; 03453 03454 // Count 03455 template <> 03456 class PowerSum<0> 03457 { 03458 public: 03459 typedef Select<> Dependencies; 03460 03461 static std::string const & name() 03462 { 03463 static const std::string n("PowerSum<0>"); 03464 return n; 03465 } 03466 03467 template <class T, class BASE> 03468 struct Impl 03469 : public SumBaseImpl<BASE, T, double, double> 03470 { 03471 void update(T const & t) 03472 { 03473 ++this->value_; 03474 } 03475 03476 void update(T const & t, double weight) 03477 { 03478 this->value_ += weight; 03479 } 03480 }; 03481 }; 03482 03483 // Sum 03484 template <> 03485 class PowerSum<1> 03486 { 03487 public: 03488 typedef Select<> Dependencies; 03489 03490 static std::string const & name() 03491 { 03492 static const std::string n("PowerSum<1>"); 03493 return n; 03494 } 03495 03496 template <class U, class BASE> 03497 struct Impl 03498 : public SumBaseImpl<BASE, U> 03499 { 03500 void update(U const & t) 03501 { 03502 this->value_ += t; 03503 } 03504 03505 void update(U const & t, double weight) 03506 { 03507 this->value_ += weight*t; 03508 } 03509 }; 03510 }; 03511 03512 /** \brief Basic statistic. PowerSum<N> =@f$ \sum_i x_i^N @f$ 03513 03514 Works in pass 1, %operator+=() supported (merging supported). 03515 */ 03516 template <unsigned N> 03517 class PowerSum 03518 { 03519 public: 03520 typedef Select<> Dependencies; 03521 03522 static std::string const & name() 03523 { 03524 static const std::string n = std::string("PowerSum<") + asString(N) + ">"; 03525 return n; 03526 } 03527 03528 template <class U, class BASE> 03529 struct Impl 03530 : public SumBaseImpl<BASE, U> 03531 { 03532 void update(U const & t) 03533 { 03534 using namespace vigra::multi_math; 03535 this->value_ += pow(t, (int)N); 03536 } 03537 03538 void update(U const & t, double weight) 03539 { 03540 using namespace vigra::multi_math; 03541 this->value_ += weight*pow(t, (int)N); 03542 } 03543 }; 03544 }; 03545 03546 template <> 03547 class AbsPowerSum<1> 03548 { 03549 public: 03550 typedef Select<> Dependencies; 03551 03552 static std::string const & name() 03553 { 03554 static const std::string n("AbsPowerSum<1>"); 03555 return n; 03556 } 03557 03558 template <class U, class BASE> 03559 struct Impl 03560 : public SumBaseImpl<BASE, U> 03561 { 03562 void update(U const & t) 03563 { 03564 using namespace vigra::multi_math; 03565 this->value_ += abs(t); 03566 } 03567 03568 void update(U const & t, double weight) 03569 { 03570 using namespace vigra::multi_math; 03571 this->value_ += weight*abs(t); 03572 } 03573 }; 03574 }; 03575 03576 /** \brief Basic statistic. AbsPowerSum<N> =@f$ \sum_i |x_i|^N @f$ 03577 03578 Works in pass 1, %operator+=() supported (merging supported). 03579 */ 03580 template <unsigned N> 03581 class AbsPowerSum 03582 { 03583 public: 03584 typedef Select<> Dependencies; 03585 03586 static std::string const & name() 03587 { 03588 static const std::string n = std::string("AbsPowerSum<") + asString(N) + ">"; 03589 return n; 03590 } 03591 03592 template <class U, class BASE> 03593 struct Impl 03594 : public SumBaseImpl<BASE, U> 03595 { 03596 void update(U const & t) 03597 { 03598 using namespace vigra::multi_math; 03599 this->value_ += pow(abs(t), (int)N); 03600 } 03601 03602 void update(U const & t, double weight) 03603 { 03604 using namespace vigra::multi_math; 03605 this->value_ += weight*pow(abs(t), (int)N); 03606 } 03607 }; 03608 }; 03609 03610 template <class BASE, class VALUE_TYPE, class U> 03611 struct CachedResultBase 03612 : public BASE 03613 { 03614 typedef typename AccumulatorResultTraits<U>::element_type element_type; 03615 typedef VALUE_TYPE value_type; 03616 typedef value_type const & result_type; 03617 03618 mutable value_type value_; 03619 03620 CachedResultBase() 03621 : value_() // call default constructor explicitly to ensure zero initialization 03622 {} 03623 03624 void reset() 03625 { 03626 value_ = element_type(); 03627 this->setClean(); 03628 } 03629 03630 template <class Shape> 03631 void reshape(Shape const & s) 03632 { 03633 detail::reshapeImpl(value_, s); 03634 } 03635 03636 void operator+=(CachedResultBase const &) 03637 { 03638 this->setDirty(); 03639 } 03640 03641 void update(U const &) 03642 { 03643 this->setDirty(); 03644 } 03645 03646 void update(U const &, double) 03647 { 03648 this->setDirty(); 03649 } 03650 }; 03651 03652 // cached Mean and Variance 03653 /** \brief Modifier. Divide statistic by Count: DivideByCount<TAG> = TAG / Count . 03654 */ 03655 template <class TAG> 03656 class DivideByCount 03657 { 03658 public: 03659 typedef typename StandardizeTag<TAG>::type TargetTag; 03660 typedef Select<TargetTag, Count> Dependencies; 03661 03662 static std::string const & name() 03663 { 03664 static const std::string n = std::string("DivideByCount<") + TargetTag::name() + " >"; 03665 return n; 03666 } 03667 03668 template <class U, class BASE> 03669 struct Impl 03670 : public CachedResultBase<BASE, typename LookupDependency<TargetTag, BASE>::value_type, U> 03671 { 03672 typedef typename CachedResultBase<BASE, typename LookupDependency<TargetTag, BASE>::value_type, U>::result_type result_type; 03673 03674 result_type operator()() const 03675 { 03676 if(this->isDirty()) 03677 { 03678 using namespace multi_math; 03679 this->value_ = getDependency<TargetTag>(*this) / getDependency<Count>(*this); 03680 this->setClean(); 03681 } 03682 return this->value_; 03683 } 03684 }; 03685 }; 03686 03687 // UnbiasedVariance 03688 /** \brief Modifier. Divide statistics by Count-1: DivideUnbiased<TAG> = TAG / (Count-1) 03689 */ 03690 template <class TAG> 03691 class DivideUnbiased 03692 { 03693 public: 03694 typedef typename StandardizeTag<TAG>::type TargetTag; 03695 typedef Select<TargetTag, Count> Dependencies; 03696 03697 static std::string const & name() 03698 { 03699 static const std::string n = std::string("DivideUnbiased<") + TargetTag::name() + " >"; 03700 return n; 03701 } 03702 03703 template <class U, class BASE> 03704 struct Impl 03705 : public BASE 03706 { 03707 typedef typename LookupDependency<TargetTag, BASE>::value_type value_type; 03708 typedef value_type result_type; 03709 03710 result_type operator()() const 03711 { 03712 using namespace multi_math; 03713 return getDependency<TargetTag>(*this) / (getDependency<Count>(*this) - 1.0); 03714 } 03715 }; 03716 }; 03717 03718 // RootMeanSquares and StdDev 03719 /** \brief Modifier. RootDivideByCount<TAG> = sqrt( TAG/Count ) 03720 */ 03721 template <class TAG> 03722 class RootDivideByCount 03723 { 03724 public: 03725 typedef typename StandardizeTag<DivideByCount<TAG> >::type TargetTag; 03726 typedef Select<TargetTag> Dependencies; 03727 03728 static std::string const & name() 03729 { 03730 typedef typename StandardizeTag<TAG>::type InnerTag; 03731 static const std::string n = std::string("RootDivideByCount<") + InnerTag::name() + " >"; 03732 return n; 03733 } 03734 03735 template <class U, class BASE> 03736 struct Impl 03737 : public BASE 03738 { 03739 typedef typename LookupDependency<TargetTag, BASE>::value_type value_type; 03740 typedef value_type result_type; 03741 03742 result_type operator()() const 03743 { 03744 using namespace multi_math; 03745 return sqrt(getDependency<TargetTag>(*this)); 03746 } 03747 }; 03748 }; 03749 03750 // UnbiasedStdDev 03751 /** \brief Modifier. RootDivideUnbiased<TAG> = sqrt( TAG / (Count-1) ) 03752 */ 03753 template <class TAG> 03754 class RootDivideUnbiased 03755 { 03756 public: 03757 typedef typename StandardizeTag<DivideUnbiased<TAG> >::type TargetTag; 03758 typedef Select<TargetTag> Dependencies; 03759 03760 static std::string const & name() 03761 { 03762 typedef typename StandardizeTag<TAG>::type InnerTag; 03763 static const std::string n = std::string("RootDivideUnbiased<") + InnerTag::name() + " >"; 03764 return n; 03765 } 03766 03767 template <class U, class BASE> 03768 struct Impl 03769 : public BASE 03770 { 03771 typedef typename LookupDependency<TargetTag, BASE>::value_type value_type; 03772 typedef value_type result_type; 03773 03774 result_type operator()() const 03775 { 03776 using namespace multi_math; 03777 return sqrt(getDependency<TargetTag>(*this)); 03778 } 03779 }; 03780 }; 03781 03782 /** \brief Spezialization: works in pass 1, %operator+=() supported (merging supported). 03783 */ 03784 template <> 03785 class Central<PowerSum<2> > 03786 { 03787 public: 03788 typedef Select<Mean, Count> Dependencies; 03789 03790 static std::string const & name() 03791 { 03792 static const std::string n("Central<PowerSum<2> >"); 03793 return n; 03794 } 03795 03796 template <class U, class BASE> 03797 struct Impl 03798 : public SumBaseImpl<BASE, U> 03799 { 03800 void operator+=(Impl const & o) 03801 { 03802 using namespace vigra::multi_math; 03803 double n1 = getDependency<Count>(*this), n2 = getDependency<Count>(o); 03804 if(n1 == 0.0) 03805 { 03806 this->value_ = o.value_; 03807 } 03808 else if(n2 != 0.0) 03809 { 03810 this->value_ += o.value_ + n1 * n2 / (n1 + n2) * sq(getDependency<Mean>(*this) - getDependency<Mean>(o)); 03811 } 03812 } 03813 03814 void update(U const & t) 03815 { 03816 double n = getDependency<Count>(*this); 03817 if(n > 1.0) 03818 { 03819 using namespace vigra::multi_math; 03820 this->value_ += n / (n - 1.0) * sq(getDependency<Mean>(*this) - t); 03821 } 03822 } 03823 03824 void update(U const & t, double weight) 03825 { 03826 double n = getDependency<Count>(*this); 03827 if(n > weight) 03828 { 03829 using namespace vigra::multi_math; 03830 this->value_ += n / (n - weight) * sq(getDependency<Mean>(*this) - t); 03831 } 03832 } 03833 }; 03834 }; 03835 03836 /** \brief Specialization: works in pass 2, %operator+=() supported (merging supported). 03837 */ 03838 template <> 03839 class Central<PowerSum<3> > 03840 { 03841 public: 03842 typedef Select<Centralize, Count, Mean, Central<PowerSum<2> > > Dependencies; 03843 03844 static std::string const & name() 03845 { 03846 static const std::string n("Central<PowerSum<3> >"); 03847 return n; 03848 } 03849 03850 template <class U, class BASE> 03851 struct Impl 03852 : public SumBaseImpl<BASE, U> 03853 { 03854 typedef typename SumBaseImpl<BASE, U>::value_type value_type; 03855 03856 static const unsigned int workInPass = 2; 03857 03858 void operator+=(Impl const & o) 03859 { 03860 typedef Central<PowerSum<2> > Sum2Tag; 03861 03862 using namespace vigra::multi_math; 03863 double n1 = getDependency<Count>(*this), n2 = getDependency<Count>(o); 03864 if(n1 == 0.0) 03865 { 03866 this->value_ = o.value_; 03867 } 03868 else if(n2 != 0.0) 03869 { 03870 double n = n1 + n2; 03871 double weight = n1 * n2 * (n1 - n2) / sq(n); 03872 value_type delta = getDependency<Mean>(o) - getDependency<Mean>(*this); 03873 this->value_ += o.value_ + weight * pow(delta, 3) + 03874 3.0 / n * delta * (n1 * getDependency<Sum2Tag>(o) - n2 * getDependency<Sum2Tag>(*this)); 03875 } 03876 } 03877 03878 void update(U const & t) 03879 { 03880 using namespace vigra::multi_math; 03881 this->value_ += pow(getDependency<Centralize>(*this), 3); 03882 } 03883 03884 void update(U const & t, double weight) 03885 { 03886 using namespace vigra::multi_math; 03887 this->value_ += weight*pow(getDependency<Centralize>(*this), 3); 03888 } 03889 }; 03890 }; 03891 /** \brief Specialization: works in pass 2, %operator+=() supported (merging supported). 03892 */ 03893 template <> 03894 class Central<PowerSum<4> > 03895 { 03896 public: 03897 typedef Select<Centralize, Central<PowerSum<3> > > Dependencies; 03898 03899 static std::string const & name() 03900 { 03901 static const std::string n("Central<PowerSum<4> >"); 03902 return n; 03903 } 03904 03905 template <class U, class BASE> 03906 struct Impl 03907 : public SumBaseImpl<BASE, U> 03908 { 03909 typedef typename SumBaseImpl<BASE, U>::value_type value_type; 03910 03911 static const unsigned int workInPass = 2; 03912 03913 void operator+=(Impl const & o) 03914 { 03915 typedef Central<PowerSum<2> > Sum2Tag; 03916 typedef Central<PowerSum<3> > Sum3Tag; 03917 03918 using namespace vigra::multi_math; 03919 double n1 = getDependency<Count>(*this), n2 = getDependency<Count>(o); 03920 if(n1 == 0.0) 03921 { 03922 this->value_ = o.value_; 03923 } 03924 else if(n2 != 0.0) 03925 { 03926 double n = n1 + n2; 03927 double n1_2 = sq(n1); 03928 double n2_2 = sq(n2); 03929 double n_2 = sq(n); 03930 double weight = n1 * n2 * (n1_2 - n1*n2 + n2_2) / n_2 / n; 03931 value_type delta = getDependency<Mean>(o) - getDependency<Mean>(*this); 03932 this->value_ += o.value_ + weight * pow(delta, 4) + 03933 6.0 / n_2 * sq(delta) * (n1_2 * getDependency<Sum2Tag>(o) + n2_2 * getDependency<Sum2Tag>(*this)) + 03934 4.0 / n * delta * (n1 * getDependency<Sum3Tag>(o) - n2 * getDependency<Sum3Tag>(*this)); 03935 } 03936 } 03937 03938 void update(U const & t) 03939 { 03940 using namespace vigra::multi_math; 03941 this->value_ += pow(getDependency<Centralize>(*this), 4); 03942 } 03943 03944 void update(U const & t, double weight) 03945 { 03946 using namespace vigra::multi_math; 03947 this->value_ += weight*pow(getDependency<Centralize>(*this), 4); 03948 } 03949 }; 03950 }; 03951 03952 /** \brief Basic statistic. Skewness. 03953 03954 %Skewness =@f$ \frac{ \frac{1}{n}\sum_i (x_i-\hat{x})^3 }{ (\frac{1}{n}\sum_i (x_i-\hat{x})^2)^{3/2} } @f$ . 03955 Works in pass 2, %operator+=() supported (merging supported). 03956 */ 03957 class Skewness 03958 { 03959 public: 03960 typedef Select<Central<PowerSum<2> >, Central<PowerSum<3> > > Dependencies; 03961 03962 static std::string const & name() 03963 { 03964 static const std::string n("Skewness"); 03965 return n; 03966 } 03967 03968 template <class U, class BASE> 03969 struct Impl 03970 : public BASE 03971 { 03972 static const unsigned int workInPass = 2; 03973 03974 typedef typename LookupDependency<Central<PowerSum<3> >, BASE>::value_type value_type; 03975 typedef value_type result_type; 03976 03977 result_type operator()() const 03978 { 03979 typedef Central<PowerSum<3> > Sum3; 03980 typedef Central<PowerSum<2> > Sum2; 03981 03982 using namespace multi_math; 03983 return sqrt(getDependency<Count>(*this)) * getDependency<Sum3>(*this) / pow(getDependency<Sum2>(*this), 1.5); 03984 } 03985 }; 03986 }; 03987 03988 /** \brief Basic statistic. Unbiased Skewness. 03989 03990 Works in pass 2, %operator+=() supported (merging supported). 03991 */ 03992 class UnbiasedSkewness 03993 { 03994 public: 03995 typedef Select<Skewness> Dependencies; 03996 03997 static std::string const & name() 03998 { 03999 static const std::string n("UnbiasedSkewness"); 04000 return n; 04001 } 04002 04003 template <class U, class BASE> 04004 struct Impl 04005 : public BASE 04006 { 04007 static const unsigned int workInPass = 2; 04008 04009 typedef typename LookupDependency<Central<PowerSum<3> >, BASE>::value_type value_type; 04010 typedef value_type result_type; 04011 04012 result_type operator()() const 04013 { 04014 using namespace multi_math; 04015 double n = getDependency<Count>(*this); 04016 return sqrt(n*(n-1.0)) / (n - 2.0) * getDependency<Skewness>(*this); 04017 } 04018 }; 04019 }; 04020 04021 /** \brief Basic statistic. Kurtosis. 04022 04023 %Kurtosis = @f$ \frac{ \frac{1}{n}\sum_i (x_i-\bar{x})^4 }{ 04024 (\frac{1}{n} \sum_i(x_i-\bar{x})^2)^2 } - 3 @f$ . 04025 Works in pass 2, %operator+=() supported (merging supported). 04026 */ 04027 class Kurtosis 04028 { 04029 public: 04030 typedef Select<Central<PowerSum<2> >, Central<PowerSum<4> > > Dependencies; 04031 04032 static std::string const & name() 04033 { 04034 static const std::string n("Kurtosis"); 04035 return n; 04036 } 04037 04038 template <class U, class BASE> 04039 struct Impl 04040 : public BASE 04041 { 04042 static const unsigned int workInPass = 2; 04043 04044 typedef typename LookupDependency<Central<PowerSum<4> >, BASE>::value_type value_type; 04045 typedef value_type result_type; 04046 04047 result_type operator()() const 04048 { 04049 typedef Central<PowerSum<4> > Sum4; 04050 typedef Central<PowerSum<2> > Sum2; 04051 04052 using namespace multi_math; 04053 return getDependency<Count>(*this) * getDependency<Sum4>(*this) / sq(getDependency<Sum2>(*this)) - value_type(3.0); 04054 } 04055 }; 04056 }; 04057 04058 /** \brief Basic statistic. Unbiased Kurtosis. 04059 04060 Works in pass 2, %operator+=() supported (merging supported). 04061 */ 04062 class UnbiasedKurtosis 04063 { 04064 public: 04065 typedef Select<Kurtosis> Dependencies; 04066 04067 static std::string const & name() 04068 { 04069 static const std::string n("UnbiasedKurtosis"); 04070 return n; 04071 } 04072 04073 template <class U, class BASE> 04074 struct Impl 04075 : public BASE 04076 { 04077 static const unsigned int workInPass = 2; 04078 04079 typedef typename LookupDependency<Central<PowerSum<4> >, BASE>::value_type value_type; 04080 typedef value_type result_type; 04081 04082 result_type operator()() const 04083 { 04084 using namespace multi_math; 04085 double n = getDependency<Count>(*this); 04086 return (n-1.0)/((n-2.0)*(n-3.0))*((n+1.0)*getDependency<Kurtosis>(*this) + value_type(6.0)); 04087 } 04088 }; 04089 }; 04090 04091 namespace detail { 04092 04093 template <class Scatter, class Sum> 04094 void updateFlatScatterMatrix(Scatter & sc, Sum const & s, double w) 04095 { 04096 int size = s.size(); 04097 for(MultiArrayIndex j=0, k=0; j<size; ++j) 04098 for(MultiArrayIndex i=j; i<size; ++i, ++k) 04099 sc[k] += w*s[i]*s[j]; 04100 } 04101 04102 template <class Sum> 04103 void updateFlatScatterMatrix(double & sc, Sum const & s, double w) 04104 { 04105 sc += w*s*s; 04106 } 04107 04108 template <class Cov, class Scatter> 04109 void flatScatterMatrixToScatterMatrix(Cov & cov, Scatter const & sc) 04110 { 04111 int size = cov.shape(0), k=0; 04112 for(MultiArrayIndex j=0; j<size; ++j) 04113 { 04114 cov(j,j) = sc[k++]; 04115 for(MultiArrayIndex i=j+1; i<size; ++i) 04116 { 04117 cov(i,j) = sc[k++]; 04118 cov(j,i) = cov(i,j); 04119 } 04120 } 04121 } 04122 04123 template <class Scatter> 04124 void flatScatterMatrixToScatterMatrix(double & cov, Scatter const & sc) 04125 { 04126 cov = sc; 04127 } 04128 04129 template <class Cov, class Scatter> 04130 void flatScatterMatrixToCovariance(Cov & cov, Scatter const & sc, double n) 04131 { 04132 int size = cov.shape(0), k=0; 04133 for(MultiArrayIndex j=0; j<size; ++j) 04134 { 04135 cov(j,j) = sc[k++] / n; 04136 for(MultiArrayIndex i=j+1; i<size; ++i) 04137 { 04138 cov(i,j) = sc[k++] / n; 04139 cov(j,i) = cov(i,j); 04140 } 04141 } 04142 } 04143 04144 template <class Scatter> 04145 void flatScatterMatrixToCovariance(double & cov, Scatter const & sc, double n) 04146 { 04147 cov = sc / n; 04148 } 04149 04150 } // namespace detail 04151 04152 // we only store the flattened upper triangular part of the scatter matrix 04153 /** \brief Basic statistic. Flattened uppter-triangular part of scatter matrix. 04154 04155 Works in pass 1, %operator+=() supported (merging supported). 04156 */ 04157 class FlatScatterMatrix 04158 { 04159 public: 04160 typedef Select<Mean, Count> Dependencies; 04161 04162 static std::string const & name() 04163 { 04164 static const std::string n("FlatScatterMatrix"); 04165 return n; 04166 } 04167 04168 template <class U, class BASE> 04169 struct Impl 04170 : public BASE 04171 { 04172 typedef typename AccumulatorResultTraits<U>::element_promote_type element_type; 04173 typedef typename AccumulatorResultTraits<U>::FlatCovarianceType value_type; 04174 typedef value_type const & result_type; 04175 04176 typedef typename AccumulatorResultTraits<U>::SumType SumType; 04177 04178 value_type value_; 04179 SumType diff_; 04180 04181 Impl() 04182 : value_(), // call default constructor explicitly to ensure zero initialization 04183 diff_() 04184 {} 04185 04186 void reset() 04187 { 04188 value_ = element_type(); 04189 } 04190 04191 template <class Shape> 04192 void reshape(Shape const & s) 04193 { 04194 int size = prod(s); 04195 detail::reshapeImpl(value_, Shape1(size*(size+1)/2)); 04196 detail::reshapeImpl(diff_, s); 04197 } 04198 04199 void operator+=(Impl const & o) 04200 { 04201 double n1 = getDependency<Count>(*this), n2 = getDependency<Count>(o); 04202 if(n1 == 0.0) 04203 { 04204 value_ = o.value_; 04205 } 04206 else if(n2 != 0.0) 04207 { 04208 using namespace vigra::multi_math; 04209 diff_ = getDependency<Mean>(*this) - getDependency<Mean>(o); 04210 detail::updateFlatScatterMatrix(value_, diff_, n1 * n2 / (n1 + n2)); 04211 value_ += o.value_; 04212 } 04213 } 04214 04215 void update(U const & t) 04216 { 04217 compute(t); 04218 } 04219 04220 void update(U const & t, double weight) 04221 { 04222 compute(t, weight); 04223 } 04224 04225 result_type operator()() const 04226 { 04227 return value_; 04228 } 04229 04230 private: 04231 void compute(U const & t, double weight = 1.0) 04232 { 04233 double n = getDependency<Count>(*this); 04234 if(n > weight) 04235 { 04236 using namespace vigra::multi_math; 04237 diff_ = getDependency<Mean>(*this) - t; 04238 detail::updateFlatScatterMatrix(value_, diff_, n * weight / (n - weight)); 04239 } 04240 } 04241 }; 04242 }; 04243 04244 // Covariance 04245 template <> 04246 class DivideByCount<FlatScatterMatrix> 04247 { 04248 public: 04249 typedef Select<FlatScatterMatrix, Count> Dependencies; 04250 04251 static std::string const & name() 04252 { 04253 static const std::string n("DivideByCount<FlatScatterMatrix>"); 04254 return n; 04255 } 04256 04257 template <class U, class BASE> 04258 struct Impl 04259 : public CachedResultBase<BASE, typename AccumulatorResultTraits<U>::CovarianceType, U> 04260 { 04261 typedef CachedResultBase<BASE, typename AccumulatorResultTraits<U>::CovarianceType, U> BaseType; 04262 typedef typename BaseType::result_type result_type; 04263 04264 template <class Shape> 04265 void reshape(Shape const & s) 04266 { 04267 int size = prod(s); 04268 detail::reshapeImpl(this->value_, Shape2(size,size)); 04269 } 04270 04271 result_type operator()() const 04272 { 04273 if(this->isDirty()) 04274 { 04275 detail::flatScatterMatrixToCovariance(this->value_, getDependency<FlatScatterMatrix>(*this), getDependency<Count>(*this)); 04276 this->setClean(); 04277 } 04278 return this->value_; 04279 } 04280 }; 04281 }; 04282 04283 // UnbiasedCovariance 04284 template <> 04285 class DivideUnbiased<FlatScatterMatrix> 04286 { 04287 public: 04288 typedef Select<FlatScatterMatrix, Count> Dependencies; 04289 04290 static std::string const & name() 04291 { 04292 static const std::string n("DivideUnbiased<FlatScatterMatrix>"); 04293 return n; 04294 } 04295 04296 template <class U, class BASE> 04297 struct Impl 04298 : public CachedResultBase<BASE, typename AccumulatorResultTraits<U>::CovarianceType, U> 04299 { 04300 typedef CachedResultBase<BASE, typename AccumulatorResultTraits<U>::CovarianceType, U> BaseType; 04301 typedef typename BaseType::result_type result_type; 04302 04303 template <class Shape> 04304 void reshape(Shape const & s) 04305 { 04306 int size = prod(s); 04307 detail::reshapeImpl(this->value_, Shape2(size,size)); 04308 } 04309 04310 result_type operator()() const 04311 { 04312 if(this->isDirty()) 04313 { 04314 detail::flatScatterMatrixToCovariance(this->value_, getDependency<FlatScatterMatrix>(*this), getDependency<Count>(*this) - 1.0); 04315 this->setClean(); 04316 } 04317 return this->value_; 04318 } 04319 }; 04320 }; 04321 04322 /** Basic statistic. ScatterMatrixEigensystem (?) 04323 */ 04324 class ScatterMatrixEigensystem 04325 { 04326 public: 04327 typedef Select<FlatScatterMatrix> Dependencies; 04328 04329 static std::string const & name() 04330 { 04331 static const std::string n("ScatterMatrixEigensystem"); 04332 return n; 04333 } 04334 04335 template <class U, class BASE> 04336 struct Impl 04337 : public BASE 04338 { 04339 typedef typename AccumulatorResultTraits<U>::element_promote_type element_type; 04340 typedef typename AccumulatorResultTraits<U>::SumType EigenvalueType; 04341 typedef typename AccumulatorResultTraits<U>::CovarianceType EigenvectorType; 04342 typedef std::pair<EigenvalueType, EigenvectorType> value_type; 04343 typedef value_type const & result_type; 04344 04345 mutable value_type value_; 04346 04347 Impl() 04348 : value_() 04349 {} 04350 04351 void operator+=(Impl const &) 04352 { 04353 this->setDirty(); 04354 } 04355 04356 void update(U const &) 04357 { 04358 this->setDirty(); 04359 } 04360 04361 void update(U const &, double) 04362 { 04363 this->setDirty(); 04364 } 04365 04366 void reset() 04367 { 04368 value_.first = element_type(); 04369 value_.second = element_type(); 04370 this->setClean(); 04371 } 04372 04373 template <class Shape> 04374 void reshape(Shape const & s) 04375 { 04376 int size = prod(s); 04377 detail::reshapeImpl(value_.first, Shape1(size)); 04378 detail::reshapeImpl(value_.second, Shape2(size,size)); 04379 } 04380 04381 result_type operator()() const 04382 { 04383 if(this->isDirty()) 04384 { 04385 compute(getDependency<FlatScatterMatrix>(*this), value_.first, value_.second); 04386 this->setClean(); 04387 } 04388 return value_; 04389 } 04390 04391 private: 04392 template <class Flat, class EW, class EV> 04393 static void compute(Flat const & flatScatter, EW & ew, EV & ev) 04394 { 04395 EigenvectorType scatter(ev.shape()); 04396 detail::flatScatterMatrixToScatterMatrix(scatter, flatScatter); 04397 // create a view because EW could be a TinyVector 04398 MultiArrayView<2, element_type> ewview(Shape2(ev.shape(0), 1), &ew[0]); 04399 symmetricEigensystem(scatter, ewview, ev); 04400 } 04401 04402 static void compute(double v, double & ew, double & ev) 04403 { 04404 ew = v; 04405 ev = 1.0; 04406 } 04407 }; 04408 }; 04409 04410 // CovarianceEigensystem 04411 template <> 04412 class DivideByCount<ScatterMatrixEigensystem> 04413 { 04414 public: 04415 typedef Select<ScatterMatrixEigensystem, Count> Dependencies; 04416 04417 static std::string const & name() 04418 { 04419 static const std::string n("DivideByCount<ScatterMatrixEigensystem>"); 04420 return n; 04421 } 04422 04423 template <class U, class BASE> 04424 struct Impl 04425 : public BASE 04426 { 04427 typedef typename LookupDependency<ScatterMatrixEigensystem, BASE>::type SMImpl; 04428 typedef typename SMImpl::element_type element_type; 04429 typedef typename SMImpl::EigenvalueType EigenvalueType; 04430 typedef typename SMImpl::EigenvectorType EigenvectorType; 04431 typedef std::pair<EigenvalueType, EigenvectorType const &> value_type; 04432 typedef value_type const & result_type; 04433 04434 mutable value_type value_; 04435 04436 Impl() 04437 : value_(EigenvalueType(), BASE::template call_getDependency<ScatterMatrixEigensystem>().second) 04438 {} 04439 04440 void operator+=(Impl const &) 04441 { 04442 this->setDirty(); 04443 } 04444 04445 void update(U const &) 04446 { 04447 this->setDirty(); 04448 } 04449 04450 void update(U const &, double) 04451 { 04452 this->setDirty(); 04453 } 04454 04455 void reset() 04456 { 04457 value_.first = element_type(); 04458 this->setClean(); 04459 } 04460 04461 template <class Shape> 04462 void reshape(Shape const & s) 04463 { 04464 int size = prod(s); 04465 detail::reshapeImpl(value_.first, Shape2(size,1)); 04466 } 04467 04468 result_type operator()() const 04469 { 04470 if(this->isDirty()) 04471 { 04472 value_.first = getDependency<ScatterMatrixEigensystem>(*this).first / getDependency<Count>(*this); 04473 this->setClean(); 04474 } 04475 return value_; 04476 } 04477 }; 04478 }; 04479 04480 // alternative implementation of CovarianceEigensystem - solve eigensystem directly 04481 // 04482 // template <> 04483 // class DivideByCount<ScatterMatrixEigensystem> 04484 // { 04485 // public: 04486 // typedef Select<Covariance> Dependencies; 04487 04488 // template <class U, class BASE> 04489 // struct Impl 04490 // : public BASE 04491 // { 04492 // typedef typename AccumulatorResultTraits<U>::element_promote_type element_type; 04493 // typedef typename AccumulatorResultTraits<U>::SumType EigenvalueType; 04494 // typedef typename AccumulatorResultTraits<U>::CovarianceType EigenvectorType; 04495 // typedef std::pair<EigenvalueType, EigenvectorType> value_type; 04496 // typedef value_type const & result_type; 04497 04498 // mutable value_type value_; 04499 04500 // Impl() 04501 // : value_() 04502 // {} 04503 04504 // void operator+=(Impl const &) 04505 // { 04506 // this->setDirty(); 04507 // } 04508 04509 // void update(U const &) 04510 // { 04511 // this->setDirty(); 04512 // } 04513 04514 // void update(U const &, double) 04515 // { 04516 // this->setDirty(); 04517 // } 04518 04519 // void reset() 04520 // { 04521 // value_.first = element_type(); 04522 // value_.second = element_type(); 04523 // this->setClean(); 04524 // } 04525 04526 // template <class Shape> 04527 // void reshape(Shape const & s) 04528 // { 04529 // int size = prod(s); 04530 // detail::reshapeImpl(value_.first, Shape2(size,1)); 04531 // detail::reshapeImpl(value_.second, Shape2(size,size)); 04532 // } 04533 04534 // result_type operator()() const 04535 // { 04536 // if(this->isDirty()) 04537 // { 04538 // compute(getDependency<Covariance>(*this), value_.first, value_.second); 04539 // this->setClean(); 04540 // } 04541 // return value_; 04542 // } 04543 04544 // private: 04545 // template <class Cov, class EW, class EV> 04546 // static void compute(Cov const & cov, EW & ew, EV & ev) 04547 // { 04548 // // create a view because EW could be a TinyVector 04549 // MultiArrayView<2, element_type> ewview(Shape2(cov.shape(0), 1), &ew[0]); 04550 // symmetricEigensystem(cov, ewview, ev); 04551 // } 04552 04553 // static void compute(double cov, double & ew, double & ev) 04554 // { 04555 // ew = cov; 04556 // ev = 1.0; 04557 // } 04558 // }; 04559 // }; 04560 04561 // covariance eigenvalues 04562 /** \brief Specialization (covariance eigenvalues): works in pass 1, %operator+=() supported (merging). 04563 */ 04564 template <> 04565 class Principal<PowerSum<2> > 04566 { 04567 public: 04568 typedef Select<ScatterMatrixEigensystem> Dependencies; 04569 04570 static std::string const & name() 04571 { 04572 static const std::string n("Principal<PowerSum<2> >"); 04573 return n; 04574 } 04575 04576 template <class U, class BASE> 04577 struct Impl 04578 : public BASE 04579 { 04580 typedef typename LookupDependency<ScatterMatrixEigensystem, BASE>::type::EigenvalueType value_type; 04581 typedef value_type const & result_type; 04582 04583 result_type operator()() const 04584 { 04585 return getDependency<ScatterMatrixEigensystem>(*this).first; 04586 } 04587 }; 04588 }; 04589 04590 04591 // Principal<CoordinateSystem> == covariance eigenvectors 04592 /** \brief Specialization (covariance eigenvectors): works in pass 1, %operator+=() supported (merging). 04593 */ 04594 template <> 04595 class Principal<CoordinateSystem> 04596 { 04597 public: 04598 typedef Select<ScatterMatrixEigensystem> Dependencies; 04599 04600 static std::string const & name() 04601 { 04602 static const std::string n("Principal<CoordinateSystem>"); 04603 return n; 04604 } 04605 04606 template <class U, class BASE> 04607 struct Impl 04608 : public BASE 04609 { 04610 typedef typename LookupDependency<ScatterMatrixEigensystem, BASE>::type::EigenvectorType value_type; 04611 typedef value_type const & result_type; 04612 04613 result_type operator()() const 04614 { 04615 return getDependency<ScatterMatrixEigensystem>(*this).second; 04616 } 04617 }; 04618 }; 04619 04620 /** \brief Basic statistic. %Minimum value. 04621 04622 Works in pass 1, %operator+=() supported (merging supported). 04623 */ 04624 class Minimum 04625 { 04626 public: 04627 typedef Select<> Dependencies; 04628 04629 static std::string const & name() 04630 { 04631 static const std::string n("Minimum"); 04632 return n; 04633 } 04634 04635 template <class U, class BASE> 04636 struct Impl 04637 : public BASE 04638 { 04639 typedef typename AccumulatorResultTraits<U>::element_type element_type; 04640 typedef typename AccumulatorResultTraits<U>::MinmaxType value_type; 04641 typedef value_type const & result_type; 04642 04643 value_type value_; 04644 04645 Impl() 04646 { 04647 value_ = NumericTraits<element_type>::max(); 04648 } 04649 04650 void reset() 04651 { 04652 value_ = NumericTraits<element_type>::max(); 04653 } 04654 04655 template <class Shape> 04656 void reshape(Shape const & s) 04657 { 04658 detail::reshapeImpl(value_, s, NumericTraits<element_type>::max()); 04659 } 04660 04661 void operator+=(Impl const & o) 04662 { 04663 updateImpl(o.value_); // necessary because std::min causes ambiguous overload 04664 } 04665 04666 void update(U const & t) 04667 { 04668 updateImpl(t); 04669 } 04670 04671 void update(U const & t, double) 04672 { 04673 updateImpl(t); 04674 } 04675 04676 result_type operator()() const 04677 { 04678 return value_; 04679 } 04680 04681 private: 04682 template <class T> 04683 void updateImpl(T const & o) 04684 { 04685 using namespace multi_math; 04686 value_ = min(value_, o); 04687 } 04688 04689 template <class T, class Alloc> 04690 void updateImpl(MultiArray<1, T, Alloc> const & o) 04691 { 04692 value_ = multi_math::min(value_, o); 04693 } 04694 }; 04695 }; 04696 04697 /** \brief Basic statistic. %Maximum value. 04698 04699 Works in pass 1, %operator+=() supported (merging supported). 04700 */ 04701 class Maximum 04702 { 04703 public: 04704 typedef Select<> Dependencies; 04705 04706 static std::string const & name() 04707 { 04708 static const std::string n("Maximum"); 04709 return n; 04710 } 04711 04712 template <class U, class BASE> 04713 struct Impl 04714 : public BASE 04715 { 04716 typedef typename AccumulatorResultTraits<U>::element_type element_type; 04717 typedef typename AccumulatorResultTraits<U>::MinmaxType value_type; 04718 typedef value_type const & result_type; 04719 04720 value_type value_; 04721 04722 Impl() 04723 { 04724 value_ = NumericTraits<element_type>::min(); 04725 } 04726 04727 void reset() 04728 { 04729 value_ = NumericTraits<element_type>::min(); 04730 } 04731 04732 template <class Shape> 04733 void reshape(Shape const & s) 04734 { 04735 detail::reshapeImpl(value_, s, NumericTraits<element_type>::min()); 04736 } 04737 04738 void operator+=(Impl const & o) 04739 { 04740 updateImpl(o.value_); // necessary because std::max causes ambiguous overload 04741 } 04742 04743 void update(U const & t) 04744 { 04745 updateImpl(t); 04746 } 04747 04748 void update(U const & t, double) 04749 { 04750 updateImpl(t); 04751 } 04752 04753 result_type operator()() const 04754 { 04755 return value_; 04756 } 04757 04758 private: 04759 template <class T> 04760 void updateImpl(T const & o) 04761 { 04762 using namespace multi_math; 04763 value_ = max(value_, o); 04764 } 04765 04766 template <class T, class Alloc> 04767 void updateImpl(MultiArray<1, T, Alloc> const & o) 04768 { 04769 value_ = multi_math::max(value_, o); 04770 } 04771 }; 04772 }; 04773 04774 /** \brief Basic statistic. Data value where weight assumes its minimal value. 04775 04776 Weights must be given. Coord<ArgMinWeight> gives coordinate where weight assumes its minimal value. Works in pass 1, %operator+=() supported (merging supported). 04777 */ 04778 class ArgMinWeight 04779 { 04780 public: 04781 typedef Select<> Dependencies; 04782 04783 static std::string const & name() 04784 { 04785 static const std::string n("ArgMinWeight"); 04786 return n; 04787 } 04788 04789 template <class U, class BASE> 04790 struct Impl 04791 : public BASE 04792 { 04793 typedef typename AccumulatorResultTraits<U>::element_type element_type; 04794 typedef typename AccumulatorResultTraits<U>::MinmaxType value_type; 04795 typedef value_type const & result_type; 04796 04797 double min_weight_; 04798 value_type value_; 04799 04800 Impl() 04801 : min_weight_(NumericTraits<double>::max()), 04802 value_() 04803 {} 04804 04805 void reset() 04806 { 04807 min_weight_ = NumericTraits<double>::max(); 04808 value_ = element_type(); 04809 } 04810 04811 template <class Shape> 04812 void reshape(Shape const & s) 04813 { 04814 detail::reshapeImpl(value_, s); 04815 } 04816 04817 void operator+=(Impl const & o) 04818 { 04819 using namespace multi_math; 04820 if(o.min_weight_ < min_weight_) 04821 { 04822 min_weight_ = o.min_weight_; 04823 value_ = o.value_; 04824 } 04825 } 04826 04827 void update(U const & t) 04828 { 04829 vigra_precondition(false, "ArgMinWeight::update() needs weights."); 04830 } 04831 04832 void update(U const & t, double weight) 04833 { 04834 if(weight < min_weight_) 04835 { 04836 min_weight_ = weight; 04837 value_ = t; 04838 } 04839 } 04840 04841 result_type operator()() const 04842 { 04843 return value_; 04844 } 04845 }; 04846 }; 04847 04848 /** \brief Basic statistic. Data where weight assumes its maximal value. 04849 04850 Weights must be given. Coord<ArgMinWeight> gives coordinate where weight assumes its maximal value. Works in pass 1, %operator+=() supported (merging supported). 04851 */ 04852 class ArgMaxWeight 04853 { 04854 public: 04855 typedef Select<> Dependencies; 04856 04857 static std::string const & name() 04858 { 04859 static const std::string n("ArgMaxWeight"); 04860 return n; 04861 } 04862 04863 template <class U, class BASE> 04864 struct Impl 04865 : public BASE 04866 { 04867 typedef typename AccumulatorResultTraits<U>::element_type element_type; 04868 typedef typename AccumulatorResultTraits<U>::MinmaxType value_type; 04869 typedef value_type const & result_type; 04870 04871 double max_weight_; 04872 value_type value_; 04873 04874 Impl() 04875 : max_weight_(NumericTraits<double>::min()), 04876 value_() 04877 {} 04878 04879 void reset() 04880 { 04881 max_weight_ = NumericTraits<double>::min(); 04882 value_ = element_type(); 04883 } 04884 04885 template <class Shape> 04886 void reshape(Shape const & s) 04887 { 04888 detail::reshapeImpl(value_, s); 04889 } 04890 04891 void operator+=(Impl const & o) 04892 { 04893 using namespace multi_math; 04894 if(o.max_weight_ > max_weight_) 04895 { 04896 max_weight_ = o.max_weight_; 04897 value_ = o.value_; 04898 } 04899 } 04900 04901 void update(U const & t) 04902 { 04903 vigra_precondition(false, "ArgMaxWeight::update() needs weights."); 04904 } 04905 04906 void update(U const & t, double weight) 04907 { 04908 if(weight > max_weight_) 04909 { 04910 max_weight_ = weight; 04911 value_ = t; 04912 } 04913 } 04914 04915 result_type operator()() const 04916 { 04917 return value_; 04918 } 04919 }; 04920 }; 04921 04922 04923 template <class BASE, int BinCount> 04924 class HistogramBase 04925 : public BASE 04926 { 04927 public: 04928 04929 typedef double element_type; 04930 typedef TinyVector<double, BinCount> value_type; 04931 typedef value_type const & result_type; 04932 04933 value_type value_; 04934 double left_outliers, right_outliers; 04935 04936 HistogramBase() 04937 : value_(), 04938 left_outliers(), 04939 right_outliers() 04940 {} 04941 04942 void reset() 04943 { 04944 value_ = element_type(); 04945 left_outliers = 0.0; 04946 right_outliers = 0.0; 04947 } 04948 04949 void operator+=(HistogramBase const & o) 04950 { 04951 value_ += o.value_; 04952 left_outliers += o.left_outliers; 04953 right_outliers += o.right_outliers; 04954 } 04955 04956 result_type operator()() const 04957 { 04958 return value_; 04959 } 04960 }; 04961 04962 template <class BASE> 04963 class HistogramBase<BASE, 0> 04964 : public BASE 04965 { 04966 public: 04967 04968 typedef double element_type; 04969 typedef MultiArray<1, double> value_type; 04970 typedef value_type const & result_type; 04971 04972 value_type value_; 04973 double left_outliers, right_outliers; 04974 04975 HistogramBase() 04976 : value_(), 04977 left_outliers(), 04978 right_outliers() 04979 {} 04980 04981 void reset() 04982 { 04983 value_ = element_type(); 04984 left_outliers = 0.0; 04985 right_outliers = 0.0; 04986 } 04987 04988 void operator+=(HistogramBase const & o) 04989 { 04990 value_ += o.value_; 04991 left_outliers += o.left_outliers; 04992 right_outliers += o.right_outliers; 04993 } 04994 04995 void setBinCount(int binCount) 04996 { 04997 vigra_precondition(binCount > 0, 04998 "HistogramBase:.setBinCount(): binCount > 0 required."); 04999 value_type(Shape1(binCount)).swap(value_); 05000 } 05001 05002 result_type operator()() const 05003 { 05004 return value_; 05005 } 05006 }; 05007 05008 template <class BASE, int BinCount, class U=typename BASE::input_type> 05009 class RangeHistogramBase 05010 : public HistogramBase<BASE, BinCount> 05011 { 05012 public: 05013 double scale_, offset_, inverse_scale_; 05014 05015 RangeHistogramBase() 05016 : scale_(), 05017 offset_(), 05018 inverse_scale_() 05019 {} 05020 05021 void reset() 05022 { 05023 scale_ = 0.0; 05024 offset_ = 0.0; 05025 inverse_scale_ = 0.0; 05026 HistogramBase<BASE, BinCount>::reset(); 05027 } 05028 05029 void operator+=(RangeHistogramBase const & o) 05030 { 05031 vigra_precondition(scale_ == 0.0 || o.scale_ == 0.0 || (scale_ == o.scale_ && offset_ == o.offset_), 05032 "RangeHistogramBase::operator+=(): cannot merge histograms with different data mapping."); 05033 05034 HistogramBase<BASE, BinCount>::operator+=(o); 05035 if(scale_ == 0.0) 05036 { 05037 scale_ = o.scale_; 05038 offset_ = o.offset_; 05039 inverse_scale_ = o.inverse_scale_; 05040 } 05041 } 05042 05043 void update(U const & t) 05044 { 05045 update(t, 1.0); 05046 } 05047 05048 void update(U const & t, double weight) 05049 { 05050 double m = mapItem(t); 05051 int index = (m == (double)this->value_.size()) 05052 ? (int)m - 1 05053 : (int)m; 05054 if(index < 0) 05055 this->left_outliers += weight; 05056 else if(index >= (int)this->value_.size()) 05057 this->right_outliers += weight; 05058 else 05059 this->value_[index] += weight; 05060 } 05061 05062 void setMinMax(double mi, double ma) 05063 { 05064 vigra_precondition(this->value_.size() > 0, 05065 "RangeHistogramBase::setMinMax(...): setBinCount(...) has not been called."); 05066 vigra_precondition(mi < ma, 05067 "RangeHistogramBase::setMinMax(...): min < max required."); 05068 offset_ = mi; 05069 scale_ = (double)this->value_.size() / (ma - mi); 05070 inverse_scale_ = 1.0 / scale_; 05071 } 05072 05073 double mapItem(double t) const 05074 { 05075 return scale_ * (t - offset_); 05076 } 05077 05078 double mapItemInverse(double t) const 05079 { 05080 return inverse_scale_ * t + offset_; 05081 } 05082 05083 template <class ArrayLike> 05084 void computeStandardQuantiles(double minimum, double maximum, double count, 05085 ArrayLike const & desiredQuantiles, ArrayLike & res) const 05086 { 05087 ArrayVector<double> keypoints, cumhist; 05088 double mappedMinimum = mapItem(minimum); 05089 double mappedMaximum = mapItem(maximum); 05090 05091 keypoints.push_back(mappedMinimum); 05092 cumhist.push_back(0.0); 05093 05094 if(this->left_outliers > 0.0) 05095 { 05096 keypoints.push_back(0.0); 05097 cumhist.push_back(this->left_outliers); 05098 } 05099 05100 int size = (int)this->value_.size(); 05101 double cumulative = this->left_outliers; 05102 for(int k=0; k<size; ++k) 05103 { 05104 if(this->value_[k] > 0.0) 05105 { 05106 if(keypoints.back() <= k) 05107 { 05108 keypoints.push_back(k); 05109 cumhist.push_back(cumulative); 05110 } 05111 cumulative += this->value_[k]; 05112 keypoints.push_back(k+1); 05113 cumhist.push_back(cumulative); 05114 } 05115 } 05116 05117 if(this->right_outliers > 0.0) 05118 { 05119 if(keypoints.back() != size) 05120 { 05121 keypoints.push_back(size); 05122 cumhist.push_back(cumulative); 05123 } 05124 keypoints.push_back(mappedMaximum); 05125 cumhist.push_back(count); 05126 } 05127 else 05128 { 05129 keypoints.back() = mappedMaximum; 05130 cumhist.back() = count; 05131 } 05132 05133 int quantile = 0, end = (int)desiredQuantiles.size(); 05134 05135 if(desiredQuantiles[0] == 0.0) 05136 { 05137 res[0] = minimum; 05138 ++quantile; 05139 } 05140 if(desiredQuantiles[end-1] == 1.0) 05141 { 05142 res[end-1] = maximum; 05143 --end; 05144 } 05145 05146 int point = 0; 05147 double qcount = count * desiredQuantiles[quantile]; 05148 while(quantile < end) 05149 { 05150 if(cumhist[point] < qcount && cumhist[point+1] >= qcount) 05151 { 05152 double t = (qcount - cumhist[point]) / (cumhist[point+1] - cumhist[point]) * (keypoints[point+1] - keypoints[point]); 05153 res[quantile] = mapItemInverse(t + keypoints[point]); 05154 ++quantile; 05155 qcount = count * desiredQuantiles[quantile]; 05156 } 05157 else 05158 { 05159 ++point; 05160 } 05161 } 05162 } 05163 }; 05164 05165 /** \brief Histogram where data values are equal to bin indices. 05166 05167 - If BinCount != 0, the return type of the accumulator is TinyVector<double, BinCount> . 05168 - If BinCount == 0, the return type of the accumulator is MultiArray<1, double> . BinCount can be set by calling getAccumulator<IntegerHistogram<0> >(acc_chain).setBinCount(bincount). 05169 - Outliers can be accessed via getAccumulator<IntegerHistogram<Bincount>>(a).left_outliers and getAccumulator<...>(acc_chain).right_outliers. 05170 - Note that histogram options (for all histograms in the accumulator chain) can also be set by passing an instance of HistogramOptions to the accumulator chain via acc_chain.setHistogramOptions(). 05171 Works in pass 1, %operator+=() supported (merging supported). 05172 */ 05173 template <int BinCount> 05174 class IntegerHistogram 05175 { 05176 public: 05177 05178 typedef Select<> Dependencies; 05179 05180 static std::string const & name() 05181 { 05182 static const std::string n = std::string("IntegerHistogram<") + asString(BinCount) + ">"; 05183 return n; 05184 } 05185 05186 template <class U, class BASE> 05187 struct Impl 05188 : public HistogramBase<BASE, BinCount> 05189 { 05190 void update(int index) 05191 { 05192 if(index < 0) 05193 ++this->left_outliers; 05194 else if(index >= (int)this->value_.size()) 05195 ++this->right_outliers; 05196 else 05197 ++this->value_[index]; 05198 } 05199 05200 void update(int index, double weight) 05201 { 05202 // cannot compute quantile from weighted integer histograms, 05203 // so force people to use UserRangeHistogram or AutoRangeHistogram 05204 vigra_precondition(false, "IntegerHistogram::update(): weighted histograms not supported, use another histogram type."); 05205 } 05206 05207 template <class ArrayLike> 05208 void computeStandardQuantiles(double minimum, double maximum, double count, 05209 ArrayLike const & desiredQuantiles, ArrayLike & res) const 05210 { 05211 int quantile = 0, end = (int)desiredQuantiles.size(); 05212 05213 if(desiredQuantiles[0] == 0.0) 05214 { 05215 res[0] = minimum; 05216 ++quantile; 05217 } 05218 if(desiredQuantiles[end-1] == 1.0) 05219 { 05220 res[end-1] = maximum; 05221 --end; 05222 } 05223 05224 count -= 1.0; 05225 int currentBin = 0, size = (int)this->value_.size(); 05226 double cumulative1 = this->left_outliers, 05227 cumulative2 = this->value_[currentBin] + cumulative1; 05228 05229 // add a to the quantiles to account for the fact that counting 05230 // corresponds to 1-based indexing (one element == index 1) 05231 double qcount = desiredQuantiles[quantile]*count + 1.0; 05232 05233 while(quantile < end) 05234 { 05235 if(cumulative2 == qcount) 05236 { 05237 res[quantile] = currentBin; 05238 ++quantile; 05239 qcount = desiredQuantiles[quantile]*count + 1.0; 05240 } 05241 else if(cumulative2 > qcount) 05242 { 05243 if(cumulative1 > qcount) // in left_outlier bin 05244 { 05245 res[quantile] = minimum; 05246 } 05247 if(cumulative1 + 1.0 > qcount) // between bins 05248 { 05249 res[quantile] = currentBin - 1 + qcount - std::floor(qcount); 05250 } 05251 else // standard case 05252 { 05253 res[quantile] = currentBin; 05254 } 05255 ++quantile; 05256 qcount = desiredQuantiles[quantile]*count + 1.0; 05257 } 05258 else if(currentBin == size-1) // in right outlier bin 05259 { 05260 res[quantile] = maximum; 05261 ++quantile; 05262 qcount = desiredQuantiles[quantile]*count + 1.0; 05263 } 05264 else 05265 { 05266 ++currentBin; 05267 cumulative1 = cumulative2; 05268 cumulative2 += this->value_[currentBin]; 05269 } 05270 } 05271 } 05272 }; 05273 }; 05274 05275 /** \brief Histogram where user provides bounds for linear range mapping from values to indices. 05276 05277 - If BinCount != 0, the return type of the accumulator is TinyVector<double, BinCount> . 05278 - If BinCount == 0, the return type of the accumulator is MultiArray<1, double> . BinCount can be set by calling getAccumulator<UserRangeHistogram<0> >(acc_chain).setBinCount(bincount). 05279 - Bounds for the mapping (min/max) must be set before seeing data by calling getAccumulator<UserRangeHistogram<BinCount> >.setMinMax(min, max). 05280 - Options can also be passed to the accumulator chain via an instance of HistogramOptions . 05281 - Works in pass 1, %operator+=() is supported (merging) if both histograms have the same data mapping. 05282 - Outliers can be accessed via getAccumulator<...>(a).left_outliers and getAccumulator<...>(a).right_outliers. 05283 - Note that histogram options (for all histograms in the accumulator chain) can also be set by passing an instance of HistogramOptions to the accumulator chain via acc_chain.setHistogramOptions(). 05284 */ 05285 template <int BinCount> 05286 class UserRangeHistogram 05287 { 05288 public: 05289 05290 typedef Select<> Dependencies; 05291 05292 static std::string const & name() 05293 { 05294 static const std::string n = std::string("UserRangeHistogram<") + asString(BinCount) + ">"; 05295 return n; 05296 } 05297 05298 template <class U, class BASE> 05299 struct Impl 05300 : public RangeHistogramBase<BASE, BinCount, U> 05301 { 05302 void update(U const & t) 05303 { 05304 update(t, 1.0); 05305 } 05306 05307 void update(U const & t, double weight) 05308 { 05309 vigra_precondition(this->scale_ != 0.0, 05310 "UserRangeHistogram::update(): setMinMax(...) has not been called."); 05311 05312 RangeHistogramBase<BASE, BinCount, U>::update(t, weight); 05313 } 05314 }; 05315 }; 05316 05317 /** \brief Histogram where range mapping bounds are defined by minimum and maximum of data. 05318 05319 - If BinCount != 0, the return type of the accumulator is TinyVector<double, BinCount> . 05320 - If BinCount == 0, the return type of the accumulator is MultiArray<1, double> . BinCount can be set by calling getAccumulator<AutoRangeHistogram>(acc_chain).setBinCount(bincount). 05321 - Becomes a UserRangeHistogram if min/max is set. 05322 - Works in pass 2, %operator+=() is supported (merging) if both histograms have the same data mapping. 05323 - Outliers can be accessed via getAccumulator<...>(acc_chain).left_outliers and getAccumulator<...>(acc_chain).right_outliers . 05324 - Note that histogram options (for all histograms in the accumulator chain) can also be set by passing an instance of HistogramOptions to the accumulator chain via acc_chain.setHistogramOptions(). 05325 */ 05326 template <int BinCount> 05327 class AutoRangeHistogram 05328 { 05329 public: 05330 05331 typedef Select<Minimum, Maximum> Dependencies; 05332 05333 static std::string const & name() 05334 { 05335 static const std::string n = std::string("AutoRangeHistogram<") + asString(BinCount) + ">"; 05336 return n; 05337 } 05338 05339 template <class U, class BASE> 05340 struct Impl 05341 : public RangeHistogramBase<BASE, BinCount, U> 05342 { 05343 static const unsigned int workInPass = LookupDependency<Minimum, BASE>::type::workInPass + 1; 05344 05345 void update(U const & t) 05346 { 05347 update(t, 1.0); 05348 } 05349 05350 void update(U const & t, double weight) 05351 { 05352 if(this->scale_ == 0.0) 05353 this->setMinMax(getDependency<Minimum>(*this), getDependency<Maximum>(*this)); 05354 05355 RangeHistogramBase<BASE, BinCount, U>::update(t, weight); 05356 } 05357 }; 05358 }; 05359 05360 /** \brief Like AutoRangeHistogram, but use global min/max rather than region min/max. 05361 05362 - If BinCount != 0, the return type of the accumulator is TinyVector<double, BinCount> . 05363 - If BinCount == 0, the return type of the accumulator is MultiArray<1, double> . BinCount can be set by calling getAccumulator<GlobalRangeHistogram<0>>(acc_chain).setBinCount(bincount). 05364 - Becomes a UserRangeHistogram if min/max is set. 05365 - Works in pass 2, %operator+=() is supported (merging) if both histograms have the same data mapping. 05366 - Outliers can be accessed via getAccumulator<GlobalRangeHistogram<Bincount>>(acc_chain).left_outliers and getAccumulator<...>(acc_chain).right_outliers . 05367 - Histogram options (for all histograms in the accumulator chain) can also be set by passing an instance of HistogramOptions to the accumulator chain via acc_chain.setHistogramOptions(). 05368 */ 05369 template <int BinCount> 05370 class GlobalRangeHistogram 05371 { 05372 public: 05373 05374 typedef Select<Global<Minimum>, Global<Maximum>, Minimum, Maximum> Dependencies; 05375 05376 static std::string const & name() 05377 { 05378 static const std::string n = std::string("GlobalRangeHistogram<") + asString(BinCount) + ">"; 05379 return n; 05380 } 05381 05382 template <class U, class BASE> 05383 struct Impl 05384 : public RangeHistogramBase<BASE, BinCount, U> 05385 { 05386 static const unsigned int workInPass = LookupDependency<Minimum, BASE>::type::workInPass + 1; 05387 05388 bool useLocalMinimax_; 05389 05390 Impl() 05391 : useLocalMinimax_(false) 05392 {} 05393 05394 void setRegionAutoInit(bool locally) 05395 { 05396 this->scale_ = 0.0; 05397 useLocalMinimax_ = locally; 05398 } 05399 05400 void update(U const & t) 05401 { 05402 update(t, 1.0); 05403 } 05404 05405 void update(U const & t, double weight) 05406 { 05407 if(this->scale_ == 0.0) 05408 { 05409 if(useLocalMinimax_) 05410 this->setMinMax(getDependency<Minimum>(*this), getDependency<Maximum>(*this)); 05411 else 05412 this->setMinMax(getDependency<Global<Minimum> >(*this), getDependency<Global<Maximum> >(*this)); 05413 } 05414 05415 RangeHistogramBase<BASE, BinCount, U>::update(t, weight); 05416 } 05417 }; 05418 }; 05419 05420 /** \brief Compute (0%, 10%, 25%, 50%, 75%, 90%, 100%) quantiles from given histogram. 05421 05422 Return type is TinyVector<double, 7> . 05423 */ 05424 template <class HistogramAccumulator> 05425 class StandardQuantiles 05426 { 05427 public: 05428 05429 typedef typename StandardizeTag<HistogramAccumulator>::type HistogramTag; 05430 typedef Select<HistogramTag, Minimum, Maximum, Count> Dependencies; 05431 05432 static std::string const & name() 05433 { 05434 static const std::string n = std::string("StandardQuantiles<") + HistogramTag::name() + " >"; 05435 return n; 05436 } 05437 05438 template <class U, class BASE> 05439 struct Impl 05440 : public CachedResultBase<BASE, TinyVector<double, 7>, U> 05441 { 05442 typedef typename CachedResultBase<BASE, TinyVector<double, 7>, U>::result_type result_type; 05443 typedef typename CachedResultBase<BASE, TinyVector<double, 7>, U>::value_type value_type; 05444 05445 static const unsigned int workInPass = LookupDependency<HistogramTag, BASE>::type::workInPass; 05446 05447 result_type operator()() const 05448 { 05449 if(this->isDirty()) 05450 { 05451 static const double desiredQuantiles[] = {0.0, 0.1, 0.25, 0.5, 0.75, 0.9, 1.0 }; 05452 getAccumulator<HistogramTag>(*this).computeStandardQuantiles(getDependency<Minimum>(*this), getDependency<Maximum>(*this), 05453 getDependency<Count>(*this), value_type(desiredQuantiles), 05454 this->value_); 05455 this->setClean(); 05456 } 05457 return this->value_; 05458 } 05459 }; 05460 }; 05461 05462 }} // namespace vigra::acc 05463 05464 #endif // VIGRA_ACCUMULATOR_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|