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

vigra/accumulator.hxx VIGRA

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)