vigra - vigra/accumulator.hxx Source File

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

vigra/accumulator.hxx |

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)