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

accumulator.hxx VIGRA

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

© 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.10.0 (Mon Apr 28 2014)