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