[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

accumulator.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2011-2012 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_ACCUMULATOR_HXX
37 #define VIGRA_ACCUMULATOR_HXX
38 
39 #ifdef _MSC_VER
40 #pragma warning (disable: 4503)
41 #endif
42 
43 #include "accumulator-grammar.hxx"
44 #include "config.hxx"
45 #include "metaprogramming.hxx"
46 #include "bit_array.hxx"
47 #include "static_assert.hxx"
48 #include "mathutil.hxx"
49 #include "utilities.hxx"
50 #include "multi_iterator_coupled.hxx"
51 #include "matrix.hxx"
52 #include "multi_math.hxx"
53 #include "eigensystem.hxx"
54 #include "histogram.hxx"
55 #include <algorithm>
56 #include <iostream>
57 
58 namespace vigra {
59 
60 /** \defgroup FeatureAccumulators Feature Accumulators
61 
62 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.
63 
64 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".
65 
66 <b>\#include</b> <vigra/accumulator.hxx>
67 
68  <b>Basic statistics:</b>
69  - PowerSum<N> (@f$ \sum_i x_i^N @f$)
70  - AbsPowerSum<N> (@f$ \sum_i |x_i|^N @f$)
71  - Skewness, UnbiasedSkewness
72  - Kurtosis, UnbiasedKurtosis
73  - Minimum, Maximum
74  - FlatScatterMatrix (flattened upper-triangular part of scatter matrix)
75  - 4 histogram classes (see \ref histogram "below")
76  - StandardQuantiles (0%, 10%, 25%, 50%, 75%, 90%, 100%)
77  - ArgMinWeight, ArgMaxWeight (store data or coordinate where weight assumes its minimal or maximal value)
78  - CoordinateSystem (identity matrix of appropriate size)
79 
80  <b>Modifiers:</b> (S is the statistc to be modified)
81  - Normalization
82  <table border="0">
83  <tr><td> DivideByCount<S> </td><td> S/Count </td></tr>
84  <tr><td> RootDivideByCount<S> </td><td> sqrt( S/Count ) </td></tr>
85  <tr><td> DivideUnbiased<S> </td><td> S/(Count-1) </td></tr>
86  <tr><td> RootDivideUnbiased<S> &nbsp; &nbsp; </td><td> sqrt( S/(Count-1) ) </td></tr>
87  </table>
88  - Data preparation:
89  <table border="0">
90  <tr><td> Central<S> </td><td> substract mean before computing S </td></tr>
91  <tr><td> Principal<S> </td><td> project onto PCA eigenvectors </td></tr>
92  <tr><td> Whitened<S> &nbsp; &nbsp; </td><td> scale to unit variance after PCA </td></tr>
93  <tr><td> Coord<S> </td><td> compute S from pixel coordinates rather than from pixel values </td></tr>
94  <tr><td> Weighted<S> </td><td> compute weighted version of S </td></tr>
95  <tr><td> Global<S> </td><td> compute S globally rather than per region (per region is default if labels are given) </td></tr>
96  </table>
97 
98  Aliases for a couple of important features are implemented (mainly as <tt>typedef FullName Alias</tt>). The alias names are equivalent to full names.
99  Here are some examples for supported alias names (these examples also show how to compose statistics from the fundamental statistics and modifiers):
100 
101  <table border="0">
102  <tr><th> Alias </th><th> Full Name </th></tr>
103  <tr><td> Count </td><td> PowerSum<0> </td></tr>
104  <tr><td> Sum </td><td> PowerSum<1> </td></tr>
105  <tr><td> SumOfSquares </td><td> PowerSum<2> </td></tr>
106  <tr><td> Mean </td><td> DivideByCount<PowerSum<1>> </td></tr>
107  <tr><td> RootMeanSquares &nbsp; </td><td> RootDivideByCount<PowerSum<2>> </td></tr>
108  <tr><td> Moment<N> </td><td> DivideByCount<PowerSum<N>> </td></tr>
109  <tr><td> Variance </td><td> DivideByCount<Central<PowerSum<2>>> </td></tr>
110  <tr><td> StdDev </td><td> RootDivideByCount<Central<PowerSum<2>>> </td></tr>
111  <tr><td> Covariance </td><td> DivideByCount<FlatScatterMatrix> </td></tr>
112  <tr><td> RegionCenter </td><td> Coord<Mean> </td></tr>
113  <tr><td> CenterOfMass </td><td> Weighted<Coord<Mean>> </td></tr>
114  </table>
115 
116  There are a few <b>rules for composing statistics</b>:
117  - modifiers can be specified in any order, but are internally transformed to standard order: Global<Weighted<Coord<normalization<data preparation<basic statistic
118  - only one normalization modifier and one data preparation modifier (Central or Principal or Whitened) is permitted
119  - Count ignores all modifiers except Global and Weighted
120  - Sum ignores Central and Principal, because sum would be zero
121  - ArgMinWeight and ArgMaxWeight are automatically Weighted
122 
123 
124  Here is an example how to use \ref acc::AccumulatorChain to compute statistics. (To use Weighted<> or Coord<> modifiers, see below):
125 
126  \code
127  #include <vigra/multi_array.hxx>
128  #include <vigra/impex.hxx>
129  #include <vigra/accumulator.hxx>
130  using namespace vigra::acc;
131  typedef double DataType;
132  int size = 1000;
133  vigra::MultiArray<2, DataType> data(vigra::Shape2(size, size));
134 
135  AccumulatorChain<DataType,
136  Select<Variance, Mean, StdDev, Minimum, Maximum, RootMeanSquares, Skewness, Covariance> >
137  a;
138 
139  std::cout << "passes required: " << a.passesRequired() << std::endl;
140  extractFeatures(data.begin(), data.end(), a);
141 
142  std::cout << "Mean: " << get<Mean>(a) << std::endl;
143  std::cout << "Variance: " << get<Variance>(a) << std::endl;
144  \endcode
145 
146  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 .
147 
148  Rules and notes:
149  - order of statistics in Select<> is arbitrary
150  - up to 20 statistics in Select<>, but Select<> can be nested
151  - dependencies are automatically inserted
152  - duplicates are automatically removed
153  - extractFeatures() does as many passes through the data as necessary
154  - each accumulator only sees data in the appropriate pass (its "working pass")
155 
156  The Accumulators can also be used with vector-valued data (vigra::RGBValue, vigra::TinyVector, vigra::MultiArray or vigra::MultiArrayView):
157 
158  \code
159  typedef vigra::RGBValue<double> DataType;
160  AccumulatorChain<DataType, Select<...> > a;
161  ...
162  \endcode
163 
164  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.
165 
166 These <b>index specifiers</b> are: (INDEX is of type int)
167 
168  - DataArg<INDEX>: CoupledHandle holds data at index 'INDEX' (default INDEX=1)
169  - LabelArg<INDEX>: CoupledHandle holds labels at index 'INDEX' (default INDEX=2)
170  - WeightArg<INDEX>: CoupledHandle holds weights at index 'INDEX' (default INDEX=outermost index)
171 
172 Pixel coordinates are always at index 0.
173 
174  \code
175  using namespace vigra::acc;
176  vigra::MultiArray<3, double> data(...), weights(...);
177  typedef vigra::CoupledIteratorType<3, double, double>::type Iterator; //type of the CoupledScanOrderIterator
178  typedef Iterator::value_type Handle; //type of the corresponding CoupledHandle
179 
180  AccumulatorChain<Handle,
181  Select<DataArg<1>, WeightArg<2>, //where to look in the Handle (coordinates are always arg 0)
182  Mean, Variance, //statistics over values
183  Coord<Mean>, Coord<Variance>, //statistics over coordinates,
184  Weighted<Mean>, Weighted<Variance>, //weighted values,
185  Weighted<Coord<Mean> > > > //weighted coordinates.
186  a;
187 
188  Iterator start = createCoupledIterator(data, weights); //coord->index 0, data->index 1, weights->index 2
189  Iterator end = start.getEndIterator();
190 
191  extractFeatures(start,end,a);
192  \endcode
193 
194  To compute <b>region statistics</b>, use \ref acc::AccumulatorChainArray :
195 
196  \code
197  using namespace vigra::acc;
198  vigra::MultiArray<3, double> data(...);
199  vigra::MultiArray<3, int> labels(...);
200  typedef vigra::CoupledIteratorType<3, double, int>::type Iterator;
201  typedef Iterator::value_type Handle;
202 
203  AccumulatorChainArray<Handle,
204  Select<DataArg<1>, LabelArg<2>, //where to look in the Handle (coordinates are always arg 0)
205  Mean, Variance, //per-region statistics over values
206  Coord<Mean>, Coord<Variance>, //per-region statistics over coordinates
207  Global<Mean>, Global<Variance> > > //global statistics
208  a;
209 
210  Iterator start = createCoupledIterator(data, labels);
211  Iterator end = start.getEndIterator();
212 
213  a.ignoreLabel(0); //statistics will not be computed for region 0 (e.g. background)
214 
215  extractFeatures(start,end,a);
216 
217  int regionlabel = ...;
218  std::cout << get<Mean>(a, regionlabel) << std::endl; //get Mean of region with label 'regionlabel'
219  \endcode
220 
221 
222  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:
223 
224  \code
225  using namespace vigra::acc;
226  vigra::MultiArray<2, double> data(...);
227  DynamicAccumulatorChain<double,
228  Select<Mean, Minimum, Maximum, Variance, StdDev> > a; // at compile-time
229  activate<Mean>(a); //at run-time
230  a.activate("Minimum"); //same as activate<Minimum>(a) (alias names are not recognized)
231 
232  extractFeatures(data.begin(), data.end(), a);
233  std::cout << "Mean: " << get<Mean>(a) << std::endl; //ok
234  //std::cout << "Maximum: " << get<Maximum>(a) << std::endl; // run-time error because Maximum not activated
235  \endcode
236 
237  Likewise, for run-time activation of region statistics, use \ref acc::DynamicAccumulatorChainArray.
238 
239  <b>Accumulator merging</b> (e.g. for parallelization or hierarchical segmentation) is possible for many accumulators:
240 
241  \code
242  using namespace vigra::acc;
243  vigra::MultiArray<2, double> data(...);
244  AccumulatorChain<double, Select<Mean, Variance, Skewness> > a, a1, a2;
245 
246  extractFeatures(data.begin(), data.end(), a); //process entire data set at once
247  extractFeatures(data.begin(), data.begin()+data.size()/2, a1); //process first half
248  extractFeatures(data.begin()+data.size()/2, data.end(), a2); //process second half
249  a1 += a2; // merge: a1 now equals a0 (with numerical tolerances)
250  \endcode
251 
252  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:
253 
254  \code
255  using namespace vigra::acc;
256  vigra::MultiArray<2, double> data(...);
257  AccumulatorChain<double, Select<Mean, Variance> > a;
258 
259  extractFeatures(data.begin(), data.begin()+data.size()/2, a); // this works because
260  extractFeatures(data.begin()+data.size()/2, data.end(), a); // all statistics only work in pass 1
261 
262  \endcode
263 
264 
265  \anchor histogram
266  Four kinds of <b>histograms</b> are currently implemented:
267 
268  <table border="0">
269  <tr><td> IntegerHistogram </td><td> Data values are equal to bin indices </td></tr>
270  <tr><td> UserRangeHistogram </td><td> User provides lower and upper bounds for linear range mapping from values to indices. </td></tr>
271  <tr><td> AutoRangeHistogram </td><td> Range mapping bounds are defiend by minimum and maximum of the data (2 passes needed!) </td></tr>
272  <tr><td> GlobalRangeHistogram &nbsp; </td><td> Likewise, but use global min/max rather than region min/max as AutoRangeHistogram will </td></tr>
273  </table>
274 
275 
276 
277  - 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).
278  - If UserRangeHistogram is used, the bounds for the linear range mapping from values to indices must be set before seeing data (see below).
279  - 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.
280  - Merging is supported if the range mapping of the histograms is the same.
281  - Histogram accumulators have two members for outliers (left_outliers, right_outliers).
282 
283  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> .
284 
285  \anchor acc_hist_options Usage:
286  \code
287  using namespace vigra::acc;
288  typedef double DataType;
289  vigra::MultiArray<2, DataType> data(...);
290 
291  typedef UserRangeHistogram<40> SomeHistogram; //binCount set at compile time
292  typedef UserRangeHistogram<0> SomeHistogram2; // binCount must be set at run-time
293  typedef AutoRangeHistogram<0> SomeHistogram3;
294  typedef StandardQuantiles<SomeHistogram3> Quantiles3;
295 
296  AccumulatorChain<DataType, Select<SomeHistogram, SomeHistogram2, SomeHistogram3, Quantiles3> > a;
297 
298  //set options for all histograms in the accumulator chain:
299  vigra::HistogramOptions histogram_opt;
300  histogram_opt = histogram_opt.setBinCount(50);
301  //histogram_opt = histogram_opt.setMinMax(0.1, 0.9); // this would set min/max for all three histograms, but range bounds
302  // shall be set automatically by min/max of data for SomeHistogram3
303  a.setHistogramOptions(histogram_opt);
304 
305  // set options for a specific histogram in the accumulator chain:
306  getAccumulator<SomeHistogram>(a).setMinMax(0.1, 0.9); // number of bins must be set before setting min/max
307  getAccumulator<SomeHistogram2>(a).setMinMax(0.0, 1.0);
308 
309  extractFeatures(data.begin(), data.end(), a);
310 
311  vigra::TinyVector<double, 40> hist = get<SomeHistogram>(a);
312  vigra::MultiArray<1, double> hist2 = get<SomeHistogram2>(a);
313  vigra::TinyVector<double, 7> quant = get<Quantiles3>(a);
314  double right_outliers = getAccumulator<SomeHistogram>(a).right_outliers;
315  \endcode
316 
317 
318 
319 */
320 
321 
322 /** This namespace contains the accumulator classes, fundamental statistics and modifiers. See \ref FeatureAccumulators for examples of usage.
323 */
324 namespace acc {
325 
326 /****************************************************************************/
327 /* */
328 /* infrastructure */
329 /* */
330 /****************************************************************************/
331 
332  /// \brief Wrapper for MakeTypeList that additionally performs tag standardization.
333 
334 template <class T01=void, class T02=void, class T03=void, class T04=void, class T05=void,
335  class T06=void, class T07=void, class T08=void, class T09=void, class T10=void,
336  class T11=void, class T12=void, class T13=void, class T14=void, class T15=void,
337  class T16=void, class T17=void, class T18=void, class T19=void, class T20=void>
338 struct Select
339 : public MakeTypeList<
340  typename StandardizeTag<T01>::type, typename StandardizeTag<T02>::type, typename StandardizeTag<T03>::type,
341  typename StandardizeTag<T04>::type, typename StandardizeTag<T05>::type, typename StandardizeTag<T06>::type,
342  typename StandardizeTag<T07>::type, typename StandardizeTag<T08>::type, typename StandardizeTag<T09>::type,
343  typename StandardizeTag<T10>::type, typename StandardizeTag<T11>::type, typename StandardizeTag<T12>::type,
344  typename StandardizeTag<T13>::type, typename StandardizeTag<T14>::type, typename StandardizeTag<T15>::type,
345  typename StandardizeTag<T16>::type, typename StandardizeTag<T17>::type, typename StandardizeTag<T18>::type,
346  typename StandardizeTag<T19>::type, typename StandardizeTag<T20>::type
347  >
348 {};
349 
350  // enable nesting of Select<> expressions
351 template <class T01, class T02, class T03, class T04, class T05,
352  class T06, class T07, class T08, class T09, class T10,
353  class T11, class T12, class T13, class T14, class T15,
354  class T16, class T17, class T18, class T19, class T20>
355 struct StandardizeTag<Select<T01, T02, T03, T04, T05,
356  T06, T07, T08, T09, T10,
357  T11, T12, T13, T14, T15,
358  T16, T17, T18, T19, T20>,
359  Select<T01, T02, T03, T04, T05,
360  T06, T07, T08, T09, T10,
361  T11, T12, T13, T14, T15,
362  T16, T17, T18, T19, T20> >
363 {
364  typedef typename Select<T01, T02, T03, T04, T05,
365  T06, T07, T08, T09, T10,
366  T11, T12, T13, T14, T15,
367  T16, T17, T18, T19, T20>::type type;
368 };
369 
370 struct AccumulatorBegin
371 {
372  typedef Select<> Dependencies;
373 
374  static std::string const & name()
375  {
376  static const std::string n("AccumulatorBegin (internal)");
377  return n;
378  }
379 
380  template <class T, class BASE>
381  struct Impl
382  : public BASE
383  {};
384 };
385 
386 
387 struct AccumulatorEnd;
388 struct DataArgTag;
389 struct WeightArgTag;
390 struct LabelArgTag;
391 struct CoordArgTag;
392 struct LabelDispatchTag;
393 
394 struct Error__Global_statistics_are_only_defined_for_AccumulatorChainArray;
395 
396 /** \brief Specifies index of labels in CoupledHandle.
397 
398  LabelArg<INDEX> tells the acc::AccumulatorChainArray which index of the Handle contains the labels. (Note that coordinates are always index 0)
399  */
400 template <int INDEX>
401 class LabelArg
402 {
403  public:
404  typedef Select<> Dependencies;
405 
406  static std::string const & name()
407  {
408  static const std::string n = std::string("LabelArg<") + asString(INDEX) + "> (internal)";
409  return n;
410  }
411 
412  template <class T, class BASE>
413  struct Impl
414  : public BASE
415  {
416  typedef LabelArgTag Tag;
417  typedef void value_type;
418  typedef void result_type;
419 
420  static const int value = INDEX;
421  static const unsigned int workInPass = 0;
422  };
423 };
424 
425 template <int INDEX>
426 class CoordArg
427 {
428  public:
429  typedef Select<> Dependencies;
430 
431  static std::string const & name()
432  {
433  static const std::string n = std::string("CoordArg<") + asString(INDEX) + "> (internal)";
434  return n;
435  }
436 
437  template <class T, class BASE>
438  struct Impl
439  : public BASE
440  {
441  typedef CoordArgTag Tag;
442  typedef void value_type;
443  typedef void result_type;
444 
445  static const int value = INDEX;
446  static const unsigned int workInPass = 0;
447  };
448 };
449 
450 template <class T, class TAG, class NEXT=AccumulatorEnd>
451 struct AccumulatorBase;
452 
453 template <class Tag, class A>
454 struct LookupTag;
455 
456 template <class Tag, class A, class TargetTag=typename A::Tag>
457 struct LookupDependency;
458 
459 #ifndef _MSC_VER // compiler bug? (causes 'ambiguous overload error')
460 
461 template <class TAG, class A>
462 typename LookupTag<TAG, A>::reference
463 getAccumulator(A & a);
464 
465 template <class TAG, class A>
466 typename LookupDependency<TAG, A>::result_type
467 getDependency(A const & a);
468 
469 #endif
470 
471 namespace detail {
472 
473 /****************************************************************************/
474 /* */
475 /* internal tag handling meta-functions */
476 /* */
477 /****************************************************************************/
478 
479  // we must make sure that Arg<INDEX> tags are at the end of the chain because
480  // all other tags potentially depend on them
481 template <class T>
482 struct PushArgTagToTail
483 {
484  typedef T type;
485 };
486 
487 #define VIGRA_PUSHARGTAG(TAG) \
488 template <int INDEX, class TAIL> \
489 struct PushArgTagToTail<TypeList<TAG<INDEX>, TAIL> > \
490 { \
491  typedef typename Push<TAIL, TypeList<TAG<INDEX> > >::type type; \
492 };
493 
494 VIGRA_PUSHARGTAG(DataArg)
495 VIGRA_PUSHARGTAG(WeightArg)
496 VIGRA_PUSHARGTAG(CoordArg)
497 VIGRA_PUSHARGTAG(LabelArg)
498 
499 #undef VIGRA_PUSHARGTAG
500 
501  // Insert the dependencies of the selected functors into the TypeList and sort
502  // the list such that dependencies come after the functors using them. Make sure
503  // that each functor is contained only once.
504 template <class T>
505 struct AddDependencies;
506 
507 template <class HEAD, class TAIL>
508 struct AddDependencies<TypeList<HEAD, TAIL> >
509 {
510  typedef typename AddDependencies<TAIL>::type TailWithDependencies;
511  typedef typename StandardizeDependencies<HEAD>::type HeadDependencies;
512  typedef typename AddDependencies<HeadDependencies>::type TransitiveHeadDependencies;
513  typedef TypeList<HEAD, TransitiveHeadDependencies> HeadWithDependencies;
514  typedef typename PushUnique<HeadWithDependencies, TailWithDependencies>::type UnsortedDependencies;
515  typedef typename PushArgTagToTail<UnsortedDependencies>::type type;
516 };
517 
518 template <>
519 struct AddDependencies<void>
520 {
521  typedef void type;
522 };
523 
524  // Helper class to activate dependencies at runtime (i.e. when activate<Tag>(accu) is called,
525  // activate() must also be called for Tag's dependencies).
526 template <class Dependencies>
527 struct ActivateDependencies;
528 
529 template <class HEAD, class TAIL>
530 struct ActivateDependencies<TypeList<HEAD, TAIL> >
531 {
532  template <class Chain, class ActiveFlags>
533  static void exec(ActiveFlags & flags)
534  {
535  LookupTag<HEAD, Chain>::type::activateImpl(flags);
536  ActivateDependencies<TAIL>::template exec<Chain>(flags);
537  }
538 
539  template <class Chain, class ActiveFlags, class GlobalFlags>
540  static void exec(ActiveFlags & flags, GlobalFlags & gflags)
541  {
542  LookupTag<HEAD, Chain>::type::template activateImpl<Chain>(flags, gflags);
543  ActivateDependencies<TAIL>::template exec<Chain>(flags, gflags);
544  }
545 };
546 
547 template <class HEAD, class TAIL>
548 struct ActivateDependencies<TypeList<Global<HEAD>, TAIL> >
549 {
550  template <class Chain, class ActiveFlags, class GlobalFlags>
551  static void exec(ActiveFlags & flags, GlobalFlags & gflags)
552  {
553  LookupTag<Global<HEAD>, Chain>::type::activateImpl(gflags);
554  ActivateDependencies<TAIL>::template exec<Chain>(flags, gflags);
555  }
556 };
557 
558 template <>
559 struct ActivateDependencies<void>
560 {
561  template <class Chain, class ActiveFlags>
562  static void exec(ActiveFlags &)
563  {}
564 
565  template <class Chain, class ActiveFlags, class GlobalFlags>
566  static void exec(ActiveFlags &, GlobalFlags &)
567  {}
568 };
569 
570 template <class List>
571 struct SeparateGlobalAndRegionTags;
572 
573 template <class HEAD, class TAIL>
574 struct SeparateGlobalAndRegionTags<TypeList<HEAD, TAIL> >
575 {
576  typedef SeparateGlobalAndRegionTags<TAIL> Inner;
577  typedef TypeList<HEAD, typename Inner::RegionTags> RegionTags;
578  typedef typename Inner::GlobalTags GlobalTags;
579 };
580 
581 template <class HEAD, class TAIL>
582 struct SeparateGlobalAndRegionTags<TypeList<Global<HEAD>, TAIL> >
583 {
584  typedef SeparateGlobalAndRegionTags<TAIL> Inner;
585  typedef typename Inner::RegionTags RegionTags;
586  typedef TypeList<HEAD, typename Inner::GlobalTags> GlobalTags;
587 };
588 
589 template <int INDEX, class TAIL>
590 struct SeparateGlobalAndRegionTags<TypeList<DataArg<INDEX>, TAIL> >
591 {
592  typedef SeparateGlobalAndRegionTags<TAIL> Inner;
593  typedef TypeList<DataArg<INDEX>, typename Inner::RegionTags> RegionTags;
594  typedef TypeList<DataArg<INDEX>, typename Inner::GlobalTags> GlobalTags;
595 };
596 
597 template <int INDEX, class TAIL>
598 struct SeparateGlobalAndRegionTags<TypeList<LabelArg<INDEX>, TAIL> >
599 {
600  typedef SeparateGlobalAndRegionTags<TAIL> Inner;
601  typedef typename Inner::RegionTags RegionTags;
602  typedef TypeList<LabelArg<INDEX>, typename Inner::GlobalTags> GlobalTags;
603 };
604 
605 template <int INDEX, class TAIL>
606 struct SeparateGlobalAndRegionTags<TypeList<WeightArg<INDEX>, TAIL> >
607 {
608  typedef SeparateGlobalAndRegionTags<TAIL> Inner;
609  typedef TypeList<WeightArg<INDEX>, typename Inner::RegionTags> RegionTags;
610  typedef TypeList<WeightArg<INDEX>, typename Inner::GlobalTags> GlobalTags;
611 };
612 
613 template <int INDEX, class TAIL>
614 struct SeparateGlobalAndRegionTags<TypeList<CoordArg<INDEX>, TAIL> >
615 {
616  typedef SeparateGlobalAndRegionTags<TAIL> Inner;
617  typedef TypeList<CoordArg<INDEX>, typename Inner::RegionTags> RegionTags;
618  typedef TypeList<CoordArg<INDEX>, typename Inner::GlobalTags> GlobalTags;
619 };
620 
621 template <>
622 struct SeparateGlobalAndRegionTags<void>
623 {
624  typedef void RegionTags;
625  typedef void GlobalTags;
626 };
627 
628 /****************************************************************************/
629 /* */
630 /* helper classes to handle tags at runtime via strings */
631 /* */
632 /****************************************************************************/
633 
634 template <class Accumulators>
635 struct CollectAccumulatorNames;
636 
637 template <class HEAD, class TAIL>
638 struct CollectAccumulatorNames<TypeList<HEAD, TAIL> >
639 {
640  template <class BackInsertable>
641  static void exec(BackInsertable & a, bool skipInternals=true)
642  {
643  if(!skipInternals || HEAD::name().find("internal") == std::string::npos)
644  a.push_back(HEAD::name());
645  CollectAccumulatorNames<TAIL>::exec(a, skipInternals);
646  }
647 };
648 
649 template <>
650 struct CollectAccumulatorNames<void>
651 {
652  template <class BackInsertable>
653  static void exec(BackInsertable & a, bool skipInternals=true)
654  {}
655 };
656 
657 template <class T>
658 struct ApplyVisitorToTag;
659 
660 template <class HEAD, class TAIL>
661 struct ApplyVisitorToTag<TypeList<HEAD, TAIL> >
662 {
663  template <class Accu, class Visitor>
664  static bool exec(Accu & a, std::string const & tag, Visitor const & v)
665  {
666  static const std::string name = normalizeString(HEAD::name());
667  if(name == tag)
668  {
669  v.template exec<HEAD>(a);
670  return true;
671  }
672  else
673  {
674  return ApplyVisitorToTag<TAIL>::exec(a, tag, v);
675  }
676  }
677 };
678 
679 template <>
680 struct ApplyVisitorToTag<void>
681 {
682  template <class Accu, class Visitor>
683  static bool exec(Accu & a, std::string const & tag, Visitor const & v)
684  {
685  return false;
686  }
687 };
688 
689 struct ActivateTag_Visitor
690 {
691  template <class TAG, class Accu>
692  void exec(Accu & a) const
693  {
694  a.template activate<TAG>();
695  }
696 };
697 
698 struct TagIsActive_Visitor
699 {
700  mutable bool result;
701 
702  template <class TAG, class Accu>
703  void exec(Accu & a) const
704  {
705  result = a.template isActive<TAG>();
706  }
707 };
708 
709 /****************************************************************************/
710 /* */
711 /* histogram initialization functors */
712 /* */
713 /****************************************************************************/
714 
715 template <class TAG>
716 struct SetHistogramBincount
717 {
718  template <class Accu>
719  static void exec(Accu & a, HistogramOptions const & options)
720  {}
721 };
722 
723 template <template <int> class Histogram>
724 struct SetHistogramBincount<Histogram<0> >
725 {
726  template <class Accu>
727  static void exec(Accu & a, HistogramOptions const & options)
728  {
729  a.setBinCount(options.binCount);
730  }
731 };
732 
733 template <class TAG>
734 struct ApplyHistogramOptions
735 {
736  template <class Accu>
737  static void exec(Accu & a, HistogramOptions const & options)
738  {}
739 };
740 
741 template <class TAG>
742 struct ApplyHistogramOptions<StandardQuantiles<TAG> >
743 {
744  template <class Accu>
745  static void exec(Accu & a, HistogramOptions const & options)
746  {}
747 };
748 
749 template <class TAG, template <class> class MODIFIER>
750 struct ApplyHistogramOptions<MODIFIER<TAG> >
751 : public ApplyHistogramOptions<TAG>
752 {};
753 
754 template <>
755 struct ApplyHistogramOptions<IntegerHistogram<0> >
756 {
757  template <class Accu>
758  static void exec(Accu & a, HistogramOptions const & options)
759  {
760  SetHistogramBincount<IntegerHistogram<0> >::exec(a, options);
761  }
762 };
763 
764 template <int BinCount>
765 struct ApplyHistogramOptions<UserRangeHistogram<BinCount> >
766 {
767  template <class Accu>
768  static void exec(Accu & a, HistogramOptions const & options)
769  {
770  SetHistogramBincount<UserRangeHistogram<BinCount> >::exec(a, options);
771  if(a.scale_ == 0.0 && options.validMinMax())
772  a.setMinMax(options.minimum, options.maximum);
773  }
774 };
775 
776 template <int BinCount>
777 struct ApplyHistogramOptions<AutoRangeHistogram<BinCount> >
778 {
779  template <class Accu>
780  static void exec(Accu & a, HistogramOptions const & options)
781  {
782  SetHistogramBincount<AutoRangeHistogram<BinCount> >::exec(a, options);
783  if(a.scale_ == 0.0 && options.validMinMax())
784  a.setMinMax(options.minimum, options.maximum);
785  }
786 };
787 
788 template <int BinCount>
789 struct ApplyHistogramOptions<GlobalRangeHistogram<BinCount> >
790 {
791  template <class Accu>
792  static void exec(Accu & a, HistogramOptions const & options)
793  {
794  SetHistogramBincount<GlobalRangeHistogram<BinCount> >::exec(a, options);
795  if(a.scale_ == 0.0)
796  {
797  if(options.validMinMax())
798  a.setMinMax(options.minimum, options.maximum);
799  else
800  a.setRegionAutoInit(options.local_auto_init);
801  }
802  }
803 };
804 
805 /****************************************************************************/
806 /* */
807 /* internal accumulator chain classes */
808 /* */
809 /****************************************************************************/
810 
811  // AccumulatorEndImpl has the following functionalities:
812  // * marks end of accumulator chain by the AccumulatorEnd tag
813  // * provides empty implementation of standard accumulator functions
814  // * provides active_accumulators_ flags for run-time activation of dynamic accumulators
815  // * provides is_dirty_ flags for caching accumulators
816  // * hold the GlobalAccumulatorHandle for global accumulator lookup from region accumulators
817 template <unsigned LEVEL, class GlobalAccumulatorHandle>
818 struct AccumulatorEndImpl
819 {
820  typedef typename GlobalAccumulatorHandle::type GlobalAccumulatorType;
821 
822  typedef AccumulatorEnd Tag;
823  typedef void value_type;
824  typedef bool result_type;
825  typedef BitArray<LEVEL> AccumulatorFlags;
826 
827  static const unsigned int workInPass = 0;
828  static const int index = -1;
829  static const unsigned level = LEVEL;
830 
831  AccumulatorFlags active_accumulators_;
832  mutable AccumulatorFlags is_dirty_;
833  GlobalAccumulatorHandle globalAccumulator_;
834 
835  template <class GlobalAccumulator>
836  void setGlobalAccumulator(GlobalAccumulator const * a)
837  {
838  globalAccumulator_.pointer_ = a;
839  }
840 
841  static std::string name()
842  {
843  return "AccumulatorEnd (internal)";
844  }
845 
846  bool operator()() const { return false; }
847  bool get() const { return false; }
848 
849  template <unsigned, class U>
850  void pass(U const &)
851  {}
852 
853  template <unsigned, class U>
854  void pass(U const &, double)
855  {}
856 
857  template <class U>
858  void merge(U const &)
859  {}
860 
861  template <class U>
862  void resize(U const &)
863  {}
864 
865  void activate()
866  {}
867 
868  bool isActive() const
869  {
870  return false;
871  }
872 
873  template <class Flags>
874  static void activateImpl(Flags &)
875  {}
876 
877  template <class Accu, class Flags1, class Flags2>
878  static void activateImpl(Flags1 &, Flags2 &)
879  {}
880 
881  template <class Flags>
882  static bool isActiveImpl(Flags const &)
883  {
884  return true;
885  }
886 
887  void applyHistogramOptions(HistogramOptions const &)
888  {}
889 
890  static unsigned int passesRequired()
891  {
892  return 0;
893  }
894 
895  static unsigned int passesRequired(AccumulatorFlags const &)
896  {
897  return 0;
898  }
899 
900  void reset()
901  {
902  active_accumulators_.clear();
903  is_dirty_.clear();
904  }
905 
906  template <int which>
907  void setDirtyImpl() const
908  {
909  is_dirty_.template set<which>();
910  }
911 
912  template <int which>
913  void setCleanImpl() const
914  {
915  is_dirty_.template reset<which>();
916  }
917 
918  template <int which>
919  bool isDirtyImpl() const
920  {
921  return is_dirty_.template test<which>();
922  }
923 };
924 
925  // DecoratorImpl implement the functionality of Decorator below
926 template <class A, unsigned CurrentPass, bool allowRuntimeActivation, unsigned WorkPass=A::workInPass>
927 struct DecoratorImpl
928 {
929  template <class T>
930  static void exec(A & a, T const & t)
931  {}
932 
933  template <class T>
934  static void exec(A & a, T const & t, double weight)
935  {}
936 };
937 
938 template <class A, unsigned CurrentPass>
939 struct DecoratorImpl<A, CurrentPass, false, CurrentPass>
940 {
941  template <class T>
942  static void exec(A & a, T const & t)
943  {
944  a.update(t);
945  }
946 
947  template <class T>
948  static void exec(A & a, T const & t, double weight)
949  {
950  a.update(t, weight);
951  }
952 
953  static typename A::result_type get(A const & a)
954  {
955  return a();
956  }
957 
958  static void merge(A & a, A const & o)
959  {
960  a += o;
961  }
962 
963  template <class T>
964  static void resize(A & a, T const & t)
965  {
966  a.reshape(t);
967  }
968 
969  static void applyHistogramOptions(A & a, HistogramOptions const & options)
970  {
971  ApplyHistogramOptions<typename A::Tag>::exec(a, options);
972  }
973 
974  static unsigned int passesRequired()
975  {
976  static const unsigned int A_workInPass = A::workInPass;
977  return std::max(A_workInPass, A::InternalBaseType::passesRequired());
978  }
979 };
980 
981 template <class A, unsigned CurrentPass>
982 struct DecoratorImpl<A, CurrentPass, true, CurrentPass>
983 {
984  static bool isActive(A const & a)
985  {
986  return A::isActiveImpl(getAccumulator<AccumulatorEnd>(a).active_accumulators_);
987  }
988 
989  template <class T>
990  static void exec(A & a, T const & t)
991  {
992  if(isActive(a))
993  a.update(t);
994  }
995 
996  template <class T>
997  static void exec(A & a, T const & t, double weight)
998  {
999  if(isActive(a))
1000  a.update(t, weight);
1001  }
1002 
1003  static typename A::result_type get(A const & a)
1004  {
1005  static const std::string message = std::string("get(accumulator): attempt to access inactive statistic '") +
1006  typeid(typename A::Tag).name() + "'.";
1007  vigra_precondition(isActive(a), message);
1008  return a();
1009  }
1010 
1011  static void merge(A & a, A const & o)
1012  {
1013  if(isActive(a))
1014  a += o;
1015  }
1016 
1017  template <class T>
1018  static void resize(A & a, T const & t)
1019  {
1020  if(isActive(a))
1021  a.reshape(t);
1022  }
1023 
1024  static void applyHistogramOptions(A & a, HistogramOptions const & options)
1025  {
1026  if(isActive(a))
1027  ApplyHistogramOptions<typename A::Tag>::exec(a, options);
1028  }
1029 
1030  template <class ActiveFlags>
1031  static unsigned int passesRequired(ActiveFlags const & flags)
1032  {
1033  static const unsigned int A_workInPass = A::workInPass;
1034  return A::isActiveImpl(flags)
1035  ? std::max(A_workInPass, A::InternalBaseType::passesRequired(flags))
1036  : A::InternalBaseType::passesRequired(flags);
1037  }
1038 };
1039 
1040  // Generic reshape function (expands to a no-op when T has fixed shape, and to
1041  // the appropriate specialized call otherwise). Shape is an instance of MultiArrayShape<N>::type.
1042 template <class T, class Shape>
1043 void reshapeImpl(T &, Shape const &)
1044 {}
1045 
1046 template <class T, class Shape, class Initial>
1047 void reshapeImpl(T &, Shape const &, Initial const & = T())
1048 {}
1049 
1050 template <unsigned int N, class T, class Alloc, class Shape>
1051 void reshapeImpl(MultiArray<N, T, Alloc> & a, Shape const & s, T const & initial = T())
1052 {
1053  MultiArray<N, T, Alloc>(s, initial).swap(a);
1054 }
1055 
1056 template <class T, class Alloc, class Shape>
1057 void reshapeImpl(Matrix<T, Alloc> & a, Shape const & s, T const & initial = T())
1058 {
1059  Matrix<T, Alloc>(s, initial).swap(a);
1060 }
1061 
1062  // generic functions to create suitable shape objects from various input data types
1063 template <unsigned int N, class T, class Stride>
1064 inline typename MultiArrayShape<N>::type
1065 shapeOf(MultiArrayView<N, T, Stride> const & a)
1066 {
1067  return a.shape();
1068 }
1069 
1070 template <class T, int N>
1071 inline Shape1
1072 shapeOf(TinyVector<T, N> const &)
1073 {
1074  return Shape1(N);
1075 }
1076 
1077 template <class T, class NEXT>
1078 inline CoupledHandle<T, NEXT> const &
1079 shapeOf(CoupledHandle<T, NEXT> const & t)
1080 {
1081  return t;
1082 }
1083 
1084 #define VIGRA_SHAPE_OF(type) \
1085 inline Shape1 \
1086 shapeOf(type) \
1087 { \
1088  return Shape1(1); \
1089 }
1090 
1091 VIGRA_SHAPE_OF(unsigned char)
1092 VIGRA_SHAPE_OF(signed char)
1093 VIGRA_SHAPE_OF(unsigned short)
1094 VIGRA_SHAPE_OF(short)
1095 VIGRA_SHAPE_OF(unsigned int)
1096 VIGRA_SHAPE_OF(int)
1097 VIGRA_SHAPE_OF(unsigned long)
1098 VIGRA_SHAPE_OF(long)
1099 VIGRA_SHAPE_OF(unsigned long long)
1100 VIGRA_SHAPE_OF(long long)
1101 VIGRA_SHAPE_OF(float)
1102 VIGRA_SHAPE_OF(double)
1103 VIGRA_SHAPE_OF(long double)
1104 
1105 #undef VIGRA_SHAPE_OF
1106 
1107  // LabelDispatch is only used in AccumulatorChainArrays and has the following functionalities:
1108  // * hold an accumulator chain for global statistics
1109  // * hold an array of accumulator chains (one per region) for region statistics
1110  // * forward data to the appropriate chains
1111  // * allocate the region array with appropriate size
1112  // * store and forward activation requests
1113  // * compute required number of passes as maximum from global and region accumulators
1114 template <class T, class GlobalAccumulators, class RegionAccumulators>
1115 struct LabelDispatch
1116 {
1117  typedef LabelDispatchTag Tag;
1118  typedef GlobalAccumulators GlobalAccumulatorChain;
1119  typedef RegionAccumulators RegionAccumulatorChain;
1120  typedef typename LookupTag<AccumulatorEnd, RegionAccumulatorChain>::type::AccumulatorFlags ActiveFlagsType;
1121  typedef ArrayVector<RegionAccumulatorChain> RegionAccumulatorArray;
1122 
1123  typedef LabelDispatch type;
1124  typedef LabelDispatch & reference;
1125  typedef LabelDispatch const & const_reference;
1126  typedef GlobalAccumulatorChain InternalBaseType;
1127 
1128  typedef T const & argument_type;
1129  typedef argument_type first_argument_type;
1130  typedef double second_argument_type;
1131  typedef RegionAccumulatorChain & result_type;
1132 
1133  static const int index = GlobalAccumulatorChain::index + 1;
1134 
1135  GlobalAccumulatorChain next_;
1136  RegionAccumulatorArray regions_;
1137  HistogramOptions region_histogram_options_;
1138  MultiArrayIndex ignore_label_;
1139  ActiveFlagsType active_region_accumulators_;
1140 
1141  template <class IndexDefinition, class TagFound=typename IndexDefinition::Tag>
1142  struct LabelIndexSelector
1143  {
1144  static const int value = 2; // default: CoupledHandle holds labels at index 2
1145 
1146  template <class U, class NEXT>
1147  static MultiArrayIndex exec(CoupledHandle<U, NEXT> const & t)
1148  {
1149  return (MultiArrayIndex)get<value>(t);
1150  }
1151  };
1152 
1153  template <class IndexDefinition>
1154  struct LabelIndexSelector<IndexDefinition, LabelArgTag>
1155  {
1156  static const int value = IndexDefinition::value;
1157 
1158  template <class U, class NEXT>
1159  static MultiArrayIndex exec(CoupledHandle<U, NEXT> const & t)
1160  {
1161  return (MultiArrayIndex)get<value>(t);
1162  }
1163  };
1164 
1165  template <class TAG>
1166  struct ActivateImpl
1167  {
1168  typedef typename LookupTag<TAG, type>::type TargetAccumulator;
1169 
1170  static void activate(GlobalAccumulatorChain & globals, RegionAccumulatorArray & regions,
1171  ActiveFlagsType & flags)
1172  {
1173  TargetAccumulator::template activateImpl<LabelDispatch>(
1174  flags, getAccumulator<AccumulatorEnd>(globals).active_accumulators_);
1175  for(unsigned int k=0; k<regions.size(); ++k)
1176  getAccumulator<AccumulatorEnd>(regions[k]).active_accumulators_ = flags;
1177  }
1178 
1179  static bool isActive(GlobalAccumulatorChain const &, ActiveFlagsType const & flags)
1180  {
1181  return TargetAccumulator::isActiveImpl(flags);
1182  }
1183  };
1184 
1185  template <class TAG>
1186  struct ActivateImpl<Global<TAG> >
1187  {
1188  static void activate(GlobalAccumulatorChain & globals, RegionAccumulatorArray &, ActiveFlagsType &)
1189  {
1190  LookupTag<TAG, GlobalAccumulatorChain>::type::activateImpl(getAccumulator<AccumulatorEnd>(globals).active_accumulators_);
1191  }
1192 
1193  static bool isActive(GlobalAccumulatorChain const & globals, ActiveFlagsType const &)
1194  {
1195  return LookupTag<TAG, GlobalAccumulatorChain>::type::isActiveImpl(getAccumulator<AccumulatorEnd>(globals).active_accumulators_);
1196  }
1197  };
1198 
1199  template <int INDEX>
1200  struct ActivateImpl<LabelArg<INDEX> >
1201  {
1202  static void activate(GlobalAccumulatorChain &, RegionAccumulatorArray &, ActiveFlagsType &)
1203  {}
1204 
1205  static bool isActive(GlobalAccumulatorChain const & globals, ActiveFlagsType const &)
1206  {
1207  return getAccumulator<LabelArg<INDEX> >(globals).isActive();
1208  }
1209  };
1210 
1211  typedef typename LookupTag<LabelArgTag, GlobalAccumulatorChain>::type FindLabelIndex;
1212 
1213  LabelDispatch()
1214  : next_(),
1215  regions_(),
1216  region_histogram_options_(),
1217  ignore_label_(-1),
1218  active_region_accumulators_()
1219  {}
1220 
1221  LabelDispatch(LabelDispatch const & o)
1222  : next_(o.next_),
1223  regions_(o.regions_),
1224  region_histogram_options_(o.region_histogram_options_),
1225  ignore_label_(o.ignore_label_),
1226  active_region_accumulators_(o.active_region_accumulators_)
1227  {
1228  for(unsigned int k=0; k<regions_.size(); ++k)
1229  {
1230  getAccumulator<AccumulatorEnd>(regions_[k]).setGlobalAccumulator(&next_);
1231  }
1232  }
1233 
1234  MultiArrayIndex maxRegionLabel() const
1235  {
1236  return (MultiArrayIndex)regions_.size() - 1;
1237  }
1238 
1239  void setMaxRegionLabel(unsigned maxlabel)
1240  {
1241  if(maxRegionLabel() == (MultiArrayIndex)maxlabel)
1242  return;
1243  unsigned int oldSize = regions_.size();
1244  regions_.resize(maxlabel + 1);
1245  for(unsigned int k=oldSize; k<regions_.size(); ++k)
1246  {
1247  getAccumulator<AccumulatorEnd>(regions_[k]).setGlobalAccumulator(&next_);
1248  getAccumulator<AccumulatorEnd>(regions_[k]).active_accumulators_ = active_region_accumulators_;
1249  regions_[k].applyHistogramOptions(region_histogram_options_);
1250  }
1251  }
1252 
1253  void ignoreLabel(MultiArrayIndex l)
1254  {
1255  ignore_label_ = l;
1256  }
1257 
1258  void applyHistogramOptions(HistogramOptions const & options)
1259  {
1260  applyHistogramOptions(options, options);
1261  }
1262 
1263  void applyHistogramOptions(HistogramOptions const & regionoptions, HistogramOptions const & globaloptions)
1264  {
1265  region_histogram_options_ = regionoptions;
1266  next_.applyHistogramOptions(globaloptions);
1267  }
1268 
1269  template <class U>
1270  void resize(U const & t)
1271  {
1272  if(regions_.size() == 0)
1273  {
1274  static const int labelIndex = LabelIndexSelector<FindLabelIndex>::value;
1275  typedef typename CoupledHandleCast<labelIndex, T>::type LabelHandle;
1276  typedef typename LabelHandle::value_type LabelType;
1277  typedef MultiArrayView<LabelHandle::dimensions, LabelType, StridedArrayTag> LabelArray;
1278  LabelArray labelArray(t.shape(), cast<labelIndex>(t).strides(), const_cast<LabelType *>(cast<labelIndex>(t).ptr()));
1279 
1280  LabelType minimum, maximum;
1281  labelArray.minmax(&minimum, &maximum);
1282  setMaxRegionLabel(maximum);
1283  }
1284  next_.resize(t);
1285  // FIXME: only call resize when label k actually exists?
1286  for(unsigned int k=0; k<regions_.size(); ++k)
1287  regions_[k].resize(t);
1288  }
1289 
1290  template <unsigned N>
1291  void pass(T const & t)
1292  {
1293  if(LabelIndexSelector<FindLabelIndex>::exec(t) != ignore_label_)
1294  {
1295  next_.template pass<N>(t);
1296  regions_[LabelIndexSelector<FindLabelIndex>::exec(t)].template pass<N>(t);
1297  }
1298  }
1299 
1300  template <unsigned N>
1301  void pass(T const & t, double weight)
1302  {
1303  if(LabelIndexSelector<FindLabelIndex>::exec(t) != ignore_label_)
1304  {
1305  next_.template pass<N>(t, weight);
1306  regions_[LabelIndexSelector<FindLabelIndex>::exec(t)].template pass<N>(t, weight);
1307  }
1308  }
1309 
1310  static unsigned int passesRequired()
1311  {
1312  return std::max(GlobalAccumulatorChain::passesRequired(), RegionAccumulatorChain::passesRequired());
1313  }
1314 
1315  unsigned int passesRequiredDynamic() const
1316  {
1317  return std::max(GlobalAccumulatorChain::passesRequired(getAccumulator<AccumulatorEnd>(next_).active_accumulators_),
1318  RegionAccumulatorChain::passesRequired(active_region_accumulators_));
1319  }
1320 
1321  void reset()
1322  {
1323  next_.reset();
1324 
1325  active_region_accumulators_.clear();
1326  RegionAccumulatorArray().swap(regions_);
1327  // FIXME: or is it better to just reset the region accumulators?
1328  // for(unsigned int k=0; k<regions_.size(); ++k)
1329  // regions_[k].reset();
1330  }
1331 
1332  template <class TAG>
1333  void activate()
1334  {
1335  ActivateImpl<TAG>::activate(next_, regions_, active_region_accumulators_);
1336  }
1337 
1338  void activateAll()
1339  {
1340  getAccumulator<AccumulatorEnd>(next_).active_accumulators_.set();
1341  active_region_accumulators_.set();
1342  for(unsigned int k=0; k<regions_.size(); ++k)
1343  getAccumulator<AccumulatorEnd>(regions_[k]).active_accumulators_.set();
1344  }
1345 
1346  template <class TAG>
1347  bool isActive() const
1348  {
1349  return ActivateImpl<TAG>::isActive(next_, active_region_accumulators_);
1350  }
1351 
1352  void merge(LabelDispatch const & o)
1353  {
1354  for(unsigned int k=0; k<regions_.size(); ++k)
1355  regions_[k].merge(o.regions_[k]);
1356  next_.merge(o.next_);
1357  }
1358 
1359  void merge(unsigned i, unsigned j)
1360  {
1361  regions_[i].merge(regions_[j]);
1362  regions_[j].reset();
1363  getAccumulator<AccumulatorEnd>(regions_[j]).active_accumulators_ = active_region_accumulators_;
1364  }
1365 
1366  template <class ArrayLike>
1367  void merge(LabelDispatch const & o, ArrayLike const & labelMapping)
1368  {
1369  MultiArrayIndex newMaxLabel = std::max<MultiArrayIndex>(maxRegionLabel(), *argMax(labelMapping.begin(), labelMapping.end()));
1370  setMaxRegionLabel(newMaxLabel);
1371  for(unsigned int k=0; k<labelMapping.size(); ++k)
1372  regions_[labelMapping[k]].merge(o.regions_[k]);
1373  next_.merge(o.next_);
1374  }
1375 };
1376 
1377 template <class TargetTag, class TagList>
1378 struct FindNextTag;
1379 
1380 template <class TargetTag, class HEAD, class TAIL>
1381 struct FindNextTag<TargetTag, TypeList<HEAD, TAIL> >
1382 {
1383  typedef typename FindNextTag<TargetTag, TAIL>::type type;
1384 };
1385 
1386 template <class TargetTag, class TAIL>
1387 struct FindNextTag<TargetTag, TypeList<TargetTag, TAIL> >
1388 {
1389  typedef typename TAIL::Head type;
1390 };
1391 
1392 template <class TargetTag>
1393 struct FindNextTag<TargetTag, TypeList<TargetTag, void> >
1394 {
1395  typedef void type;
1396 };
1397 
1398 template <class TargetTag>
1399 struct FindNextTag<TargetTag, void>
1400 {
1401  typedef void type;
1402 };
1403 
1404  // AccumulatorFactory creates the decorator hierarchy for the given TAG and configuration CONFIG
1405 template <class TAG, class CONFIG, unsigned LEVEL=0>
1406 struct AccumulatorFactory
1407 {
1408  typedef typename FindNextTag<TAG, typename CONFIG::TagList>::type NextTag;
1409  typedef typename AccumulatorFactory<NextTag, CONFIG, LEVEL+1>::type NextType;
1410  typedef typename CONFIG::InputType InputType;
1411 
1412  template <class T>
1413  struct ConfigureTag
1414  {
1415  typedef TAG type;
1416  };
1417 
1418  // When InputType is a CoupledHandle, some tags need to be wrapped into
1419  // DataFromHandle<> and/or Weighted<> modifiers. The following code does
1420  // this when appropriate.
1421  template <class T, class NEXT>
1422  struct ConfigureTag<CoupledHandle<T, NEXT> >
1423  {
1424  typedef typename StandardizeTag<DataFromHandle<TAG> >::type WrappedTag;
1425  typedef typename IfBool<(!HasModifierPriority<WrappedTag, WeightingPriority>::value && ShouldBeWeighted<WrappedTag>::value),
1426  Weighted<WrappedTag>, WrappedTag>::type type;
1427  };
1428 
1429  typedef typename ConfigureTag<InputType>::type UseTag;
1430 
1431  // base class of the decorator hierarchy: default (possibly empty)
1432  // implementations of all members
1433  struct AccumulatorBase
1434  {
1435  typedef AccumulatorBase ThisType;
1436  typedef TAG Tag;
1437  typedef NextType InternalBaseType;
1438  typedef InputType input_type;
1439  typedef input_type const & argument_type;
1440  typedef argument_type first_argument_type;
1441  typedef double second_argument_type;
1442  typedef void result_type;
1443 
1444  static const unsigned int workInPass = 1;
1445  static const int index = InternalBaseType::index + 1;
1446 
1447  InternalBaseType next_;
1448 
1449  static std::string name()
1450  {
1451  return TAG::name();
1452  }
1453 
1454  template <class ActiveFlags>
1455  static void activateImpl(ActiveFlags & flags)
1456  {
1457  flags.template set<index>();
1458  typedef typename StandardizeDependencies<Tag>::type StdDeps;
1459  detail::ActivateDependencies<StdDeps>::template exec<ThisType>(flags);
1460  }
1461 
1462  template <class Accu, class ActiveFlags, class GlobalFlags>
1463  static void activateImpl(ActiveFlags & flags, GlobalFlags & gflags)
1464  {
1465  flags.template set<index>();
1466  typedef typename StandardizeDependencies<Tag>::type StdDeps;
1467  detail::ActivateDependencies<StdDeps>::template exec<Accu>(flags, gflags);
1468  }
1469 
1470  template <class ActiveFlags>
1471  static bool isActiveImpl(ActiveFlags & flags)
1472  {
1473  return flags.template test<index>();
1474  }
1475 
1476  void setDirty() const
1477  {
1478  next_.template setDirtyImpl<index>();
1479  }
1480 
1481  template <int INDEX>
1482  void setDirtyImpl() const
1483  {
1484  next_.template setDirtyImpl<INDEX>();
1485  }
1486 
1487  void setClean() const
1488  {
1489  next_.template setCleanImpl<index>();
1490  }
1491 
1492  template <int INDEX>
1493  void setCleanImpl() const
1494  {
1495  next_.template setCleanImpl<INDEX>();
1496  }
1497 
1498  bool isDirty() const
1499  {
1500  return next_.template isDirtyImpl<index>();
1501  }
1502 
1503  template <int INDEX>
1504  bool isDirtyImpl() const
1505  {
1506  return next_.template isDirtyImpl<INDEX>();
1507  }
1508 
1509  void reset()
1510  {}
1511 
1512  template <class Shape>
1513  void reshape(Shape const &)
1514  {}
1515 
1516  void operator+=(AccumulatorBase const &)
1517  {}
1518 
1519  template <class U>
1520  void update(U const &)
1521  {}
1522 
1523  template <class U>
1524  void update(U const &, double)
1525  {}
1526 
1527  template <class TargetTag>
1528  typename LookupDependency<TargetTag, ThisType>::result_type
1529  call_getDependency() const
1530  {
1531  return getDependency<TargetTag>(*this);
1532  }
1533  };
1534 
1535  // The middle class(es) of the decorator hierarchy implement the actual feature computation.
1536  typedef typename UseTag::template Impl<InputType, AccumulatorBase> AccumulatorImpl;
1537 
1538  // outer class of the decorator hierarchy. It has the following functionalities
1539  // * ensure that only active accumulators are called in a dynamic accumulator chain
1540  // * ensure that each accumulator is only called in its desired pass as defined in A::workInPass
1541  // * determine how many passes through the data are required
1542  struct Accumulator
1543  : public AccumulatorImpl
1544  {
1545  typedef Accumulator type;
1546  typedef Accumulator & reference;
1547  typedef Accumulator const & const_reference;
1548  typedef AccumulatorImpl A;
1549 
1550  static const unsigned int workInPass = A::workInPass;
1551  static const bool allowRuntimeActivation = CONFIG::allowRuntimeActivation;
1552 
1553  template <class T>
1554  void resize(T const & t)
1555  {
1556  this->next_.resize(t);
1557  DecoratorImpl<Accumulator, workInPass, allowRuntimeActivation>::resize(*this, t);
1558  }
1559 
1560  void reset()
1561  {
1562  this->next_.reset();
1563  A::reset();
1564  }
1565 
1566  typename A::result_type get() const
1567  {
1569  }
1570 
1571  template <unsigned N, class T>
1572  void pass(T const & t)
1573  {
1574  this->next_.template pass<N>(t);
1575  DecoratorImpl<Accumulator, N, allowRuntimeActivation>::exec(*this, t);
1576  }
1577 
1578  template <unsigned N, class T>
1579  void pass(T const & t, double weight)
1580  {
1581  this->next_.template pass<N>(t, weight);
1582  DecoratorImpl<Accumulator, N, allowRuntimeActivation>::exec(*this, t, weight);
1583  }
1584 
1585  void merge(Accumulator const & o)
1586  {
1587  DecoratorImpl<Accumulator, Accumulator::workInPass, allowRuntimeActivation>::merge(*this, o);
1588  this->next_.merge(o.next_);
1589  }
1590 
1591  void applyHistogramOptions(HistogramOptions const & options)
1592  {
1593  DecoratorImpl<Accumulator, workInPass, allowRuntimeActivation>::applyHistogramOptions(*this, options);
1594  this->next_.applyHistogramOptions(options);
1595  }
1596 
1597  static unsigned int passesRequired()
1598  {
1599  return DecoratorImpl<Accumulator, workInPass, allowRuntimeActivation>::passesRequired();
1600  }
1601 
1602  template <class ActiveFlags>
1603  static unsigned int passesRequired(ActiveFlags const & flags)
1604  {
1605  return DecoratorImpl<Accumulator, workInPass, allowRuntimeActivation>::passesRequired(flags);
1606  }
1607  };
1608 
1609  typedef Accumulator type;
1610 };
1611 
1612 template <class CONFIG, unsigned LEVEL>
1613 struct AccumulatorFactory<void, CONFIG, LEVEL>
1614 {
1615  typedef AccumulatorEndImpl<LEVEL, typename CONFIG::GlobalAccumulatorHandle> type;
1616 };
1617 
1618 struct InvalidGlobalAccumulatorHandle
1619 {
1620  typedef Error__Global_statistics_are_only_defined_for_AccumulatorChainArray type;
1621 
1622  InvalidGlobalAccumulatorHandle()
1623  : pointer_(0)
1624  {}
1625 
1626  type const * pointer_;
1627 };
1628 
1629  // helper classes to create an accumulator chain from a TypeList
1630  // if dynamic=true, a dynamic accumulator will be created
1631  // if dynamic=false, a plain accumulator will be created
1632 template <class T, class Selected, bool dynamic=false, class GlobalHandle=InvalidGlobalAccumulatorHandle>
1633 struct ConfigureAccumulatorChain
1634 #ifndef DOXYGEN
1635 : public ConfigureAccumulatorChain<T, typename AddDependencies<typename Selected::type>::type, dynamic>
1636 #endif
1637 {};
1638 
1639 template <class T, class HEAD, class TAIL, bool dynamic, class GlobalHandle>
1640 struct ConfigureAccumulatorChain<T, TypeList<HEAD, TAIL>, dynamic, GlobalHandle>
1641 {
1642  typedef TypeList<HEAD, TAIL> TagList;
1643  typedef T InputType;
1644  static const bool allowRuntimeActivation = dynamic;
1645  typedef GlobalHandle GlobalAccumulatorHandle;
1646 
1647  typedef typename AccumulatorFactory<HEAD, ConfigureAccumulatorChain>::type type;
1648 };
1649 
1650 template <class T, class Selected, bool dynamic=false>
1651 struct ConfigureAccumulatorChainArray
1652 #ifndef DOXYGEN
1653 : public ConfigureAccumulatorChainArray<T, typename AddDependencies<typename Selected::type>::type, dynamic>
1654 #endif
1655 {};
1656 
1657 template <class T, class HEAD, class TAIL, bool dynamic>
1658 struct ConfigureAccumulatorChainArray<T, TypeList<HEAD, TAIL>, dynamic>
1659 {
1660  typedef TypeList<HEAD, TAIL> TagList;
1661  typedef SeparateGlobalAndRegionTags<TagList> TagSeparator;
1662  typedef typename TagSeparator::GlobalTags GlobalTags;
1663  typedef typename TagSeparator::RegionTags RegionTags;
1664  typedef typename ConfigureAccumulatorChain<T, GlobalTags, dynamic>::type GlobalAccumulatorChain;
1665 
1666  struct GlobalAccumulatorHandle
1667  {
1668  typedef GlobalAccumulatorChain type;
1669 
1670  GlobalAccumulatorHandle()
1671  : pointer_(0)
1672  {}
1673 
1674  type const * pointer_;
1675  };
1676 
1677  typedef typename ConfigureAccumulatorChain<T, RegionTags, dynamic, GlobalAccumulatorHandle>::type RegionAccumulatorChain;
1678 
1679  typedef LabelDispatch<T, GlobalAccumulatorChain, RegionAccumulatorChain> type;
1680 };
1681 
1682 } // namespace detail
1683 
1684 /****************************************************************************/
1685 /* */
1686 /* accumulator chain */
1687 /* */
1688 /****************************************************************************/
1689 
1690 // Implement the high-level interface of an accumulator chain
1691 template <class T, class NEXT>
1692 class AccumulatorChainImpl
1693 {
1694  public:
1695  typedef NEXT InternalBaseType;
1696  typedef AccumulatorBegin Tag;
1697  typedef typename InternalBaseType::argument_type argument_type;
1698  typedef typename InternalBaseType::first_argument_type first_argument_type;
1699  typedef typename InternalBaseType::second_argument_type second_argument_type;
1700  typedef void value_type;
1701  typedef typename InternalBaseType::result_type result_type;
1702 
1703  static const int staticSize = InternalBaseType::index;
1704 
1705  InternalBaseType next_;
1706 
1707  /** \brief Current pass of the accumulator chain.
1708  */
1709  unsigned int current_pass_;
1710 
1711  AccumulatorChainImpl()
1712  : current_pass_(0)
1713  {}
1714 
1715  /** 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.
1716  */
1717  void setHistogramOptions(HistogramOptions const & options)
1718  {
1719  next_.applyHistogramOptions(options);
1720  }
1721 
1722 
1723  /** Set regional and global options for all histograms in the accumulator chain.
1724  */
1725  void setHistogramOptions(HistogramOptions const & regionoptions, HistogramOptions const & globaloptions)
1726  {
1727  next_.applyHistogramOptions(regionoptions, globaloptions);
1728  }
1729 
1730  /** Reset current_pass_ of the accumulator chain to 'reset_to_pass'.
1731  */
1732  void reset(unsigned int reset_to_pass = 0)
1733  {
1734  current_pass_ = reset_to_pass;
1735  if(reset_to_pass == 0)
1736  next_.reset();
1737  }
1738 
1739  template <unsigned N>
1740  void update(T const & t)
1741  {
1742  if(current_pass_ == N)
1743  {
1744  next_.template pass<N>(t);
1745  }
1746  else if(current_pass_ < N)
1747  {
1748  current_pass_ = N;
1749  if(N == 1)
1750  next_.resize(detail::shapeOf(t));
1751  next_.template pass<N>(t);
1752  }
1753  else
1754  {
1755  std::string message("AccumulatorChain::update(): cannot return to pass ");
1756  message << N << " after working on pass " << current_pass_ << ".";
1757  vigra_precondition(false, message);
1758  }
1759  }
1760 
1761  template <unsigned N>
1762  void update(T const & t, double weight)
1763  {
1764  if(current_pass_ == N)
1765  {
1766  next_.template pass<N>(t, weight);
1767  }
1768  else if(current_pass_ < N)
1769  {
1770  current_pass_ = N;
1771  if(N == 1)
1772  next_.resize(detail::shapeOf(t));
1773  next_.template pass<N>(t, weight);
1774  }
1775  else
1776  {
1777  std::string message("AccumulatorChain::update(): cannot return to pass ");
1778  message << N << " after working on pass " << current_pass_ << ".";
1779  vigra_precondition(false, message);
1780  }
1781  }
1782 
1783  /** Equivalent to merge(o) .
1784  */
1785  void operator+=(AccumulatorChainImpl const & o)
1786  {
1787  merge(o);
1788  }
1789 
1790  /** 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.
1791  */
1792  void merge(AccumulatorChainImpl const & o)
1793  {
1794  next_.merge(o.next_);
1795  }
1796 
1797  result_type operator()() const
1798  {
1799  return next_.get();
1800  }
1801 
1802  void operator()(T const & t)
1803  {
1804  update<1>(t);
1805  }
1806 
1807  void operator()(T const & t, double weight)
1808  {
1809  update<1>(t, weight);
1810  }
1811 
1812  void updatePass2(T const & t)
1813  {
1814  update<2>(t);
1815  }
1816 
1817  void updatePass2(T const & t, double weight)
1818  {
1819  update<2>(t, weight);
1820  }
1821 
1822  /** 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.
1823  */
1824  void updatePassN(T const & t, unsigned int N)
1825  {
1826  switch (N)
1827  {
1828  case 1: update<1>(t); break;
1829  case 2: update<2>(t); break;
1830  case 3: update<3>(t); break;
1831  case 4: update<4>(t); break;
1832  case 5: update<5>(t); break;
1833  default:
1834  vigra_precondition(false,
1835  "AccumulatorChain::updatePassN(): 0 < N < 6 required.");
1836  }
1837  }
1838 
1839  /** 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.
1840  */
1841  void updatePassN(T const & t, double weight, unsigned int N)
1842  {
1843  switch (N)
1844  {
1845  case 1: update<1>(t, weight); break;
1846  case 2: update<2>(t, weight); break;
1847  case 3: update<3>(t, weight); break;
1848  case 4: update<4>(t, weight); break;
1849  case 5: update<5>(t, weight); break;
1850  default:
1851  vigra_precondition(false,
1852  "AccumulatorChain::updatePassN(): 0 < N < 6 required.");
1853  }
1854  }
1855 
1856  /** Return the number of passes required to compute all statistics in the accumulator chain.
1857  */
1858  unsigned int passesRequired() const
1859  {
1860  return InternalBaseType::passesRequired();
1861  }
1862 };
1863 
1864 
1865 
1866  // Create an accumulator chain containing the Selected statistics and their dependencies.
1867 
1868 /** \brief Create an accumulator chain containing the selected statistics and their dependencies.
1869 
1870  AccumulatorChain is used to compute global statistics which have to be selected at compile time.
1871 
1872  The template parameters are as follows:
1873  - T: The input type
1874  - either element type of the data(e.g. double, int, RGBValue, ...)
1875  - or type of CoupledHandle (for simultaneous access to coordinates and/or weights)
1876  - Selected: statistics to be computed and index specifier for the CoupledHandle, wrapped with Select
1877 
1878  Usage:
1879  \code
1880  typedef double DataType;
1881  AccumulatorChain<DataType, Select<Variance, Mean, Minimum, ...> > accumulator;
1882  \endcode
1883 
1884  Usage, using CoupledHandle:
1885  \code
1886  const int dim = 3; //dimension of MultiArray
1887  typedef double DataType;
1888  typedef double WeightType;
1889  typedef vigra::CoupledIteratorType<dim, DataType, WeightType>::HandleType Handle;
1890  AccumulatorChain<Handle, Select<DataArg<1>, WeightArg<2>, Mean,...> > a;
1891  \endcode
1892 
1893  See \ref FeatureAccumulators for more information and examples of use.
1894  */
1895 template <class T, class Selected, bool dynamic=false>
1897 #ifndef DOXYGEN // hide AccumulatorChainImpl from documentation
1898 : public AccumulatorChainImpl<T, typename detail::ConfigureAccumulatorChain<T, Selected, dynamic>::type>
1899 #endif
1900 {
1901  public:
1902  // \brief TypeList of Tags in the accumulator chain (?).
1903  typedef typename detail::ConfigureAccumulatorChain<T, Selected, dynamic>::TagList AccumulatorTags;
1904 
1905  /** Before having seen data (current_pass_==0), the shape of the data can be changed... (?)
1906  */
1907  template <class U, int N>
1908  void reshape(TinyVector<U, N> const & s)
1909  {
1910  vigra_precondition(this->current_pass_ == 0,
1911  "AccumulatorChain::reshape(): cannot reshape after seeing data. Call AccumulatorChain::reset() first.");
1912  this->next_.resize(s);
1913  this->current_pass_ = 1;
1914  }
1915 
1916  /** Return the names of all tags in the accumulator chain (selected statistics and their dependencies).
1917  */
1919  {
1920  static const ArrayVector<std::string> n = collectTagNames();
1921  return n;
1922  }
1923 
1924 
1925 #ifdef DOXYGEN // hide AccumulatorChainImpl from documentation
1926 
1927  /** 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.
1928  */
1929  void setHistogramOptions(HistogramOptions const & options);
1930 
1931  /** Reset current_pass_ of the accumulator chain to 'reset_to_pass'. */
1932  void reset(unsigned int reset_to_pass = 0);
1933 
1934  /** Equivalent to merge(o) . */
1935  void operator+=(AccumulatorChainImpl const & o);
1936 
1937  /** 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.
1938  */
1939  void merge(AccumulatorChainImpl const & o);
1940 
1941  /** 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.
1942  */
1943  void updatePassN(T const & t, unsigned int N);
1944 
1945  /** 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.
1946  */
1947  void updatePassN(T const & t, double weight, unsigned int N);
1948 
1949  /** Return the number of passes required to compute all statistics in the accumulator chain.
1950  */
1951  unsigned int passesRequired() const;
1952 
1953 #endif
1954 
1955  private:
1956  static ArrayVector<std::string> collectTagNames()
1957  {
1959  detail::CollectAccumulatorNames<AccumulatorTags>::exec(n);
1960  std::sort(n.begin(), n.end());
1961  return n;
1962  }
1963 };
1964 
1965 
1966  // Create a dynamic accumulator chain containing the Selected statistics and their dependencies.
1967  // Statistics will only be computed if activate<Tag>() is called at runtime.
1968 /** \brief Create a dynamic accumulator chain containing the selected statistics and their dependencies.
1969 
1970  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).
1971 
1972  The template parameters are as follows:
1973  - T: The input type
1974  - either element type of the data(e.g. double, int, RGBValue, ...)
1975  - or type of CoupledHandle (for access to coordinates and/or weights)
1976  - Selected: statistics to be computed and index specifier for the CoupledHandle, wrapped with Select
1977 
1978  Usage:
1979  \code
1980  typedef double DataType;
1981  DynamicAccumulatorChain<DataType, Select<Variance, Mean, Minimum, ...> > accumulator;
1982  \endcode
1983 
1984  Usage, using CoupledHandle:
1985  \code
1986  const int dim = 3; //dimension of MultiArray
1987  typedef double DataType;
1988  typedef double WeightType;
1989  typedef vigra::CoupledIteratorType<dim, DataType, WeightType>::HandleType Handle;
1990  DynamicAccumulatorChain<Handle, Select<DataArg<1>, WeightArg<2>, Mean,...> > a;
1991  \endcode
1992 
1993  See \ref FeatureAccumulators for more information and examples of use.
1994  */
1995 template <class T, class Selected>
1997 : public AccumulatorChain<T, Selected, true>
1998 {
1999  public:
2000  typedef typename AccumulatorChain<T, Selected, true>::InternalBaseType InternalBaseType;
2001  typedef typename DynamicAccumulatorChain::AccumulatorTags AccumulatorTags;
2002 
2003  /** Activate statistic 'tag'. Alias names are not recognized. If the statistic is not in the accumulator chain a PreconditionViolation is thrown.
2004  */
2005  void activate(std::string tag)
2006  {
2007  vigra_precondition(activateImpl(tag),
2008  std::string("DynamicAccumulatorChain::activate(): Tag '") + tag + "' not found.");
2009  }
2010 
2011  /** %activate<TAG>() activates statistic 'TAG'. If the statistic is not in the accumulator chain it is ignored. (?)
2012  */
2013  template <class TAG>
2014  void activate()
2015  {
2016  LookupTag<TAG, DynamicAccumulatorChain>::type::activateImpl(getAccumulator<AccumulatorEnd>(*this).active_accumulators_);
2017  }
2018 
2019  /** Activate all statistics in the accumulator chain.
2020  */
2022  {
2023  getAccumulator<AccumulatorEnd>(*this).active_accumulators_.set();
2024  }
2025  /** 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.)
2026  */
2027  bool isActive(std::string tag) const
2028  {
2029  detail::TagIsActive_Visitor v;
2030  vigra_precondition(isActiveImpl(tag, v),
2031  std::string("DynamicAccumulatorChain::isActive(): Tag '") + tag + "' not found.");
2032  return v.result;
2033  }
2034 
2035  /** %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. (?)
2036  */
2037  template <class TAG>
2038  bool isActive() const
2039  {
2040  return LookupTag<TAG, DynamicAccumulatorChain>::type::isActiveImpl(getAccumulator<AccumulatorEnd>(*this).active_accumulators_);
2041  }
2042 
2043  /** Return names of all statistics in the accumulator chain that are active.
2044  */
2046  {
2048  for(unsigned k=0; k<DynamicAccumulatorChain::tagNames().size(); ++k)
2050  res.push_back(DynamicAccumulatorChain::tagNames()[k]);
2051  return res;
2052  }
2053 
2054  /** Return number of passes required to compute the active statistics in the accumulator chain.
2055  */
2056  unsigned int passesRequired() const
2057  {
2058  return InternalBaseType::passesRequired(getAccumulator<AccumulatorEnd>(*this).active_accumulators_);
2059  }
2060 
2061  protected:
2062 
2063  bool activateImpl(std::string tag)
2064  {
2065  return detail::ApplyVisitorToTag<AccumulatorTags>::exec(*this,
2066  normalizeString(tag), detail::ActivateTag_Visitor());
2067  }
2068 
2069  bool isActiveImpl(std::string tag, detail::TagIsActive_Visitor & v) const
2070  {
2071  return detail::ApplyVisitorToTag<AccumulatorTags>::exec(*this, normalizeString(tag), v);
2072  }
2073 };
2074 
2075 /** \brief Create an array of accumulator chains containing the selected per-region and global statistics and their dependencies.
2076 
2077  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).
2078 
2079  The template parameters are as follows:
2080  - T: The input type, type of CoupledHandle (for access to coordinates, labels and weights)
2081  - Selected: statistics to be computed and index specifier for the CoupledHandle, wrapped with Select
2082 
2083  Usage:
2084  \code
2085  const int dim = 3; //dimension of MultiArray
2086  typedef double DataType;
2087  typedef double WeightType;
2088  typedef unsigned int LabelType;
2089  typedef vigra::CoupledIteratorType<dim, DataType, WeightType, LabelType>::HandleType Handle;
2090  AccumulatorChainArray<Handle, Select<DataArg<1>, WeightArg<2>, LabelArg<3>, Mean, Variance, ...> > a;
2091  \endcode
2092 
2093  See \ref FeatureAccumulators for more information and examples of use.
2094 */
2095 template <class T, class Selected, bool dynamic=false>
2097 #ifndef DOXYGEN //hide AccumulatorChainImpl vom documentation
2098 : public AccumulatorChainImpl<T, typename detail::ConfigureAccumulatorChainArray<T, Selected, dynamic>::type>
2099 #endif
2100 {
2101  public:
2102  typedef typename detail::ConfigureAccumulatorChainArray<T, Selected, dynamic> Creator;
2103  typedef typename Creator::TagList AccumulatorTags;
2104  typedef typename Creator::GlobalTags GlobalTags;
2105  typedef typename Creator::RegionTags RegionTags;
2106 
2107  /** Statistics will not be computed for label l. Note that only one label can be ignored.
2108  */
2110  {
2111  this->next_.ignoreLabel(l);
2112  }
2113 
2114  /** Set the maximum region label (e.g. for merging two accumulator chains).
2115  */
2116  void setMaxRegionLabel(unsigned label)
2117  {
2118  this->next_.setMaxRegionLabel(label);
2119  }
2120 
2121  /** %Maximum region label. (equal to regionCount() - 1)
2122  */
2124  {
2125  return this->next_.maxRegionLabel();
2126  }
2127 
2128  /** Number of Regions. (equal to maxRegionLabel() + 1)
2129  */
2130  unsigned int regionCount() const
2131  {
2132  return this->next_.regions_.size();
2133  }
2134 
2135  /** Merge region i with region j.
2136  */
2137  void merge(unsigned i, unsigned j)
2138  {
2139  vigra_precondition(i <= maxRegionLabel() && j <= maxRegionLabel(),
2140  "AccumulatorChainArray::merge(): region labels out of range.");
2141  this->next_.merge(i, j);
2142  }
2143 
2144  /** Merge with accumulator chain o. maxRegionLabel() of the two accumulators must be equal.
2145  */
2147  {
2148  if(maxRegionLabel() == -1)
2150  vigra_precondition(maxRegionLabel() == o.maxRegionLabel(),
2151  "AccumulatorChainArray::merge(): maxRegionLabel must be equal.");
2152  this->next_.merge(o.next_);
2153  }
2154 
2155  /** 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().
2156  */
2157  template <class ArrayLike>
2158  void merge(AccumulatorChainArray const & o, ArrayLike const & labelMapping)
2159  {
2160  vigra_precondition(labelMapping.size() == o.regionCount(),
2161  "AccumulatorChainArray::merge(): labelMapping.size() must match regionCount() of RHS.");
2162  this->next_.merge(o.next_, labelMapping);
2163  }
2164 
2165  /** Return names of all tags in the accumulator chain (selected statistics and their dependencies).
2166  */
2168  {
2169  static const ArrayVector<std::string> n = collectTagNames();
2170  return n;
2171  }
2172 
2173 
2174 #ifdef DOXYGEN // hide AccumulatorChainImpl from documentation
2175 
2176  /** \copydoc AccumulatorChain::setHistogramOptions(HistogramOptions const &) */
2177  void setHistogramOptions(HistogramOptions const & options);
2178 
2179  /** Set regional and global options for all histograms in the accumulator chain.
2180  */
2181  void setHistogramOptions(HistogramOptions const & regionoptions, HistogramOptions const & globaloptions);
2182 
2183  /** \copydoc AccumulatorChain::reset() */
2184  void reset(unsigned int reset_to_pass = 0);
2185 
2186  /** \copydoc AccumulatorChain::operator+=() */
2187  void operator+=(AccumulatorChainImpl const & o);
2188 
2189  /** \copydoc AccumulatorChain::updatePassN(T const &,unsigned int) */
2190  void updatePassN(T const & t, unsigned int N);
2191 
2192  /** \copydoc AccumulatorChain::updatePassN(T const &,double,unsigned int) */
2193  void updatePassN(T const & t, double weight, unsigned int N);
2194 
2195 #endif
2196 
2197 
2198 
2199  private:
2200  static ArrayVector<std::string> collectTagNames()
2201  {
2203  detail::CollectAccumulatorNames<AccumulatorTags>::exec(n);
2204  std::sort(n.begin(), n.end());
2205  return n;
2206  }
2207 };
2208 
2209 /** \brief Create an array of dynamic accumulator chains containing the selected per-region and global statistics and their dependencies.
2210 
2211 
2212  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).
2213 
2214  The template parameters are as follows:
2215  - T: The input type, type of CoupledHandle (for access to coordinates, labels and weights)
2216  - Selected: statistics to be computed and index specifier for the CoupledHandle, wrapped with Select
2217 
2218  Usage:
2219  \code
2220  const int dim = 3; //dimension of MultiArray
2221  typedef double DataType;
2222  typedef double WeightType;
2223  typedef unsigned int LabelType;
2224  typedef vigra::CoupledIteratorType<dim, DataType, WeightType, LabelType>::HandleType Handle;
2225  DynamicAccumulatorChainArray<Handle, Select<DataArg<1>, WeightArg<2>, LabelArg<3>, Mean, Variance, ...> > a;
2226  \endcode
2227 
2228  See \ref FeatureAccumulators for more information and examples of use.
2229 */
2230 template <class T, class Selected>
2232 : public AccumulatorChainArray<T, Selected, true>
2233 {
2234  public:
2235  typedef typename DynamicAccumulatorChainArray::AccumulatorTags AccumulatorTags;
2236 
2237  /** \copydoc DynamicAccumulatorChain::activate(std::string tag) */
2238  void activate(std::string tag)
2239  {
2240  vigra_precondition(activateImpl(tag),
2241  std::string("DynamicAccumulatorChainArray::activate(): Tag '") + tag + "' not found.");
2242  }
2243 
2244  /** \copydoc DynamicAccumulatorChain::activate() */
2245  template <class TAG>
2246  void activate()
2247  {
2248  this->next_.template activate<TAG>();
2249  }
2250 
2251  /** \copydoc DynamicAccumulatorChain::activateAll() */
2253  {
2254  this->next_.activateAll();
2255  }
2256 
2257  /** 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.)
2258  */
2259  bool isActive(std::string tag) const
2260  {
2261  detail::TagIsActive_Visitor v;
2262  vigra_precondition(isActiveImpl(tag, v),
2263  std::string("DynamicAccumulatorChainArray::isActive(): Tag '") + tag + "' not found.");
2264  return v.result;
2265  }
2266 
2267  /** %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. (?)
2268  */
2269  template <class TAG>
2270  bool isActive() const
2271  {
2272  return this->next_.template isActive<TAG>();
2273  }
2274 
2275  /** \copydoc DynamicAccumulatorChain::activeNames() */
2277  {
2279  for(unsigned k=0; k<DynamicAccumulatorChainArray::tagNames().size(); ++k)
2281  res.push_back(DynamicAccumulatorChainArray::tagNames()[k]);
2282  return res;
2283  }
2284 
2285  /** \copydoc DynamicAccumulatorChain::passesRequired() */
2286  unsigned int passesRequired() const
2287  {
2288  return this->next_.passesRequiredDynamic();
2289  }
2290 
2291  protected:
2292 
2293  bool activateImpl(std::string tag)
2294  {
2295  return detail::ApplyVisitorToTag<AccumulatorTags>::exec(this->next_,
2296  normalizeString(tag), detail::ActivateTag_Visitor());
2297  }
2298 
2299  bool isActiveImpl(std::string tag, detail::TagIsActive_Visitor & v) const
2300  {
2301  return detail::ApplyVisitorToTag<AccumulatorTags>::exec(this->next_, normalizeString(tag), v);
2302  }
2303 };
2304 
2305 /****************************************************************************/
2306 /* */
2307 /* generic access functions */
2308 /* */
2309 /****************************************************************************/
2310 
2311 template <class TAG>
2312 struct Error__Attempt_to_access_inactive_statistic;
2313 
2314 namespace detail {
2315 
2316  // accumulator lookup rules: find the accumulator that implements TAG
2317 
2318  // When A does not implement TAG, continue search in A::InternalBaseType.
2319 template <class TAG, class A, class FromTag=typename A::Tag>
2320 struct LookupTagImpl
2321 #ifndef DOXYGEN
2322 : public LookupTagImpl<TAG, typename A::InternalBaseType>
2323 #endif
2324 {};
2325 
2326  // 'const A' is treated like A, except that the reference member is now const.
2327 template <class TAG, class A, class FromTag>
2328 struct LookupTagImpl<TAG, A const, FromTag>
2329 : public LookupTagImpl<TAG, A>
2330 {
2331  typedef typename LookupTagImpl<TAG, A>::type const & reference;
2332  typedef typename LookupTagImpl<TAG, A>::type const * pointer;
2333 };
2334 
2335  // When A implements TAG, report its type and associated information.
2336 template <class TAG, class A>
2337 struct LookupTagImpl<TAG, A, TAG>
2338 {
2339  typedef TAG Tag;
2340  typedef A type;
2341  typedef A & reference;
2342  typedef A * pointer;
2343  typedef typename A::value_type value_type;
2344  typedef typename A::result_type result_type;
2345 };
2346 
2347  // Again, 'const A' is treated like A, except that the reference member is now const.
2348 template <class TAG, class A>
2349 struct LookupTagImpl<TAG, A const, TAG>
2350 : public LookupTagImpl<TAG, A, TAG>
2351 {
2352  typedef typename LookupTagImpl<TAG, A, TAG>::type const & reference;
2353  typedef typename LookupTagImpl<TAG, A, TAG>::type const * pointer;
2354 };
2355 
2356  // Recursion termination: when we end up in AccumulatorEnd without finding a
2357  // suitable A, we stop and report an error
2358 template <class TAG, class A>
2359 struct LookupTagImpl<TAG, A, AccumulatorEnd>
2360 {
2361  typedef TAG Tag;
2362  typedef A type;
2363  typedef A & reference;
2364  typedef A * pointer;
2365  typedef Error__Attempt_to_access_inactive_statistic<TAG> value_type;
2366  typedef Error__Attempt_to_access_inactive_statistic<TAG> result_type;
2367 };
2368 
2369  // ... except when we are actually looking for AccumulatorEnd
2370 template <class A>
2371 struct LookupTagImpl<AccumulatorEnd, A, AccumulatorEnd>
2372 {
2373  typedef AccumulatorEnd Tag;
2374  typedef A type;
2375  typedef A & reference;
2376  typedef A * pointer;
2377  typedef void value_type;
2378  typedef void result_type;
2379 };
2380 
2381  // ... or we are looking for a global statistic, in which case
2382  // we continue the serach via A::GlobalAccumulatorType, but remember that
2383  // we are actually looking for a global tag.
2384 template <class TAG, class A>
2385 struct LookupTagImpl<Global<TAG>, A, AccumulatorEnd>
2386 : public LookupTagImpl<TAG, typename A::GlobalAccumulatorType>
2387 {
2388  typedef Global<TAG> Tag;
2389 };
2390 
2391  // When we encounter the LabelDispatch accumulator, we continue the
2392  // search via LabelDispatch::RegionAccumulatorChain by default
2393 template <class TAG, class A>
2394 struct LookupTagImpl<TAG, A, LabelDispatchTag>
2395 : public LookupTagImpl<TAG, typename A::RegionAccumulatorChain>
2396 {};
2397 
2398  // ... except when we are looking for a global statistic, in which case
2399  // we continue via LabelDispatch::GlobalAccumulatorChain, but remember that
2400  // we are actually looking for a global tag.
2401 template <class TAG, class A>
2402 struct LookupTagImpl<Global<TAG>, A, LabelDispatchTag>
2403 : public LookupTagImpl<TAG, typename A::GlobalAccumulatorChain>
2404 {
2405  typedef Global<TAG> Tag;
2406 };
2407 
2408  // ... or we are looking for the LabelDispatch accumulator itself
2409 template <class A>
2410 struct LookupTagImpl<LabelDispatchTag, A, LabelDispatchTag>
2411 {
2412  typedef LabelDispatchTag Tag;
2413  typedef A type;
2414  typedef A & reference;
2415  typedef A * pointer;
2416  typedef void value_type;
2417  typedef void result_type;
2418 };
2419 
2420 } // namespace detail
2421 
2422  // Lookup the accumulator in the chain A that implements the given TAG.
2423 template <class Tag, class A>
2424 struct LookupTag
2425 : public detail::LookupTagImpl<typename StandardizeTag<Tag>::type, A>
2426 {};
2427 
2428  // Lookup the dependency TAG of the accumulator A.
2429  // This template ensures that dependencies are used with matching modifiers.
2430  // Specifically, if you search for Count as a dependency of Weighted<Mean>, the search
2431  // actually returns Weighted<Count>, wheras Count will be returned for plain Mean.
2432 template <class Tag, class A, class TargetTag>
2433 struct LookupDependency
2434 : public detail::LookupTagImpl<
2435  typename TransferModifiers<TargetTag, typename StandardizeTag<Tag>::type>::type, A>
2436 {};
2437 
2438 
2439 namespace detail {
2440 
2441  // CastImpl applies the same rules as LookupTagImpl, but returns a reference to an
2442  // accumulator instance rather than an accumulator type
2443 template <class Tag, class FromTag, class reference>
2444 struct CastImpl
2445 {
2446  template <class A>
2447  static reference exec(A & a)
2448  {
2449  return CastImpl<Tag, typename A::InternalBaseType::Tag, reference>::exec(a.next_);
2450  }
2451 
2452  template <class A>
2453  static reference exec(A & a, MultiArrayIndex label)
2454  {
2455  return CastImpl<Tag, typename A::InternalBaseType::Tag, reference>::exec(a.next_, label);
2456  }
2457 };
2458 
2459 template <class Tag, class reference>
2460 struct CastImpl<Tag, Tag, reference>
2461 {
2462  template <class A>
2463  static reference exec(A & a)
2464  {
2465  return const_cast<reference>(a);
2466  }
2467 
2468  template <class A>
2469  static reference exec(A & a, MultiArrayIndex)
2470  {
2471  vigra_precondition(false,
2472  "getAccumulator(): region accumulators can only be queried for AccumulatorChainArray.");
2473  return a;
2474  }
2475 };
2476 
2477 template <class Tag, class reference>
2478 struct CastImpl<Tag, AccumulatorEnd, reference>
2479 {
2480  template <class A>
2481  static reference exec(A & a)
2482  {
2483  return a;
2484  }
2485 
2486  template <class A>
2487  static reference exec(A & a, MultiArrayIndex)
2488  {
2489  return a;
2490  }
2491 };
2492 
2493 template <class Tag, class reference>
2494 struct CastImpl<Global<Tag>, AccumulatorEnd, reference>
2495 {
2496  template <class A>
2497  static reference exec(A & a)
2498  {
2499  return CastImpl<Tag, typename A::GlobalAccumulatorType::Tag, reference>::exec(*a.globalAccumulator_.pointer_);
2500  }
2501 };
2502 
2503 template <class reference>
2504 struct CastImpl<AccumulatorEnd, AccumulatorEnd, reference>
2505 {
2506  template <class A>
2507  static reference exec(A & a)
2508  {
2509  return a;
2510  }
2511 
2512  template <class A>
2513  static reference exec(A & a, MultiArrayIndex)
2514  {
2515  return a;
2516  }
2517 };
2518 
2519 template <class Tag, class reference>
2520 struct CastImpl<Tag, LabelDispatchTag, reference>
2521 {
2522  template <class A>
2523  static reference exec(A & a)
2524  {
2525  vigra_precondition(false,
2526  "getAccumulator(): a region label is required when a region accumulator is queried.");
2527  return CastImpl<Tag, typename A::RegionAccumulatorChain::Tag, reference>::exec(a.regions_[0]);
2528  }
2529 
2530  template <class A>
2531  static reference exec(A & a, MultiArrayIndex label)
2532  {
2533  return CastImpl<Tag, typename A::RegionAccumulatorChain::Tag, reference>::exec(a.regions_[label]);
2534  }
2535 };
2536 
2537 template <class Tag, class reference>
2538 struct CastImpl<Global<Tag>, LabelDispatchTag, reference>
2539 {
2540  template <class A>
2541  static reference exec(A & a)
2542  {
2543  return CastImpl<Tag, typename A::GlobalAccumulatorChain::Tag, reference>::exec(a.next_);
2544  }
2545 };
2546 
2547 template <class reference>
2548 struct CastImpl<LabelDispatchTag, LabelDispatchTag, reference>
2549 {
2550  template <class A>
2551  static reference exec(A & a)
2552  {
2553  return a;
2554  }
2555 };
2556 
2557 } // namespace detail
2558 
2559  // Get a reference to the accumulator TAG in the accumulator chain A
2560 /** 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
2561 Example of use (set options):
2562 \code
2563  vigra::MultiArray<2, double> data(...);
2564  typedef UserRangeHistogram<40> SomeHistogram; //binCount set at compile time
2565  typedef UserRangeHistogram<0> SomeHistogram2; // binCount must be set at run-time
2566  AccumulatorChain<DataType, Select<SomeHistogram, SomeHistogram2> > a;
2567 
2568  getAccumulator<SomeHistogram>(a).setMinMax(0.1, 0.9);
2569  getAccumulator<SomeHistogram2>(a).setMinMax(0.0, 1.0);
2570 
2571  extractFeatures(data.begin(), data.end(), a);
2572 \endcode
2573 
2574 Example of use (get information):
2575 \code
2576  vigra::MultiArray<2, double> data(...));
2577  AccumulatorChain<double, Select<Mean, Skewness> > a;
2578 
2579  std::cout << "passes required for all statistics: " << a.passesRequired() << std::endl; //skewness needs two passes
2580  std::cout << "passes required by Mean: " << getAccumulator<Mean>(a).passesRequired() << std::endl;
2581 \endcode
2582 See \ref FeatureAccumulators for more information about feature computation via accumulators.
2583 */
2584 template <class TAG, class A>
2585 inline typename LookupTag<TAG, A>::reference
2587 {
2588  typedef typename LookupTag<TAG, A>::Tag StandardizedTag;
2589  typedef typename LookupTag<TAG, A>::reference reference;
2590  return detail::CastImpl<StandardizedTag, typename A::Tag, reference>::exec(a);
2591 }
2592 
2593  // Get a reference to the accumulator TAG for region 'label' in the accumulator chain A
2594 /** Get a reference to the accumulator 'TAG' for region 'label' in the accumulator chain 'a'.
2595 */
2596 template <class TAG, class A>
2597 inline typename LookupTag<TAG, A>::reference
2599 {
2600  typedef typename LookupTag<TAG, A>::Tag StandardizedTag;
2601  typedef typename LookupTag<TAG, A>::reference reference;
2602  return detail::CastImpl<StandardizedTag, typename A::Tag, reference>::exec(a, label);
2603 }
2604 
2605  // get the result of the accumulator specified by TAG
2606 /** Get the result of the accumulator 'TAG' in the accumulator chain 'a'.\n
2607 Example of use:
2608 \code
2609  vigra::MultiArray<2, double> data(...);
2610  AccumulatorChain<DataType, Select<Variance, Mean, StdDev> > a;
2611  extractFeatures(data.begin(), data.end(), a);
2612  double mean = get<Mean>(a);
2613 \endcode
2614 See \ref FeatureAccumulators for more information about feature computation via accumulators.
2615 */
2616 template <class TAG, class A>
2617 inline typename LookupTag<TAG, A>::result_type
2618 get(A const & a)
2619 {
2620  return getAccumulator<TAG>(a).get();
2621 }
2622 
2623  // get the result of the accumulator TAG for region 'label'
2624 /** Get the result of the accumulator 'TAG' for region 'label' in the accumulator chain 'a'.\n
2625 Example of use:
2626 \code
2627  vigra::MultiArray<2, double> data(...);
2628  vigra::MultiArray<2, int> labels(...);
2629  typedef vigra::CoupledIteratorType<2, double, int>::type Iterator;
2630  typedef Iterator::value_type Handle;
2631 
2632  AccumulatorChainArray<Handle,
2633  Select<DataArg<1>, LabelArg<2>, Mean, Variance> > a;
2634 
2635  Iterator start = createCoupledIterator(data, labels);
2636  Iterator end = start.getEndIterator();
2637  extractFeatures(start,end,a);
2638 
2639  double mean_of_region_1 = get<Mean>(a,1);
2640  double mean_of_background = get<Mean>(a,0);
2641 \endcode
2642 See \ref FeatureAccumulators for more information about feature computation via accumulators.
2643 */
2644 template <class TAG, class A>
2645 inline typename LookupTag<TAG, A>::result_type
2646 get(A const & a, MultiArrayIndex label)
2647 {
2648  return getAccumulator<TAG>(a, label).get();
2649 }
2650 
2651  // Get the result of the accumulator specified by TAG without checking if the accumulator is active.
2652  // This must be used within an accumulator implementation to access dependencies because
2653  // it applies the approprate modifiers to the given TAG. It must not be used in other situations.
2654  // FIXME: is there a shorter name?
2655 template <class TAG, class A>
2656 inline typename LookupDependency<TAG, A>::result_type
2657 getDependency(A const & a)
2658 {
2659  typedef typename LookupDependency<TAG, A>::Tag StandardizedTag;
2660  typedef typename LookupDependency<TAG, A>::reference reference;
2661  return detail::CastImpl<StandardizedTag, typename A::Tag, reference>::exec(a)();
2662 }
2663 
2664  // activate the dynamic accumulator specified by Tag
2665 /** 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
2666 See \ref FeatureAccumulators for more information about feature computation via accumulators.
2667 */
2668 template <class Tag, class A>
2669 inline void
2670 activate(A & a)
2671 {
2672  a.template activate<Tag>();
2673 }
2674 
2675  // check if the dynamic accumulator specified by Tag is active
2676 /** 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
2677 See \ref FeatureAccumulators for more information about feature computation via accumulators.
2678 */
2679 template <class Tag, class A>
2680 inline bool
2681 isActive(A const & a)
2682 {
2683  return a.template isActive<Tag>();
2684 }
2685 
2686 /****************************************************************************/
2687 /* */
2688 /* generic loops */
2689 /* */
2690 /****************************************************************************/
2691 
2692 /** Generic loop to collect the statistics of the accumulator chain 'a' in as many passes over the data as necessary.\n
2693 
2694 Example of use:
2695 \code
2696  vigra::MultiArray<3, double> data(...);
2697  vigra::MultiArray<3, int> labels(...);
2698  typedef vigra::CoupledIteratorType<3, double, int>::type Iterator;
2699  typedef Iterator::value_type Handle;
2700 
2701  AccumulatorChainArray<Handle,
2702  Select<DataArg<1>, LabelArg<2>, Mean, Variance> > a;
2703 
2704  Iterator start = createCoupledIterator(data, labels);
2705  Iterator end = start.getEndIterator();
2706 
2707  extractFeatures(start,end,a);
2708 \endcode
2709 See \ref FeatureAccumulators for more information about feature computation via accumulators.
2710 */
2711 template <class ITERATOR, class ACCUMULATOR>
2712 void extractFeatures(ITERATOR start, ITERATOR end, ACCUMULATOR & a)
2713 {
2714  for(unsigned int k=1; k <= a.passesRequired(); ++k)
2715  for(ITERATOR i=start; i < end; ++i)
2716  a.updatePassN(*i, k);
2717 }
2718 
2719 /****************************************************************************/
2720 /* */
2721 /* AccumulatorResultTraits */
2722 /* */
2723 /****************************************************************************/
2724 
2725 template <class T>
2726 struct AccumulatorResultTraits
2727 {
2728  typedef T type;
2729  typedef T element_type;
2730  typedef double element_promote_type;
2731  typedef T MinmaxType;
2732  typedef element_promote_type SumType;
2733  typedef element_promote_type FlatCovarianceType;
2734  typedef element_promote_type CovarianceType;
2735 };
2736 
2737 template <class T, int N>
2738 struct AccumulatorResultTraits<TinyVector<T, N> >
2739 {
2740  typedef TinyVector<T, N> type;
2741  typedef T element_type;
2742  typedef double element_promote_type;
2743  typedef TinyVector<T, N> MinmaxType;
2744  typedef TinyVector<element_promote_type, N> SumType;
2745  typedef TinyVector<element_promote_type, N*(N+1)/2> FlatCovarianceType;
2746  typedef Matrix<element_promote_type> CovarianceType;
2747 };
2748 
2749 // (?) beign change
2750 template <class T, unsigned int RED_IDX, unsigned int GREEN_IDX, unsigned int BLUE_IDX>
2751 struct AccumulatorResultTraits<RGBValue<T, RED_IDX, GREEN_IDX, BLUE_IDX> >
2752 {
2753  typedef RGBValue<T> type;
2754  typedef T element_type;
2755  typedef double element_promote_type;
2756  typedef RGBValue<T> MinmaxType;
2757  typedef RGBValue<element_promote_type> SumType;
2758  typedef TinyVector<element_promote_type, 3*(3+1)/2> FlatCovarianceType;
2759  typedef Matrix<element_promote_type> CovarianceType;
2760 };
2761 // end change
2762 
2763 
2764 template <unsigned int N, class T, class Stride>
2765 struct AccumulatorResultTraits<MultiArrayView<N, T, Stride> >
2766 {
2767  typedef MultiArrayView<N, T, Stride> type;
2768  typedef T element_type;
2769  typedef double element_promote_type;
2770  typedef MultiArray<N, T> MinmaxType;
2771  typedef MultiArray<N, element_promote_type> SumType;
2772  typedef MultiArray<1, element_promote_type> FlatCovarianceType;
2773  typedef Matrix<element_promote_type> CovarianceType;
2774 };
2775 
2776 template <unsigned int N, class T, class Alloc>
2777 struct AccumulatorResultTraits<MultiArray<N, T, Alloc> >
2778 {
2779  typedef MultiArrayView<N, T, Alloc> type;
2780  typedef T element_type;
2781  typedef double element_promote_type;
2782  typedef MultiArray<N, T> MinmaxType;
2783  typedef MultiArray<N, element_promote_type> SumType;
2784  typedef MultiArray<1, element_promote_type> FlatCovarianceType;
2785  typedef Matrix<element_promote_type> CovarianceType;
2786 };
2787 
2788 /****************************************************************************/
2789 /* */
2790 /* modifier implementations */
2791 /* */
2792 /****************************************************************************/
2793 
2794 /** \brief Modifier. Compute statistic globally rather than per region.
2795 
2796 This modifier only works when labels are given (with (Dynamic)AccumulatorChainArray), in which case statistics are computed per-region by default.
2797 */
2798 template <class TAG>
2799 class Global
2800 {
2801  public:
2802  typedef typename StandardizeTag<TAG>::type TargetTag;
2803  typedef typename TargetTag::Dependencies Dependencies;
2804 
2805  static std::string const & name()
2806  {
2807  static const std::string n = std::string("Global<") + TargetTag::name() + " >";
2808  return n;
2809  }
2810 };
2811 
2812 /** \brief Specifies index of data in CoupledHandle.
2813 
2814  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)
2815 */
2816 template <int INDEX>
2817 class DataArg
2818 {
2819  public:
2820  typedef Select<> Dependencies;
2821 
2822  static std::string const & name()
2823  {
2824  static const std::string n = std::string("DataArg<") + asString(INDEX) + "> (internal)";
2825  return n;
2826  }
2827 
2828  template <class T, class BASE>
2829  struct Impl
2830  : public BASE
2831  {
2832  typedef DataArgTag Tag;
2833  typedef void value_type;
2834  typedef void result_type;
2835 
2836  static const int value = INDEX;
2837  static const unsigned int workInPass = 0;
2838  };
2839 };
2840 
2841 // Tags are automatically wrapped with DataFromHandle if CoupledHandle used
2842 template <class TAG>
2843 class DataFromHandle
2844 {
2845  public:
2846  typedef typename StandardizeTag<TAG>::type TargetTag;
2847  typedef typename TargetTag::Dependencies Dependencies;
2848 
2849  static std::string const & name()
2850  {
2851  static const std::string n = std::string("DataFromHandle<") + TargetTag::name() + " > (internal)";
2852  return n;
2853  }
2854 
2855  template <class IndexDefinition, class TagFound=typename IndexDefinition::Tag>
2856  struct DataIndexSelector
2857  {
2858  static const int value = 1; // default: CoupledHandle holds data at index 1
2859 
2860  template <class U, class NEXT>
2861  static typename CoupledHandleCast<value, CoupledHandle<U, NEXT> >::type::const_reference
2862  exec(CoupledHandle<U, NEXT> const & t)
2863  {
2864  return get<value>(t);
2865  }
2866  };
2867 
2868  template <class IndexDefinition>
2869  struct DataIndexSelector<IndexDefinition, DataArgTag>
2870  {
2871  static const int value = IndexDefinition::value;
2872 
2873  template <class U, class NEXT>
2874  static typename CoupledHandleCast<value, CoupledHandle<U, NEXT> >::type::const_reference
2875  exec(CoupledHandle<U, NEXT> const & t)
2876  {
2877  return get<value>(t);
2878  }
2879  };
2880 
2881  template <class T, class BASE>
2882  struct SelectInputType
2883  {
2884  typedef typename LookupTag<DataArgTag, BASE>::type FindDataIndex;
2885  typedef DataIndexSelector<FindDataIndex> DataIndex;
2886  typedef typename CoupledHandleCast<DataIndex::value, T>::type::value_type type;
2887  };
2888 
2889  template <class T, class BASE>
2890  struct Impl
2891  : public TargetTag::template Impl<typename SelectInputType<T, BASE>::type, BASE>
2892  {
2893  typedef SelectInputType<T, BASE> InputTypeSelector;
2894  typedef typename InputTypeSelector::DataIndex DataIndex;
2895  typedef typename InputTypeSelector::type input_type;
2896  typedef input_type const & argument_type;
2897  typedef argument_type first_argument_type;
2898 
2899  typedef typename TargetTag::template Impl<input_type, BASE> ImplType;
2900 
2901  using ImplType::reshape;
2902 
2903  template <class U, class NEXT>
2904  void reshape(CoupledHandle<U, NEXT> const & t)
2905  {
2906  ImplType::reshape(detail::shapeOf(DataIndex::exec(t)));
2907  }
2908 
2909  template <class U, class NEXT>
2910  void update(CoupledHandle<U, NEXT> const & t)
2911  {
2912  ImplType::update(DataIndex::exec(t));
2913  }
2914 
2915  template <class U, class NEXT>
2916  void update(CoupledHandle<U, NEXT> const & t, double weight)
2917  {
2918  ImplType::update(DataIndex::exec(t), weight);
2919  }
2920  };
2921 };
2922 
2923 /** \brief Modifier. Compute statistic from pixel coordinates rather than from pixel values.
2924 
2925  AccumulatorChain must be used with CoupledIterator in order to have access to pixel coordinates.
2926  */
2927 template <class TAG>
2928 class Coord
2929 {
2930  public:
2931  typedef typename StandardizeTag<TAG>::type TargetTag;
2932  typedef typename TargetTag::Dependencies Dependencies;
2933 
2934  static std::string const & name()
2935  {
2936  static const std::string n = std::string("Coord<") + TargetTag::name() + " >";
2937  return n;
2938  }
2939 
2940  template <class IndexDefinition, class TagFound=typename IndexDefinition::Tag>
2941  struct CoordIndexSelector
2942  {
2943  static const int value = 0; // default: CoupledHandle holds coordinates at index 0
2944 
2945  template <class U, class NEXT>
2946  static typename CoupledHandleCast<value, CoupledHandle<U, NEXT> >::type::const_reference
2947  exec(CoupledHandle<U, NEXT> const & t)
2948  {
2949  return get<value>(t);
2950  }
2951  };
2952 
2953  template <class IndexDefinition>
2954  struct CoordIndexSelector<IndexDefinition, CoordArgTag>
2955  {
2956  static const int value = IndexDefinition::value;
2957 
2958  template <class U, class NEXT>
2959  static typename CoupledHandleCast<value, CoupledHandle<U, NEXT> >::type::const_reference
2960  exec(CoupledHandle<U, NEXT> const & t)
2961  {
2962  return get<value>(t);
2963  }
2964  };
2965 
2966  template <class T, class BASE>
2967  struct SelectInputType
2968  {
2969  typedef typename LookupTag<CoordArgTag, BASE>::type FindDataIndex;
2970  typedef CoordIndexSelector<FindDataIndex> CoordIndex;
2971  typedef typename CoupledHandleCast<CoordIndex::value, T>::type::value_type type;
2972  };
2973 
2974  template <class T, class BASE>
2975  struct Impl
2976  : public TargetTag::template Impl<typename SelectInputType<T, BASE>::type, BASE>
2977  {
2978  typedef SelectInputType<T, BASE> InputTypeSelector;
2979  typedef typename InputTypeSelector::CoordIndex CoordIndex;
2980  typedef typename InputTypeSelector::type input_type;
2981  typedef input_type const & argument_type;
2982  typedef argument_type first_argument_type;
2983 
2984  typedef typename TargetTag::template Impl<input_type, BASE> ImplType;
2985 
2986  using ImplType::reshape;
2987 
2988  template <class U, class NEXT>
2989  void reshape(CoupledHandle<U, NEXT> const & t)
2990  {
2991  ImplType::reshape(detail::shapeOf(CoordIndex::exec(t)));
2992  }
2993 
2994  template <class U, class NEXT>
2995  void update(CoupledHandle<U, NEXT> const & t)
2996  {
2997  ImplType::update(CoordIndex::exec(t));
2998  }
2999 
3000  template <class U, class NEXT>
3001  void update(CoupledHandle<U, NEXT> const & t, double weight)
3002  {
3003  ImplType::update(CoordIndex::exec(t), weight);
3004  }
3005  };
3006 };
3007 
3008 /** \brief Specifies index of data in CoupledHandle.
3009 
3010  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.)
3011 */
3012 template <int INDEX>
3013 class WeightArg
3014 {
3015  public:
3016  typedef Select<> Dependencies;
3017 
3018  static std::string const & name()
3019  {
3020  static const std::string n = std::string("WeightArg<") + asString(INDEX) + "> (internal)";
3021  return n;
3022  }
3023 
3024  template <class T, class BASE>
3025  struct Impl
3026  : public BASE
3027  {
3028  typedef WeightArgTag Tag;
3029  typedef void value_type;
3030  typedef void result_type;
3031 
3032  static const int value = INDEX;
3033  static const unsigned int workInPass = 0;
3034  };
3035 };
3036 
3037 /** \brief Compute weighted version of the statistic.
3038 */
3039 template <class TAG>
3040 class Weighted
3041 {
3042  public:
3043  typedef typename StandardizeTag<TAG>::type TargetTag;
3044  typedef typename TargetTag::Dependencies Dependencies;
3045 
3046  static std::string const & name()
3047  {
3048  static const std::string n = std::string("Weighted<") + TargetTag::name() + " >";
3049  return n;
3050  }
3051 
3052  template <class IndexDefinition, class TagFound=typename IndexDefinition::Tag>
3053  struct WeightIndexSelector
3054  {
3055  template <class U, class NEXT>
3056  static double exec(CoupledHandle<U, NEXT> const & t)
3057  {
3058  return (double)*t; // default: CoupledHandle holds weights at the last (outermost) index
3059  }
3060  };
3061 
3062  template <class IndexDefinition>
3063  struct WeightIndexSelector<IndexDefinition, WeightArgTag>
3064  {
3065  template <class U, class NEXT>
3066  static double exec(CoupledHandle<U, NEXT> const & t)
3067  {
3068  return (double)get<IndexDefinition::value>(t);
3069  }
3070  };
3071 
3072  template <class T, class BASE>
3073  struct Impl
3074  : public TargetTag::template Impl<T, BASE>
3075  {
3076  typedef typename TargetTag::template Impl<T, BASE> ImplType;
3077 
3078  typedef typename LookupTag<WeightArgTag, BASE>::type FindWeightIndex;
3079 
3080  template <class U, class NEXT>
3081  void update(CoupledHandle<U, NEXT> const & t)
3082  {
3083  ImplType::update(t, WeightIndexSelector<FindWeightIndex>::exec(t));
3084  }
3085  };
3086 };
3087 
3088 // Centralize by subtracting the mean and cache the result
3089 class Centralize
3090 {
3091  public:
3092  typedef Select<Mean> Dependencies;
3093 
3094  static std::string const & name()
3095  {
3096  static const std::string n("Centralize (internal)");
3097  return n;
3098  }
3099 
3100  template <class U, class BASE>
3101  struct Impl
3102  : public BASE
3103  {
3104  static const unsigned int workInPass = 2;
3105 
3106  typedef typename AccumulatorResultTraits<U>::element_promote_type element_type;
3107  typedef typename AccumulatorResultTraits<U>::SumType value_type;
3108  typedef value_type const & result_type;
3109 
3110  mutable value_type value_;
3111 
3112  Impl()
3113  : value_() // call default constructor explicitly to ensure zero initialization
3114  {}
3115 
3116  void reset()
3117  {
3118  value_ = element_type();
3119  }
3120 
3121  template <class Shape>
3122  void reshape(Shape const & s)
3123  {
3124  detail::reshapeImpl(value_, s);
3125  }
3126 
3127  void update(U const & t) const
3128  {
3129  using namespace vigra::multi_math;
3130  value_ = t - getDependency<Mean>(*this);
3131  }
3132 
3133  void update(U const & t, double) const
3134  {
3135  update(t);
3136  }
3137 
3138  result_type operator()(U const & t) const
3139  {
3140  update(t);
3141  return value_;
3142  }
3143 
3144  result_type operator()() const
3145  {
3146  return value_;
3147  }
3148  };
3149 };
3150 
3151 /** \brief Modifier. Substract mean before computing statistic.
3152 
3153 Works in pass 2, %operator+=() not supported (merging not supported).
3154 */
3155 template <class TAG>
3156 class Central
3157 {
3158  public:
3159  typedef typename StandardizeTag<TAG>::type TargetTag;
3160  typedef Select<Centralize, typename TargetTag::Dependencies> Dependencies;
3161 
3162  static std::string const & name()
3163  {
3164  static const std::string n = std::string("Central<") + TargetTag::name() + " >";
3165  return n;
3166  }
3167 
3168  template <class U, class BASE>
3169  struct Impl
3170  : public TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE>
3171  {
3172  typedef typename TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE> ImplType;
3173 
3174  static const unsigned int workInPass = 2;
3175 
3176  void operator+=(Impl const & o)
3177  {
3178  vigra_precondition(false,
3179  "Central<...>::operator+=(): not supported.");
3180  }
3181 
3182  template <class T>
3183  void update(T const & t)
3184  {
3185  ImplType::update(getDependency<Centralize>(*this));
3186  }
3187 
3188  template <class T>
3189  void update(T const & t, double weight)
3190  {
3191  ImplType::update(getDependency<Centralize>(*this), weight);
3192  }
3193  };
3194 };
3195 
3196  // alternative implementation without caching
3197  //
3198 // template <class TAG>
3199 // class Central
3200 // {
3201  // public:
3202  // typedef typename StandardizeTag<TAG>::type TargetTag;
3203  // typedef TypeList<Mean, typename TransferModifiers<Central<TargetTag>, typename TargetTag::Dependencies::type>::type> Dependencies;
3204 
3205  // template <class U, class BASE>
3206  // struct Impl
3207  // : public TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE>
3208  // {
3209  // typedef typename TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE> ImplType;
3210 
3211  // static const unsigned int workInPass = 2;
3212 
3213  // void operator+=(Impl const & o)
3214  // {
3215  // vigra_precondition(false,
3216  // "Central<...>::operator+=(): not supported.");
3217  // }
3218 
3219  // template <class T>
3220  // void update(T const & t)
3221  // {
3222  // ImplType::update(t - getDependency<Mean>(*this));
3223  // }
3224 
3225  // template <class T>
3226  // void update(T const & t, double weight)
3227  // {
3228  // ImplType::update(t - getDependency<Mean>(*this), weight);
3229  // }
3230  // };
3231 // };
3232 
3233 
3234 class PrincipalProjection
3235 {
3236  public:
3237  typedef Select<Centralize, Principal<CoordinateSystem> > Dependencies;
3238 
3239  static std::string const & name()
3240  {
3241  static const std::string n("PrincipalProjection (internal)");
3242  return n;
3243  }
3244 
3245  template <class U, class BASE>
3246  struct Impl
3247  : public BASE
3248  {
3249  static const unsigned int workInPass = 2;
3250 
3251  typedef typename AccumulatorResultTraits<U>::element_promote_type element_type;
3252  typedef typename AccumulatorResultTraits<U>::SumType value_type;
3253  typedef value_type const & result_type;
3254 
3255  mutable value_type value_;
3256 
3257  Impl()
3258  : value_() // call default constructor explicitly to ensure zero initialization
3259  {}
3260 
3261  void reset()
3262  {
3263  value_ = element_type();
3264  }
3265 
3266  template <class Shape>
3267  void reshape(Shape const & s)
3268  {
3269  detail::reshapeImpl(value_, s);
3270  }
3271 
3272  void update(U const & t) const
3273  {
3274  for(unsigned int k=0; k<t.size(); ++k)
3275  {
3276  value_[k] = getDependency<Principal<CoordinateSystem> >(*this)(0, k)*getDependency<Centralize>(*this)[0];
3277  for(unsigned int d=1; d<t.size(); ++d)
3278  value_[k] += getDependency<Principal<CoordinateSystem> >(*this)(d, k)*getDependency<Centralize>(*this)[d];
3279  }
3280  }
3281 
3282  void update(U const & t, double) const
3283  {
3284  update(t);
3285  }
3286 
3287  result_type operator()(U const & t) const
3288  {
3289  getAccumulator<Centralize>(*this).update(t);
3290  update(t);
3291  return value_;
3292  }
3293 
3294  result_type operator()() const
3295  {
3296  return value_;
3297  }
3298  };
3299 };
3300 
3301 /** \brief Modifier. Project onto PCA eigenvectors.
3302 
3303  Works in pass 2, %operator+=() not supported (merging not supported).
3304 */
3305 template <class TAG>
3306 class Principal
3307 {
3308  public:
3309  typedef typename StandardizeTag<TAG>::type TargetTag;
3310  typedef Select<PrincipalProjection, typename TargetTag::Dependencies> Dependencies;
3311 
3312  static std::string const & name()
3313  {
3314  static const std::string n = std::string("Principal<") + TargetTag::name() + " >";
3315  return n;
3316  }
3317 
3318  template <class U, class BASE>
3319  struct Impl
3320  : public TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE>
3321  {
3322  typedef typename TargetTag::template Impl<typename AccumulatorResultTraits<U>::SumType, BASE> ImplType;
3323 
3324  static const unsigned int workInPass = 2;
3325 
3326  void operator+=(Impl const & o)
3327  {
3328  vigra_precondition(false,
3329  "Principal<...>::operator+=(): not supported.");
3330  }
3331 
3332  template <class T>
3333  void update(T const & t)
3334  {
3335  ImplType::update(getDependency<PrincipalProjection>(*this));
3336  }
3337 
3338  template <class T>
3339  void update(T const & t, double weight)
3340  {
3341  ImplType::update(getDependency<PrincipalProjection>(*this), weight);
3342  }
3343  };
3344 };
3345 
3346 /*
3347 important notes on modifiers:
3348  * upon accumulator creation, modifiers are reordered so that data preparation is innermost,
3349  and data access is outermost, e.g.:
3350  Coord<DivideByCount<Principal<PowerSum<2> > > >
3351  * modifiers are automatically transfered to dependencies as appropriate
3352  * modifiers for lookup (getAccumulator and get) of dependent accumulators are automatically adjusted
3353  * modifiers must adjust workInPass for the contained accumulator as appropriate
3354  * we may implement convenience versions of Select that apply a modifier to all
3355  contained tags at once
3356  * weighted accumulators have their own Count object when used together
3357  with unweighted ones (this is as yet untested - FIXME)
3358  * certain accumulators must remain unchanged when wrapped in certain modifiers:
3359  * Count: always except for Weighted<Count> and CoordWeighted<Count>
3360  * Sum: data preparation modifiers
3361  * FlatScatterMatrixImpl, CovarianceEigensystemImpl: Principal and Whitened
3362  * will it be useful to implement initPass<N>() or finalizePass<N>() ?
3363 */
3364 
3365 /****************************************************************************/
3366 /* */
3367 /* the actual accumulators */
3368 /* */
3369 /****************************************************************************/
3370 
3371 /** \brief Basic statistic. Identity matrix of appropriate size.
3372 */
3374 {
3375  public:
3376  typedef Select<> Dependencies;
3377 
3378  static std::string const & name()
3379  {
3380  static const std::string n("CoordinateSystem");
3381  return n;
3382  }
3383 
3384  template <class U, class BASE>
3385  struct Impl
3386  : public BASE
3387  {
3388  typedef double element_type;
3389  typedef Matrix<double> value_type;
3390  typedef value_type const & result_type;
3391 
3392  value_type value_;
3393 
3394  Impl()
3395  : value_() // call default constructor explicitly to ensure zero initialization
3396  {}
3397 
3398  void reset()
3399  {
3400  value_ = element_type();
3401  }
3402 
3403  template <class Shape>
3404  void reshape(Shape const & s)
3405  {
3406  detail::reshapeImpl(value_, s);
3407  }
3408 
3409  result_type operator()() const
3410  {
3411  return value_;
3412  }
3413  };
3414 };
3415 
3416 template <class BASE, class T,
3417  class ElementType=typename AccumulatorResultTraits<T>::element_promote_type,
3418  class SumType=typename AccumulatorResultTraits<T>::SumType>
3419 struct SumBaseImpl
3420 : public BASE
3421 {
3422  typedef ElementType element_type;
3423  typedef SumType value_type;
3424  typedef value_type const & result_type;
3425 
3426  value_type value_;
3427 
3428  SumBaseImpl()
3429  : value_() // call default constructor explicitly to ensure zero initialization
3430  {}
3431 
3432  void reset()
3433  {
3434  value_ = element_type();
3435  }
3436 
3437  template <class Shape>
3438  void reshape(Shape const & s)
3439  {
3440  detail::reshapeImpl(value_, s);
3441  }
3442 
3443  void operator+=(SumBaseImpl const & o)
3444  {
3445  value_ += o.value_;
3446  }
3447 
3448  result_type operator()() const
3449  {
3450  return value_;
3451  }
3452 };
3453 
3454 // Count
3455 template <>
3456 class PowerSum<0>
3457 {
3458  public:
3459  typedef Select<> Dependencies;
3460 
3461  static std::string const & name()
3462  {
3463  static const std::string n("PowerSum<0>");
3464  return n;
3465  }
3466 
3467  template <class T, class BASE>
3468  struct Impl
3469  : public SumBaseImpl<BASE, T, double, double>
3470  {
3471  void update(T const & t)
3472  {
3473  ++this->value_;
3474  }
3475 
3476  void update(T const & t, double weight)
3477  {
3478  this->value_ += weight;
3479  }
3480  };
3481 };
3482 
3483 // Sum
3484 template <>
3485 class PowerSum<1>
3486 {
3487  public:
3488  typedef Select<> Dependencies;
3489 
3490  static std::string const & name()
3491  {
3492  static const std::string n("PowerSum<1>");
3493  return n;
3494  }
3495 
3496  template <class U, class BASE>
3497  struct Impl
3498  : public SumBaseImpl<BASE, U>
3499  {
3500  void update(U const & t)
3501  {
3502  this->value_ += t;
3503  }
3504 
3505  void update(U const & t, double weight)
3506  {
3507  this->value_ += weight*t;
3508  }
3509  };
3510 };
3511 
3512 /** \brief Basic statistic. PowerSum<N> =@f$ \sum_i x_i^N @f$
3513 
3514  Works in pass 1, %operator+=() supported (merging supported).
3515 */
3516 template <unsigned N>
3517 class PowerSum
3518 {
3519  public:
3520  typedef Select<> Dependencies;
3521 
3522  static std::string const & name()
3523  {
3524  static const std::string n = std::string("PowerSum<") + asString(N) + ">";
3525  return n;
3526  }
3527 
3528  template <class U, class BASE>
3529  struct Impl
3530  : public SumBaseImpl<BASE, U>
3531  {
3532  void update(U const & t)
3533  {
3534  using namespace vigra::multi_math;
3535  this->value_ += pow(t, (int)N);
3536  }
3537 
3538  void update(U const & t, double weight)
3539  {
3540  using namespace vigra::multi_math;
3541  this->value_ += weight*pow(t, (int)N);
3542  }
3543  };
3544 };
3545 
3546 template <>
3547 class AbsPowerSum<1>
3548 {
3549  public:
3550  typedef Select<> Dependencies;
3551 
3552  static std::string const & name()
3553  {
3554  static const std::string n("AbsPowerSum<1>");
3555  return n;
3556  }
3557 
3558  template <class U, class BASE>
3559  struct Impl
3560  : public SumBaseImpl<BASE, U>
3561  {
3562  void update(U const & t)
3563  {
3564  using namespace vigra::multi_math;
3565  this->value_ += abs(t);
3566  }
3567 
3568  void update(U const & t, double weight)
3569  {
3570  using namespace vigra::multi_math;
3571  this->value_ += weight*abs(t);
3572  }
3573  };
3574 };
3575 
3576 /** \brief Basic statistic. AbsPowerSum<N> =@f$ \sum_i |x_i|^N @f$
3577 
3578  Works in pass 1, %operator+=() supported (merging supported).
3579 */
3580 template <unsigned N>
3581 class AbsPowerSum
3582 {
3583  public:
3584  typedef Select<> Dependencies;
3585 
3586  static std::string const & name()
3587  {
3588  static const std::string n = std::string("AbsPowerSum<") + asString(N) + ">";
3589  return n;
3590  }
3591 
3592  template <class U, class BASE>
3593  struct Impl
3594  : public SumBaseImpl<BASE, U>
3595  {
3596  void update(U const & t)
3597  {
3598  using namespace vigra::multi_math;
3599  this->value_ += pow(abs(t), (int)N);
3600  }
3601 
3602  void update(U const & t, double weight)
3603  {
3604  using namespace vigra::multi_math;
3605  this->value_ += weight*pow(abs(t), (int)N);
3606  }
3607  };
3608 };
3609 
3610 template <class BASE, class VALUE_TYPE, class U>
3611 struct CachedResultBase
3612 : public BASE
3613 {
3614  typedef typename AccumulatorResultTraits<U>::element_type element_type;
3615  typedef VALUE_TYPE value_type;
3616  typedef value_type const & result_type;
3617 
3618  mutable value_type value_;
3619 
3620  CachedResultBase()
3621  : value_() // call default constructor explicitly to ensure zero initialization
3622  {}
3623 
3624  void reset()
3625  {
3626  value_ = element_type();
3627  this->setClean();
3628  }
3629 
3630  template <class Shape>
3631  void reshape(Shape const & s)
3632  {
3633  detail::reshapeImpl(value_, s);
3634  }
3635 
3636  void operator+=(CachedResultBase const &)
3637  {
3638  this->setDirty();
3639  }
3640 
3641  void update(U const &)
3642  {
3643  this->setDirty();
3644  }
3645 
3646  void update(U const &, double)
3647  {
3648  this->setDirty();
3649  }
3650 };
3651 
3652 // cached Mean and Variance
3653 /** \brief Modifier. Divide statistic by Count: DivideByCount<TAG> = TAG / Count .
3654 */
3655 template <class TAG>
3656 class DivideByCount
3657 {
3658  public:
3659  typedef typename StandardizeTag<TAG>::type TargetTag;
3660  typedef Select<TargetTag, Count> Dependencies;
3661 
3662  static std::string const & name()
3663  {
3664  static const std::string n = std::string("DivideByCount<") + TargetTag::name() + " >";
3665  return n;
3666  }
3667 
3668  template <class U, class BASE>
3669  struct Impl
3670  : public CachedResultBase<BASE, typename LookupDependency<TargetTag, BASE>::value_type, U>
3671  {
3672  typedef typename CachedResultBase<BASE, typename LookupDependency<TargetTag, BASE>::value_type, U>::result_type result_type;
3673 
3674  result_type operator()() const
3675  {
3676  if(this->isDirty())
3677  {
3678  using namespace multi_math;
3679  this->value_ = getDependency<TargetTag>(*this) / getDependency<Count>(*this);
3680  this->setClean();
3681  }
3682  return this->value_;
3683  }
3684  };
3685 };
3686 
3687 // UnbiasedVariance
3688 /** \brief Modifier. Divide statistics by Count-1: DivideUnbiased<TAG> = TAG / (Count-1)
3689 */
3690 template <class TAG>
3691 class DivideUnbiased
3692 {
3693  public:
3694  typedef typename StandardizeTag<TAG>::type TargetTag;
3695  typedef Select<TargetTag, Count> Dependencies;
3696 
3697  static std::string const & name()
3698  {
3699  static const std::string n = std::string("DivideUnbiased<") + TargetTag::name() + " >";
3700  return n;
3701  }
3702 
3703  template <class U, class BASE>
3704  struct Impl
3705  : public BASE
3706  {
3707  typedef typename LookupDependency<TargetTag, BASE>::value_type value_type;
3708  typedef value_type result_type;
3709 
3710  result_type operator()() const
3711  {
3712  using namespace multi_math;
3713  return getDependency<TargetTag>(*this) / (getDependency<Count>(*this) - 1.0);
3714  }
3715  };
3716 };
3717 
3718 // RootMeanSquares and StdDev
3719 /** \brief Modifier. RootDivideByCount<TAG> = sqrt( TAG/Count )
3720 */
3721 template <class TAG>
3722 class RootDivideByCount
3723 {
3724  public:
3725  typedef typename StandardizeTag<DivideByCount<TAG> >::type TargetTag;
3726  typedef Select<TargetTag> Dependencies;
3727 
3728  static std::string const & name()
3729  {
3730  typedef typename StandardizeTag<TAG>::type InnerTag;
3731  static const std::string n = std::string("RootDivideByCount<") + InnerTag::name() + " >";
3732  return n;
3733  }
3734 
3735  template <class U, class BASE>
3736  struct Impl
3737  : public BASE
3738  {
3739  typedef typename LookupDependency<TargetTag, BASE>::value_type value_type;
3740  typedef value_type result_type;
3741 
3742  result_type operator()() const
3743  {
3744  using namespace multi_math;
3745  return sqrt(getDependency<TargetTag>(*this));
3746  }
3747  };
3748 };
3749 
3750 // UnbiasedStdDev
3751 /** \brief Modifier. RootDivideUnbiased<TAG> = sqrt( TAG / (Count-1) )
3752 */
3753 template <class TAG>
3754 class RootDivideUnbiased
3755 {
3756  public:
3757  typedef typename StandardizeTag<DivideUnbiased<TAG> >::type TargetTag;
3758  typedef Select<TargetTag> Dependencies;
3759 
3760  static std::string const & name()
3761  {
3762  typedef typename StandardizeTag<TAG>::type InnerTag;
3763  static const std::string n = std::string("RootDivideUnbiased<") + InnerTag::name() + " >";
3764  return n;
3765  }
3766 
3767  template <class U, class BASE>
3768  struct Impl
3769  : public BASE
3770  {
3771  typedef typename LookupDependency<TargetTag, BASE>::value_type value_type;
3772  typedef value_type result_type;
3773 
3774  result_type operator()() const
3775  {
3776  using namespace multi_math;
3777  return sqrt(getDependency<TargetTag>(*this));
3778  }
3779  };
3780 };
3781 
3782 /** \brief Spezialization: works in pass 1, %operator+=() supported (merging supported).
3783 */
3784 template <>
3785 class Central<PowerSum<2> >
3786 {
3787  public:
3789 
3790  static std::string const & name()
3791  {
3792  static const std::string n("Central<PowerSum<2> >");
3793  return n;
3794  }
3795 
3796  template <class U, class BASE>
3797  struct Impl
3798  : public SumBaseImpl<BASE, U>
3799  {
3800  void operator+=(Impl const & o)
3801  {
3802  using namespace vigra::multi_math;
3803  double n1 = getDependency<Count>(*this), n2 = getDependency<Count>(o);
3804  if(n1 == 0.0)
3805  {
3806  this->value_ = o.value_;
3807  }
3808  else if(n2 != 0.0)
3809  {
3810  this->value_ += o.value_ + n1 * n2 / (n1 + n2) * sq(getDependency<Mean>(*this) - getDependency<Mean>(o));
3811  }
3812  }
3813 
3814  void update(U const & t)
3815  {
3816  double n = getDependency<Count>(*this);
3817  if(n > 1.0)
3818  {
3819  using namespace vigra::multi_math;
3820  this->value_ += n / (n - 1.0) * sq(getDependency<Mean>(*this) - t);
3821  }
3822  }
3823 
3824  void update(U const & t, double weight)
3825  {
3826  double n = getDependency<Count>(*this);
3827  if(n > weight)
3828  {
3829  using namespace vigra::multi_math;
3830  this->value_ += n / (n - weight) * sq(getDependency<Mean>(*this) - t);
3831  }
3832  }
3833  };
3834 };
3835 
3836 /** \brief Specialization: works in pass 2, %operator+=() supported (merging supported).
3837 */
3838 template <>
3839 class Central<PowerSum<3> >
3840 {
3841  public:
3843 
3844  static std::string const & name()
3845  {
3846  static const std::string n("Central<PowerSum<3> >");
3847  return n;
3848  }
3849 
3850  template <class U, class BASE>
3851  struct Impl
3852  : public SumBaseImpl<BASE, U>
3853  {
3854  typedef typename SumBaseImpl<BASE, U>::value_type value_type;
3855 
3856  static const unsigned int workInPass = 2;
3857 
3858  void operator+=(Impl const & o)
3859  {
3860  typedef Central<PowerSum<2> > Sum2Tag;
3861 
3862  using namespace vigra::multi_math;
3863  double n1 = getDependency<Count>(*this), n2 = getDependency<Count>(o);
3864  if(n1 == 0.0)
3865  {
3866  this->value_ = o.value_;
3867  }
3868  else if(n2 != 0.0)
3869  {
3870  double n = n1 + n2;
3871  double weight = n1 * n2 * (n1 - n2) / sq(n);
3872  value_type delta = getDependency<Mean>(o) - getDependency<Mean>(*this);
3873  this->value_ += o.value_ + weight * pow(delta, 3) +
3874  3.0 / n * delta * (n1 * getDependency<Sum2Tag>(o) - n2 * getDependency<Sum2Tag>(*this));
3875  }
3876  }
3877 
3878  void update(U const & t)
3879  {
3880  using namespace vigra::multi_math;
3881  this->value_ += pow(getDependency<Centralize>(*this), 3);
3882  }
3883 
3884  void update(U const & t, double weight)
3885  {
3886  using namespace vigra::multi_math;
3887  this->value_ += weight*pow(getDependency<Centralize>(*this), 3);
3888  }
3889  };
3890 };
3891 /** \brief Specialization: works in pass 2, %operator+=() supported (merging supported).
3892 */
3893 template <>
3894 class Central<PowerSum<4> >
3895 {
3896  public:
3898 
3899  static std::string const & name()
3900  {
3901  static const std::string n("Central<PowerSum<4> >");
3902  return n;
3903  }
3904 
3905  template <class U, class BASE>
3906  struct Impl
3907  : public SumBaseImpl<BASE, U>
3908  {
3909  typedef typename SumBaseImpl<BASE, U>::value_type value_type;
3910 
3911  static const unsigned int workInPass = 2;
3912 
3913  void operator+=(Impl const & o)
3914  {
3915  typedef Central<PowerSum<2> > Sum2Tag;
3916  typedef Central<PowerSum<3> > Sum3Tag;
3917 
3918  using namespace vigra::multi_math;
3919  double n1 = getDependency<Count>(*this), n2 = getDependency<Count>(o);
3920  if(n1 == 0.0)
3921  {
3922  this->value_ = o.value_;
3923  }
3924  else if(n2 != 0.0)
3925  {
3926  double n = n1 + n2;
3927  double n1_2 = sq(n1);
3928  double n2_2 = sq(n2);
3929  double n_2 = sq(n);
3930  double weight = n1 * n2 * (n1_2 - n1*n2 + n2_2) / n_2 / n;
3931  value_type delta = getDependency<Mean>(o) - getDependency<Mean>(*this);
3932  this->value_ += o.value_ + weight * pow(delta, 4) +
3933  6.0 / n_2 * sq(delta) * (n1_2 * getDependency<Sum2Tag>(o) + n2_2 * getDependency<Sum2Tag>(*this)) +
3934  4.0 / n * delta * (n1 * getDependency<Sum3Tag>(o) - n2 * getDependency<Sum3Tag>(*this));
3935  }
3936  }
3937 
3938  void update(U const & t)
3939  {
3940  using namespace vigra::multi_math;
3941  this->value_ += pow(getDependency<Centralize>(*this), 4);
3942  }
3943 
3944  void update(U const & t, double weight)
3945  {
3946  using namespace vigra::multi_math;
3947  this->value_ += weight*pow(getDependency<Centralize>(*this), 4);
3948  }
3949  };
3950 };
3951 
3952 /** \brief Basic statistic. Skewness.
3953 
3954  %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$ .
3955  Works in pass 2, %operator+=() supported (merging supported).
3956 */
3958 {
3959  public:
3961 
3962  static std::string const & name()
3963  {
3964  static const std::string n("Skewness");
3965  return n;
3966  }
3967 
3968  template <class U, class BASE>
3969  struct Impl
3970  : public BASE
3971  {
3972  static const unsigned int workInPass = 2;
3973 
3974  typedef typename LookupDependency<Central<PowerSum<3> >, BASE>::value_type value_type;
3975  typedef value_type result_type;
3976 
3977  result_type operator()() const
3978  {
3979  typedef Central<PowerSum<3> > Sum3;
3980  typedef Central<PowerSum<2> > Sum2;
3981 
3982  using namespace multi_math;
3983  return sqrt(getDependency<Count>(*this)) * getDependency<Sum3>(*this) / pow(getDependency<Sum2>(*this), 1.5);
3984  }
3985  };
3986 };
3987 
3988 /** \brief Basic statistic. Unbiased Skewness.
3989 
3990  Works in pass 2, %operator+=() supported (merging supported).
3991 */
3993 {
3994  public:
3996 
3997  static std::string const & name()
3998  {
3999  static const std::string n("UnbiasedSkewness");
4000  return n;
4001  }
4002 
4003  template <class U, class BASE>
4004  struct Impl
4005  : public BASE
4006  {
4007  static const unsigned int workInPass = 2;
4008 
4009  typedef typename LookupDependency<Central<PowerSum<3> >, BASE>::value_type value_type;
4010  typedef value_type result_type;
4011 
4012  result_type operator()() const
4013  {
4014  using namespace multi_math;
4015  double n = getDependency<Count>(*this);
4016  return sqrt(n*(n-1.0)) / (n - 2.0) * getDependency<Skewness>(*this);
4017  }
4018  };
4019 };
4020 
4021 /** \brief Basic statistic. Kurtosis.
4022 
4023  %Kurtosis = @f$ \frac{ \frac{1}{n}\sum_i (x_i-\bar{x})^4 }{
4024  (\frac{1}{n} \sum_i(x_i-\bar{x})^2)^2 } - 3 @f$ .
4025  Works in pass 2, %operator+=() supported (merging supported).
4026 */
4028 {
4029  public:
4031 
4032  static std::string const & name()
4033  {
4034  static const std::string n("Kurtosis");
4035  return n;
4036  }
4037 
4038  template <class U, class BASE>
4039  struct Impl
4040  : public BASE
4041  {
4042  static const unsigned int workInPass = 2;
4043 
4044  typedef typename LookupDependency<Central<PowerSum<4> >, BASE>::value_type value_type;
4045  typedef value_type result_type;
4046 
4047  result_type operator()() const
4048  {
4049  typedef Central<PowerSum<4> > Sum4;
4050  typedef Central<PowerSum<2> > Sum2;
4051 
4052  using namespace multi_math;
4053  return getDependency<Count>(*this) * getDependency<Sum4>(*this) / sq(getDependency<Sum2>(*this)) - value_type(3.0);
4054  }
4055  };
4056 };
4057 
4058 /** \brief Basic statistic. Unbiased Kurtosis.
4059 
4060  Works in pass 2, %operator+=() supported (merging supported).
4061 */
4063 {
4064  public:
4066 
4067  static std::string const & name()
4068  {
4069  static const std::string n("UnbiasedKurtosis");
4070  return n;
4071  }
4072 
4073  template <class U, class BASE>
4074  struct Impl
4075  : public BASE
4076  {
4077  static const unsigned int workInPass = 2;
4078 
4079  typedef typename LookupDependency<Central<PowerSum<4> >, BASE>::value_type value_type;
4080  typedef value_type result_type;
4081 
4082  result_type operator()() const
4083  {
4084  using namespace multi_math;
4085  double n = getDependency<Count>(*this);
4086  return (n-1.0)/((n-2.0)*(n-3.0))*((n+1.0)*getDependency<Kurtosis>(*this) + value_type(6.0));
4087  }
4088  };
4089 };
4090 
4091 namespace detail {
4092 
4093 template <class Scatter, class Sum>
4094 void updateFlatScatterMatrix(Scatter & sc, Sum const & s, double w)
4095 {
4096  int size = s.size();
4097  for(MultiArrayIndex j=0, k=0; j<size; ++j)
4098  for(MultiArrayIndex i=j; i<size; ++i, ++k)
4099  sc[k] += w*s[i]*s[j];
4100 }
4101 
4102 template <class Sum>
4103 void updateFlatScatterMatrix(double & sc, Sum const & s, double w)
4104 {
4105  sc += w*s*s;
4106 }
4107 
4108 template <class Cov, class Scatter>
4109 void flatScatterMatrixToScatterMatrix(Cov & cov, Scatter const & sc)
4110 {
4111  int size = cov.shape(0), k=0;
4112  for(MultiArrayIndex j=0; j<size; ++j)
4113  {
4114  cov(j,j) = sc[k++];
4115  for(MultiArrayIndex i=j+1; i<size; ++i)
4116  {
4117  cov(i,j) = sc[k++];
4118  cov(j,i) = cov(i,j);
4119  }
4120  }
4121 }
4122 
4123 template <class Scatter>
4124 void flatScatterMatrixToScatterMatrix(double & cov, Scatter const & sc)
4125 {
4126  cov = sc;
4127 }
4128 
4129 template <class Cov, class Scatter>
4130 void flatScatterMatrixToCovariance(Cov & cov, Scatter const & sc, double n)
4131 {
4132  int size = cov.shape(0), k=0;
4133  for(MultiArrayIndex j=0; j<size; ++j)
4134  {
4135  cov(j,j) = sc[k++] / n;
4136  for(MultiArrayIndex i=j+1; i<size; ++i)
4137  {
4138  cov(i,j) = sc[k++] / n;
4139  cov(j,i) = cov(i,j);
4140  }
4141  }
4142 }
4143 
4144 template <class Scatter>
4145 void flatScatterMatrixToCovariance(double & cov, Scatter const & sc, double n)
4146 {
4147  cov = sc / n;
4148 }
4149 
4150 } // namespace detail
4151 
4152 // we only store the flattened upper triangular part of the scatter matrix
4153 /** \brief Basic statistic. Flattened uppter-triangular part of scatter matrix.
4154 
4155  Works in pass 1, %operator+=() supported (merging supported).
4156 */
4158 {
4159  public:
4161 
4162  static std::string const & name()
4163  {
4164  static const std::string n("FlatScatterMatrix");
4165  return n;
4166  }
4167 
4168  template <class U, class BASE>
4169  struct Impl
4170  : public BASE
4171  {
4172  typedef typename AccumulatorResultTraits<U>::element_promote_type element_type;
4173  typedef typename AccumulatorResultTraits<U>::FlatCovarianceType value_type;
4174  typedef value_type const & result_type;
4175 
4176  typedef typename AccumulatorResultTraits<U>::SumType SumType;
4177 
4178  value_type value_;
4179  SumType diff_;
4180 
4181  Impl()
4182  : value_(), // call default constructor explicitly to ensure zero initialization
4183  diff_()
4184  {}
4185 
4186  void reset()
4187  {
4188  value_ = element_type();
4189  }
4190 
4191  template <class Shape>
4192  void reshape(Shape const & s)
4193  {
4194  int size = prod(s);
4195  detail::reshapeImpl(value_, Shape1(size*(size+1)/2));
4196  detail::reshapeImpl(diff_, s);
4197  }
4198 
4199  void operator+=(Impl const & o)
4200  {
4201  double n1 = getDependency<Count>(*this), n2 = getDependency<Count>(o);
4202  if(n1 == 0.0)
4203  {
4204  value_ = o.value_;
4205  }
4206  else if(n2 != 0.0)
4207  {
4208  using namespace vigra::multi_math;
4209  diff_ = getDependency<Mean>(*this) - getDependency<Mean>(o);
4210  detail::updateFlatScatterMatrix(value_, diff_, n1 * n2 / (n1 + n2));
4211  value_ += o.value_;
4212  }
4213  }
4214 
4215  void update(U const & t)
4216  {
4217  compute(t);
4218  }
4219 
4220  void update(U const & t, double weight)
4221  {
4222  compute(t, weight);
4223  }
4224 
4225  result_type operator()() const
4226  {
4227  return value_;
4228  }
4229 
4230  private:
4231  void compute(U const & t, double weight = 1.0)
4232  {
4233  double n = getDependency<Count>(*this);
4234  if(n > weight)
4235  {
4236  using namespace vigra::multi_math;
4237  diff_ = getDependency<Mean>(*this) - t;
4238  detail::updateFlatScatterMatrix(value_, diff_, n * weight / (n - weight));
4239  }
4240  }
4241  };
4242 };
4243 
4244 // Covariance
4245 template <>
4247 {
4248  public:
4249  typedef Select<FlatScatterMatrix, Count> Dependencies;
4250 
4251  static std::string const & name()
4252  {
4253  static const std::string n("DivideByCount<FlatScatterMatrix>");
4254  return n;
4255  }
4256 
4257  template <class U, class BASE>
4258  struct Impl
4259  : public CachedResultBase<BASE, typename AccumulatorResultTraits<U>::CovarianceType, U>
4260  {
4261  typedef CachedResultBase<BASE, typename AccumulatorResultTraits<U>::CovarianceType, U> BaseType;
4262  typedef typename BaseType::result_type result_type;
4263 
4264  template <class Shape>
4265  void reshape(Shape const & s)
4266  {
4267  int size = prod(s);
4268  detail::reshapeImpl(this->value_, Shape2(size,size));
4269  }
4270 
4271  result_type operator()() const
4272  {
4273  if(this->isDirty())
4274  {
4275  detail::flatScatterMatrixToCovariance(this->value_, getDependency<FlatScatterMatrix>(*this), getDependency<Count>(*this));
4276  this->setClean();
4277  }
4278  return this->value_;
4279  }
4280  };
4281 };
4282 
4283 // UnbiasedCovariance
4284 template <>
4285 class DivideUnbiased<FlatScatterMatrix>
4286 {
4287  public:
4288  typedef Select<FlatScatterMatrix, Count> Dependencies;
4289 
4290  static std::string const & name()
4291  {
4292  static const std::string n("DivideUnbiased<FlatScatterMatrix>");
4293  return n;
4294  }
4295 
4296  template <class U, class BASE>
4297  struct Impl
4298  : public CachedResultBase<BASE, typename AccumulatorResultTraits<U>::CovarianceType, U>
4299  {
4300  typedef CachedResultBase<BASE, typename AccumulatorResultTraits<U>::CovarianceType, U> BaseType;
4301  typedef typename BaseType::result_type result_type;
4302 
4303  template <class Shape>
4304  void reshape(Shape const & s)
4305  {
4306  int size = prod(s);
4307  detail::reshapeImpl(this->value_, Shape2(size,size));
4308  }
4309 
4310  result_type operator()() const
4311  {
4312  if(this->isDirty())
4313  {
4314  detail::flatScatterMatrixToCovariance(this->value_, getDependency<FlatScatterMatrix>(*this), getDependency<Count>(*this) - 1.0);
4315  this->setClean();
4316  }
4317  return this->value_;
4318  }
4319  };
4320 };
4321 
4322 /** Basic statistic. ScatterMatrixEigensystem (?)
4323 */
4325 {
4326  public:
4328 
4329  static std::string const & name()
4330  {
4331  static const std::string n("ScatterMatrixEigensystem");
4332  return n;
4333  }
4334 
4335  template <class U, class BASE>
4336  struct Impl
4337  : public BASE
4338  {
4339  typedef typename AccumulatorResultTraits<U>::element_promote_type element_type;
4340  typedef typename AccumulatorResultTraits<U>::SumType EigenvalueType;
4341  typedef typename AccumulatorResultTraits<U>::CovarianceType EigenvectorType;
4342  typedef std::pair<EigenvalueType, EigenvectorType> value_type;
4343  typedef value_type const & result_type;
4344 
4345  mutable value_type value_;
4346 
4347  Impl()
4348  : value_()
4349  {}
4350 
4351  void operator+=(Impl const &)
4352  {
4353  this->setDirty();
4354  }
4355 
4356  void update(U const &)
4357  {
4358  this->setDirty();
4359  }
4360 
4361  void update(U const &, double)
4362  {
4363  this->setDirty();
4364  }
4365 
4366  void reset()
4367  {
4368  value_.first = element_type();
4369  value_.second = element_type();
4370  this->setClean();
4371  }
4372 
4373  template <class Shape>
4374  void reshape(Shape const & s)
4375  {
4376  int size = prod(s);
4377  detail::reshapeImpl(value_.first, Shape1(size));
4378  detail::reshapeImpl(value_.second, Shape2(size,size));
4379  }
4380 
4381  result_type operator()() const
4382  {
4383  if(this->isDirty())
4384  {
4385  compute(getDependency<FlatScatterMatrix>(*this), value_.first, value_.second);
4386  this->setClean();
4387  }
4388  return value_;
4389  }
4390 
4391  private:
4392  template <class Flat, class EW, class EV>
4393  static void compute(Flat const & flatScatter, EW & ew, EV & ev)
4394  {
4395  EigenvectorType scatter(ev.shape());
4396  detail::flatScatterMatrixToScatterMatrix(scatter, flatScatter);
4397  // create a view because EW could be a TinyVector
4398  MultiArrayView<2, element_type> ewview(Shape2(ev.shape(0), 1), &ew[0]);
4399  symmetricEigensystem(scatter, ewview, ev);
4400  }
4401 
4402  static void compute(double v, double & ew, double & ev)
4403  {
4404  ew = v;
4405  ev = 1.0;
4406  }
4407  };
4408 };
4409 
4410 // CovarianceEigensystem
4411 template <>
4413 {
4414  public:
4415  typedef Select<ScatterMatrixEigensystem, Count> Dependencies;
4416 
4417  static std::string const & name()
4418  {
4419  static const std::string n("DivideByCount<ScatterMatrixEigensystem>");
4420  return n;
4421  }
4422 
4423  template <class U, class BASE>
4424  struct Impl
4425  : public BASE
4426  {
4427  typedef typename LookupDependency<ScatterMatrixEigensystem, BASE>::type SMImpl;
4428  typedef typename SMImpl::element_type element_type;
4429  typedef typename SMImpl::EigenvalueType EigenvalueType;
4430  typedef typename SMImpl::EigenvectorType EigenvectorType;
4431  typedef std::pair<EigenvalueType, EigenvectorType const &> value_type;
4432  typedef value_type const & result_type;
4433 
4434  mutable value_type value_;
4435 
4436  Impl()
4437  : value_(EigenvalueType(), BASE::template call_getDependency<ScatterMatrixEigensystem>().second)
4438  {}
4439 
4440  void operator+=(Impl const &)
4441  {
4442  this->setDirty();
4443  }
4444 
4445  void update(U const &)
4446  {
4447  this->setDirty();
4448  }
4449 
4450  void update(U const &, double)
4451  {
4452  this->setDirty();
4453  }
4454 
4455  void reset()
4456  {
4457  value_.first = element_type();
4458  this->setClean();
4459  }
4460 
4461  template <class Shape>
4462  void reshape(Shape const & s)
4463  {
4464  int size = prod(s);
4465  detail::reshapeImpl(value_.first, Shape2(size,1));
4466  }
4467 
4468  result_type operator()() const
4469  {
4470  if(this->isDirty())
4471  {
4472  value_.first = getDependency<ScatterMatrixEigensystem>(*this).first / getDependency<Count>(*this);
4473  this->setClean();
4474  }
4475  return value_;
4476  }
4477  };
4478 };
4479 
4480 // alternative implementation of CovarianceEigensystem - solve eigensystem directly
4481 //
4482 // template <>
4483 // class DivideByCount<ScatterMatrixEigensystem>
4484 // {
4485  // public:
4486  // typedef Select<Covariance> Dependencies;
4487 
4488  // template <class U, class BASE>
4489  // struct Impl
4490  // : public BASE
4491  // {
4492  // typedef typename AccumulatorResultTraits<U>::element_promote_type element_type;
4493  // typedef typename AccumulatorResultTraits<U>::SumType EigenvalueType;
4494  // typedef typename AccumulatorResultTraits<U>::CovarianceType EigenvectorType;
4495  // typedef std::pair<EigenvalueType, EigenvectorType> value_type;
4496  // typedef value_type const & result_type;
4497 
4498  // mutable value_type value_;
4499 
4500  // Impl()
4501  // : value_()
4502  // {}
4503 
4504  // void operator+=(Impl const &)
4505  // {
4506  // this->setDirty();
4507  // }
4508 
4509  // void update(U const &)
4510  // {
4511  // this->setDirty();
4512  // }
4513 
4514  // void update(U const &, double)
4515  // {
4516  // this->setDirty();
4517  // }
4518 
4519  // void reset()
4520  // {
4521  // value_.first = element_type();
4522  // value_.second = element_type();
4523  // this->setClean();
4524  // }
4525 
4526  // template <class Shape>
4527  // void reshape(Shape const & s)
4528  // {
4529  // int size = prod(s);
4530  // detail::reshapeImpl(value_.first, Shape2(size,1));
4531  // detail::reshapeImpl(value_.second, Shape2(size,size));
4532  // }
4533 
4534  // result_type operator()() const
4535  // {
4536  // if(this->isDirty())
4537  // {
4538  // compute(getDependency<Covariance>(*this), value_.first, value_.second);
4539  // this->setClean();
4540  // }
4541  // return value_;
4542  // }
4543 
4544  // private:
4545  // template <class Cov, class EW, class EV>
4546  // static void compute(Cov const & cov, EW & ew, EV & ev)
4547  // {
4548  // // create a view because EW could be a TinyVector
4549  // MultiArrayView<2, element_type> ewview(Shape2(cov.shape(0), 1), &ew[0]);
4550  // symmetricEigensystem(cov, ewview, ev);
4551  // }
4552 
4553  // static void compute(double cov, double & ew, double & ev)
4554  // {
4555  // ew = cov;
4556  // ev = 1.0;
4557  // }
4558  // };
4559 // };
4560 
4561 // covariance eigenvalues
4562 /** \brief Specialization (covariance eigenvalues): works in pass 1, %operator+=() supported (merging).
4563 */
4564 template <>
4566 {
4567  public:
4569 
4570  static std::string const & name()
4571  {
4572  static const std::string n("Principal<PowerSum<2> >");
4573  return n;
4574  }
4575 
4576  template <class U, class BASE>
4577  struct Impl
4578  : public BASE
4579  {
4580  typedef typename LookupDependency<ScatterMatrixEigensystem, BASE>::type::EigenvalueType value_type;
4581  typedef value_type const & result_type;
4582 
4583  result_type operator()() const
4584  {
4585  return getDependency<ScatterMatrixEigensystem>(*this).first;
4586  }
4587  };
4588 };
4589 
4590 
4591 // Principal<CoordinateSystem> == covariance eigenvectors
4592 /** \brief Specialization (covariance eigenvectors): works in pass 1, %operator+=() supported (merging).
4593 */
4594 template <>
4596 {
4597  public:
4599 
4600  static std::string const & name()
4601  {
4602  static const std::string n("Principal<CoordinateSystem>");
4603  return n;
4604  }
4605 
4606  template <class U, class BASE>
4607  struct Impl
4608  : public BASE
4609  {
4610  typedef typename LookupDependency<ScatterMatrixEigensystem, BASE>::type::EigenvectorType value_type;
4611  typedef value_type const & result_type;
4612 
4613  result_type operator()() const
4614  {
4615  return getDependency<ScatterMatrixEigensystem>(*this).second;
4616  }
4617  };
4618 };
4619 
4620 /** \brief Basic statistic. %Minimum value.
4621 
4622  Works in pass 1, %operator+=() supported (merging supported).
4623 */
4624 class Minimum
4625 {
4626  public:
4627  typedef Select<> Dependencies;
4628 
4629  static std::string const & name()
4630  {
4631  static const std::string n("Minimum");
4632  return n;
4633  }
4634 
4635  template <class U, class BASE>
4636  struct Impl
4637  : public BASE
4638  {
4639  typedef typename AccumulatorResultTraits<U>::element_type element_type;
4640  typedef typename AccumulatorResultTraits<U>::MinmaxType value_type;
4641  typedef value_type const & result_type;
4642 
4643  value_type value_;
4644 
4645  Impl()
4646  {
4647  value_ = NumericTraits<element_type>::max();
4648  }
4649 
4650  void reset()
4651  {
4652  value_ = NumericTraits<element_type>::max();
4653  }
4654 
4655  template <class Shape>
4656  void reshape(Shape const & s)
4657  {
4658  detail::reshapeImpl(value_, s, NumericTraits<element_type>::max());
4659  }
4660 
4661  void operator+=(Impl const & o)
4662  {
4663  updateImpl(o.value_); // necessary because std::min causes ambiguous overload
4664  }
4665 
4666  void update(U const & t)
4667  {
4668  updateImpl(t);
4669  }
4670 
4671  void update(U const & t, double)
4672  {
4673  updateImpl(t);
4674  }
4675 
4676  result_type operator()() const
4677  {
4678  return value_;
4679  }
4680 
4681  private:
4682  template <class T>
4683  void updateImpl(T const & o)
4684  {
4685  using namespace multi_math;
4686  value_ = min(value_, o);
4687  }
4688 
4689  template <class T, class Alloc>
4690  void updateImpl(MultiArray<1, T, Alloc> const & o)
4691  {
4692  value_ = multi_math::min(value_, o);
4693  }
4694  };
4695 };
4696 
4697 /** \brief Basic statistic. %Maximum value.
4698 
4699  Works in pass 1, %operator+=() supported (merging supported).
4700 */
4701 class Maximum
4702 {
4703  public:
4704  typedef Select<> Dependencies;
4705 
4706  static std::string const & name()
4707  {
4708  static const std::string n("Maximum");
4709  return n;
4710  }
4711 
4712  template <class U, class BASE>
4713  struct Impl
4714  : public BASE
4715  {
4716  typedef typename AccumulatorResultTraits<U>::element_type element_type;
4717  typedef typename AccumulatorResultTraits<U>::MinmaxType value_type;
4718  typedef value_type const & result_type;
4719 
4720  value_type value_;
4721 
4722  Impl()
4723  {
4724  value_ = NumericTraits<element_type>::min();
4725  }
4726 
4727  void reset()
4728  {
4729  value_ = NumericTraits<element_type>::min();
4730  }
4731 
4732  template <class Shape>
4733  void reshape(Shape const & s)
4734  {
4735  detail::reshapeImpl(value_, s, NumericTraits<element_type>::min());
4736  }
4737 
4738  void operator+=(Impl const & o)
4739  {
4740  updateImpl(o.value_); // necessary because std::max causes ambiguous overload
4741  }
4742 
4743  void update(U const & t)
4744  {
4745  updateImpl(t);
4746  }
4747 
4748  void update(U const & t, double)
4749  {
4750  updateImpl(t);
4751  }
4752 
4753  result_type operator()() const
4754  {
4755  return value_;
4756  }
4757 
4758  private:
4759  template <class T>
4760  void updateImpl(T const & o)
4761  {
4762  using namespace multi_math;
4763  value_ = max(value_, o);
4764  }
4765 
4766  template <class T, class Alloc>
4767  void updateImpl(MultiArray<1, T, Alloc> const & o)
4768  {
4769  value_ = multi_math::max(value_, o);
4770  }
4771  };
4772 };
4773 
4774 /** \brief Basic statistic. Data value where weight assumes its minimal value.
4775 
4776  Weights must be given. Coord<ArgMinWeight> gives coordinate where weight assumes its minimal value. Works in pass 1, %operator+=() supported (merging supported).
4777 */
4779 {
4780  public:
4781  typedef Select<> Dependencies;
4782 
4783  static std::string const & name()
4784  {
4785  static const std::string n("ArgMinWeight");
4786  return n;
4787  }
4788 
4789  template <class U, class BASE>
4790  struct Impl
4791  : public BASE
4792  {
4793  typedef typename AccumulatorResultTraits<U>::element_type element_type;
4794  typedef typename AccumulatorResultTraits<U>::MinmaxType value_type;
4795  typedef value_type const & result_type;
4796 
4797  double min_weight_;
4798  value_type value_;
4799 
4800  Impl()
4801  : min_weight_(NumericTraits<double>::max()),
4802  value_()
4803  {}
4804 
4805  void reset()
4806  {
4807  min_weight_ = NumericTraits<double>::max();
4808  value_ = element_type();
4809  }
4810 
4811  template <class Shape>
4812  void reshape(Shape const & s)
4813  {
4814  detail::reshapeImpl(value_, s);
4815  }
4816 
4817  void operator+=(Impl const & o)
4818  {
4819  using namespace multi_math;
4820  if(o.min_weight_ < min_weight_)
4821  {
4822  min_weight_ = o.min_weight_;
4823  value_ = o.value_;
4824  }
4825  }
4826 
4827  void update(U const & t)
4828  {
4829  vigra_precondition(false, "ArgMinWeight::update() needs weights.");
4830  }
4831 
4832  void update(U const & t, double weight)
4833  {
4834  if(weight < min_weight_)
4835  {
4836  min_weight_ = weight;
4837  value_ = t;
4838  }
4839  }
4840 
4841  result_type operator()() const
4842  {
4843  return value_;
4844  }
4845  };
4846 };
4847 
4848 /** \brief Basic statistic. Data where weight assumes its maximal value.
4849 
4850  Weights must be given. Coord<ArgMinWeight> gives coordinate where weight assumes its maximal value. Works in pass 1, %operator+=() supported (merging supported).
4851 */
4853 {
4854  public:
4855  typedef Select<> Dependencies;
4856 
4857  static std::string const & name()
4858  {
4859  static const std::string n("ArgMaxWeight");
4860  return n;
4861  }
4862 
4863  template <class U, class BASE>
4864  struct Impl
4865  : public BASE
4866  {
4867  typedef typename AccumulatorResultTraits<U>::element_type element_type;
4868  typedef typename AccumulatorResultTraits<U>::MinmaxType value_type;
4869  typedef value_type const & result_type;
4870 
4871  double max_weight_;
4872  value_type value_;
4873 
4874  Impl()
4875  : max_weight_(NumericTraits<double>::min()),
4876  value_()
4877  {}
4878 
4879  void reset()
4880  {
4881  max_weight_ = NumericTraits<double>::min();
4882  value_ = element_type();
4883  }
4884 
4885  template <class Shape>
4886  void reshape(Shape const & s)
4887  {
4888  detail::reshapeImpl(value_, s);
4889  }
4890 
4891  void operator+=(Impl const & o)
4892  {
4893  using namespace multi_math;
4894  if(o.max_weight_ > max_weight_)
4895  {
4896  max_weight_ = o.max_weight_;
4897  value_ = o.value_;
4898  }
4899  }
4900 
4901  void update(U const & t)
4902  {
4903  vigra_precondition(false, "ArgMaxWeight::update() needs weights.");
4904  }
4905 
4906  void update(U const & t, double weight)
4907  {
4908  if(weight > max_weight_)
4909  {
4910  max_weight_ = weight;
4911  value_ = t;
4912  }
4913  }
4914 
4915  result_type operator()() const
4916  {
4917  return value_;
4918  }
4919  };
4920 };
4921 
4922 
4923 template <class BASE, int BinCount>
4924 class HistogramBase
4925 : public BASE
4926 {
4927  public:
4928 
4929  typedef double element_type;
4930  typedef TinyVector<double, BinCount> value_type;
4931  typedef value_type const & result_type;
4932 
4933  value_type value_;
4934  double left_outliers, right_outliers;
4935 
4936  HistogramBase()
4937  : value_(),
4938  left_outliers(),
4939  right_outliers()
4940  {}
4941 
4942  void reset()
4943  {
4944  value_ = element_type();
4945  left_outliers = 0.0;
4946  right_outliers = 0.0;
4947  }
4948 
4949  void operator+=(HistogramBase const & o)
4950  {
4951  value_ += o.value_;
4952  left_outliers += o.left_outliers;
4953  right_outliers += o.right_outliers;
4954  }
4955 
4956  result_type operator()() const
4957  {
4958  return value_;
4959  }
4960 };
4961 
4962 template <class BASE>
4963 class HistogramBase<BASE, 0>
4964 : public BASE
4965 {
4966  public:
4967 
4968  typedef double element_type;
4969  typedef MultiArray<1, double> value_type;
4970  typedef value_type const & result_type;
4971 
4972  value_type value_;
4973  double left_outliers, right_outliers;
4974 
4975  HistogramBase()
4976  : value_(),
4977  left_outliers(),
4978  right_outliers()
4979  {}
4980 
4981  void reset()
4982  {
4983  value_ = element_type();
4984  left_outliers = 0.0;
4985  right_outliers = 0.0;
4986  }
4987 
4988  void operator+=(HistogramBase const & o)
4989  {
4990  value_ += o.value_;
4991  left_outliers += o.left_outliers;
4992  right_outliers += o.right_outliers;
4993  }
4994 
4995  void setBinCount(int binCount)
4996  {
4997  vigra_precondition(binCount > 0,
4998  "HistogramBase:.setBinCount(): binCount > 0 required.");
4999  value_type(Shape1(binCount)).swap(value_);
5000  }
5001 
5002  result_type operator()() const
5003  {
5004  return value_;
5005  }
5006 };
5007 
5008 template <class BASE, int BinCount, class U=typename BASE::input_type>
5009 class RangeHistogramBase
5010 : public HistogramBase<BASE, BinCount>
5011 {
5012  public:
5013  double scale_, offset_, inverse_scale_;
5014 
5015  RangeHistogramBase()
5016  : scale_(),
5017  offset_(),
5018  inverse_scale_()
5019  {}
5020 
5021  void reset()
5022  {
5023  scale_ = 0.0;
5024  offset_ = 0.0;
5025  inverse_scale_ = 0.0;
5026  HistogramBase<BASE, BinCount>::reset();
5027  }
5028 
5029  void operator+=(RangeHistogramBase const & o)
5030  {
5031  vigra_precondition(scale_ == 0.0 || o.scale_ == 0.0 || (scale_ == o.scale_ && offset_ == o.offset_),
5032  "RangeHistogramBase::operator+=(): cannot merge histograms with different data mapping.");
5033 
5035  if(scale_ == 0.0)
5036  {
5037  scale_ = o.scale_;
5038  offset_ = o.offset_;
5039  inverse_scale_ = o.inverse_scale_;
5040  }
5041  }
5042 
5043  void update(U const & t)
5044  {
5045  update(t, 1.0);
5046  }
5047 
5048  void update(U const & t, double weight)
5049  {
5050  double m = mapItem(t);
5051  int index = (m == (double)this->value_.size())
5052  ? (int)m - 1
5053  : (int)m;
5054  if(index < 0)
5055  this->left_outliers += weight;
5056  else if(index >= (int)this->value_.size())
5057  this->right_outliers += weight;
5058  else
5059  this->value_[index] += weight;
5060  }
5061 
5062  void setMinMax(double mi, double ma)
5063  {
5064  vigra_precondition(this->value_.size() > 0,
5065  "RangeHistogramBase::setMinMax(...): setBinCount(...) has not been called.");
5066  vigra_precondition(mi < ma,
5067  "RangeHistogramBase::setMinMax(...): min < max required.");
5068  offset_ = mi;
5069  scale_ = (double)this->value_.size() / (ma - mi);
5070  inverse_scale_ = 1.0 / scale_;
5071  }
5072 
5073  double mapItem(double t) const
5074  {
5075  return scale_ * (t - offset_);
5076  }
5077 
5078  double mapItemInverse(double t) const
5079  {
5080  return inverse_scale_ * t + offset_;
5081  }
5082 
5083  template <class ArrayLike>
5084  void computeStandardQuantiles(double minimum, double maximum, double count,
5085  ArrayLike const & desiredQuantiles, ArrayLike & res) const
5086  {
5087  ArrayVector<double> keypoints, cumhist;
5088  double mappedMinimum = mapItem(minimum);
5089  double mappedMaximum = mapItem(maximum);
5090 
5091  keypoints.push_back(mappedMinimum);
5092  cumhist.push_back(0.0);
5093 
5094  if(this->left_outliers > 0.0)
5095  {
5096  keypoints.push_back(0.0);
5097  cumhist.push_back(this->left_outliers);
5098  }
5099 
5100  int size = (int)this->value_.size();
5101  double cumulative = this->left_outliers;
5102  for(int k=0; k<size; ++k)
5103  {
5104  if(this->value_[k] > 0.0)
5105  {
5106  if(keypoints.back() <= k)
5107  {
5108  keypoints.push_back(k);
5109  cumhist.push_back(cumulative);
5110  }
5111  cumulative += this->value_[k];
5112  keypoints.push_back(k+1);
5113  cumhist.push_back(cumulative);
5114  }
5115  }
5116 
5117  if(this->right_outliers > 0.0)
5118  {
5119  if(keypoints.back() != size)
5120  {
5121  keypoints.push_back(size);
5122  cumhist.push_back(cumulative);
5123  }
5124  keypoints.push_back(mappedMaximum);
5125  cumhist.push_back(count);
5126  }
5127  else
5128  {
5129  keypoints.back() = mappedMaximum;
5130  cumhist.back() = count;
5131  }
5132 
5133  int quantile = 0, end = (int)desiredQuantiles.size();
5134 
5135  if(desiredQuantiles[0] == 0.0)
5136  {
5137  res[0] = minimum;
5138  ++quantile;
5139  }
5140  if(desiredQuantiles[end-1] == 1.0)
5141  {
5142  res[end-1] = maximum;
5143  --end;
5144  }
5145 
5146  int point = 0;
5147  double qcount = count * desiredQuantiles[quantile];
5148  while(quantile < end)
5149  {
5150  if(cumhist[point] < qcount && cumhist[point+1] >= qcount)
5151  {
5152  double t = (qcount - cumhist[point]) / (cumhist[point+1] - cumhist[point]) * (keypoints[point+1] - keypoints[point]);
5153  res[quantile] = mapItemInverse(t + keypoints[point]);
5154  ++quantile;
5155  qcount = count * desiredQuantiles[quantile];
5156  }
5157  else
5158  {
5159  ++point;
5160  }
5161  }
5162  }
5163 };
5164 
5165 /** \brief Histogram where data values are equal to bin indices.
5166 
5167  - If BinCount != 0, the return type of the accumulator is TinyVector<double, BinCount> .
5168  - 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).
5169  - Outliers can be accessed via getAccumulator<IntegerHistogram<Bincount>>(a).left_outliers and getAccumulator<...>(acc_chain).right_outliers.
5170  - 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().
5171  Works in pass 1, %operator+=() supported (merging supported).
5172 */
5173 template <int BinCount>
5174 class IntegerHistogram
5175 {
5176  public:
5177 
5178  typedef Select<> Dependencies;
5179 
5180  static std::string const & name()
5181  {
5182  static const std::string n = std::string("IntegerHistogram<") + asString(BinCount) + ">";
5183  return n;
5184  }
5185 
5186  template <class U, class BASE>
5187  struct Impl
5188  : public HistogramBase<BASE, BinCount>
5189  {
5190  void update(int index)
5191  {
5192  if(index < 0)
5193  ++this->left_outliers;
5194  else if(index >= (int)this->value_.size())
5195  ++this->right_outliers;
5196  else
5197  ++this->value_[index];
5198  }
5199 
5200  void update(int index, double weight)
5201  {
5202  // cannot compute quantile from weighted integer histograms,
5203  // so force people to use UserRangeHistogram or AutoRangeHistogram
5204  vigra_precondition(false, "IntegerHistogram::update(): weighted histograms not supported, use another histogram type.");
5205  }
5206 
5207  template <class ArrayLike>
5208  void computeStandardQuantiles(double minimum, double maximum, double count,
5209  ArrayLike const & desiredQuantiles, ArrayLike & res) const
5210  {
5211  int quantile = 0, end = (int)desiredQuantiles.size();
5212 
5213  if(desiredQuantiles[0] == 0.0)
5214  {
5215  res[0] = minimum;
5216  ++quantile;
5217  }
5218  if(desiredQuantiles[end-1] == 1.0)
5219  {
5220  res[end-1] = maximum;
5221  --end;
5222  }
5223 
5224  count -= 1.0;
5225  int currentBin = 0, size = (int)this->value_.size();
5226  double cumulative1 = this->left_outliers,
5227  cumulative2 = this->value_[currentBin] + cumulative1;
5228 
5229  // add a to the quantiles to account for the fact that counting
5230  // corresponds to 1-based indexing (one element == index 1)
5231  double qcount = desiredQuantiles[quantile]*count + 1.0;
5232 
5233  while(quantile < end)
5234  {
5235  if(cumulative2 == qcount)
5236  {
5237  res[quantile] = currentBin;
5238  ++quantile;
5239  qcount = desiredQuantiles[quantile]*count + 1.0;
5240  }
5241  else if(cumulative2 > qcount)
5242  {
5243  if(cumulative1 > qcount) // in left_outlier bin
5244  {
5245  res[quantile] = minimum;
5246  }
5247  if(cumulative1 + 1.0 > qcount) // between bins
5248  {
5249  res[quantile] = currentBin - 1 + qcount - std::floor(qcount);
5250  }
5251  else // standard case
5252  {
5253  res[quantile] = currentBin;
5254  }
5255  ++quantile;
5256  qcount = desiredQuantiles[quantile]*count + 1.0;
5257  }
5258  else if(currentBin == size-1) // in right outlier bin
5259  {
5260  res[quantile] = maximum;
5261  ++quantile;
5262  qcount = desiredQuantiles[quantile]*count + 1.0;
5263  }
5264  else
5265  {
5266  ++currentBin;
5267  cumulative1 = cumulative2;
5268  cumulative2 += this->value_[currentBin];
5269  }
5270  }
5271  }
5272  };
5273 };
5274 
5275 /** \brief Histogram where user provides bounds for linear range mapping from values to indices.
5276 
5277  - If BinCount != 0, the return type of the accumulator is TinyVector<double, BinCount> .
5278  - 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).
5279  - Bounds for the mapping (min/max) must be set before seeing data by calling getAccumulator<UserRangeHistogram<BinCount> >.setMinMax(min, max).
5280  - Options can also be passed to the accumulator chain via an instance of HistogramOptions .
5281  - Works in pass 1, %operator+=() is supported (merging) if both histograms have the same data mapping.
5282  - Outliers can be accessed via getAccumulator<...>(a).left_outliers and getAccumulator<...>(a).right_outliers.
5283  - 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().
5284 */
5285 template <int BinCount>
5286 class UserRangeHistogram
5287 {
5288  public:
5289 
5290  typedef Select<> Dependencies;
5291 
5292  static std::string const & name()
5293  {
5294  static const std::string n = std::string("UserRangeHistogram<") + asString(BinCount) + ">";
5295  return n;
5296  }
5297 
5298  template <class U, class BASE>
5299  struct Impl
5300  : public RangeHistogramBase<BASE, BinCount, U>
5301  {
5302  void update(U const & t)
5303  {
5304  update(t, 1.0);
5305  }
5306 
5307  void update(U const & t, double weight)
5308  {
5309  vigra_precondition(this->scale_ != 0.0,
5310  "UserRangeHistogram::update(): setMinMax(...) has not been called.");
5311 
5312  RangeHistogramBase<BASE, BinCount, U>::update(t, weight);
5313  }
5314  };
5315 };
5316 
5317 /** \brief Histogram where range mapping bounds are defined by minimum and maximum of data.
5318 
5319  - If BinCount != 0, the return type of the accumulator is TinyVector<double, BinCount> .
5320  - 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).
5321  - Becomes a UserRangeHistogram if min/max is set.
5322  - Works in pass 2, %operator+=() is supported (merging) if both histograms have the same data mapping.
5323  - Outliers can be accessed via getAccumulator<...>(acc_chain).left_outliers and getAccumulator<...>(acc_chain).right_outliers .
5324  - 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().
5325 */
5326 template <int BinCount>
5327 class AutoRangeHistogram
5328 {
5329  public:
5330 
5331  typedef Select<Minimum, Maximum> Dependencies;
5332 
5333  static std::string const & name()
5334  {
5335  static const std::string n = std::string("AutoRangeHistogram<") + asString(BinCount) + ">";
5336  return n;
5337  }
5338 
5339  template <class U, class BASE>
5340  struct Impl
5341  : public RangeHistogramBase<BASE, BinCount, U>
5342  {
5343  static const unsigned int workInPass = LookupDependency<Minimum, BASE>::type::workInPass + 1;
5344 
5345  void update(U const & t)
5346  {
5347  update(t, 1.0);
5348  }
5349 
5350  void update(U const & t, double weight)
5351  {
5352  if(this->scale_ == 0.0)
5353  this->setMinMax(getDependency<Minimum>(*this), getDependency<Maximum>(*this));
5354 
5355  RangeHistogramBase<BASE, BinCount, U>::update(t, weight);
5356  }
5357  };
5358 };
5359 
5360 /** \brief Like AutoRangeHistogram, but use global min/max rather than region min/max.
5361 
5362  - If BinCount != 0, the return type of the accumulator is TinyVector<double, BinCount> .
5363  - 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).
5364  - Becomes a UserRangeHistogram if min/max is set.
5365  - Works in pass 2, %operator+=() is supported (merging) if both histograms have the same data mapping.
5366  - Outliers can be accessed via getAccumulator<GlobalRangeHistogram<Bincount>>(acc_chain).left_outliers and getAccumulator<...>(acc_chain).right_outliers .
5367  - 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().
5368 */
5369 template <int BinCount>
5370 class GlobalRangeHistogram
5371 {
5372  public:
5373 
5374  typedef Select<Global<Minimum>, Global<Maximum>, Minimum, Maximum> Dependencies;
5375 
5376  static std::string const & name()
5377  {
5378  static const std::string n = std::string("GlobalRangeHistogram<") + asString(BinCount) + ">";
5379  return n;
5380  }
5381 
5382  template <class U, class BASE>
5383  struct Impl
5384  : public RangeHistogramBase<BASE, BinCount, U>
5385  {
5386  static const unsigned int workInPass = LookupDependency<Minimum, BASE>::type::workInPass + 1;
5387 
5388  bool useLocalMinimax_;
5389 
5390  Impl()
5391  : useLocalMinimax_(false)
5392  {}
5393 
5394  void setRegionAutoInit(bool locally)
5395  {
5396  this->scale_ = 0.0;
5397  useLocalMinimax_ = locally;
5398  }
5399 
5400  void update(U const & t)
5401  {
5402  update(t, 1.0);
5403  }
5404 
5405  void update(U const & t, double weight)
5406  {
5407  if(this->scale_ == 0.0)
5408  {
5409  if(useLocalMinimax_)
5410  this->setMinMax(getDependency<Minimum>(*this), getDependency<Maximum>(*this));
5411  else
5412  this->setMinMax(getDependency<Global<Minimum> >(*this), getDependency<Global<Maximum> >(*this));
5413  }
5414 
5415  RangeHistogramBase<BASE, BinCount, U>::update(t, weight);
5416  }
5417  };
5418 };
5419 
5420 /** \brief Compute (0%, 10%, 25%, 50%, 75%, 90%, 100%) quantiles from given histogram.
5421 
5422  Return type is TinyVector<double, 7> .
5423 */
5424 template <class HistogramAccumulator>
5425 class StandardQuantiles
5426 {
5427  public:
5428 
5429  typedef typename StandardizeTag<HistogramAccumulator>::type HistogramTag;
5430  typedef Select<HistogramTag, Minimum, Maximum, Count> Dependencies;
5431 
5432  static std::string const & name()
5433  {
5434  static const std::string n = std::string("StandardQuantiles<") + HistogramTag::name() + " >";
5435  return n;
5436  }
5437 
5438  template <class U, class BASE>
5439  struct Impl
5440  : public CachedResultBase<BASE, TinyVector<double, 7>, U>
5441  {
5442  typedef typename CachedResultBase<BASE, TinyVector<double, 7>, U>::result_type result_type;
5443  typedef typename CachedResultBase<BASE, TinyVector<double, 7>, U>::value_type value_type;
5444 
5445  static const unsigned int workInPass = LookupDependency<HistogramTag, BASE>::type::workInPass;
5446 
5447  result_type operator()() const
5448  {
5449  if(this->isDirty())
5450  {
5451  static const double desiredQuantiles[] = {0.0, 0.1, 0.25, 0.5, 0.75, 0.9, 1.0 };
5452  getAccumulator<HistogramTag>(*this).computeStandardQuantiles(getDependency<Minimum>(*this), getDependency<Maximum>(*this),
5453  getDependency<Count>(*this), value_type(desiredQuantiles),
5454  this->value_);
5455  this->setClean();
5456  }
5457  return this->value_;
5458  }
5459  };
5460 };
5461 
5462 }} // namespace vigra::acc
5463 
5464 #endif // VIGRA_ACCUMULATOR_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.9.0 (Sat Oct 5 2013)