vigra - Feature Accumulators

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

Feature Accumulators |

The namespace `vigra::acc` provides the function 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 tags (see below), and a template meta-program automatically generates an efficient functor that computes exactly those statistics.

The function 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".

#include <vigra/accumulator.hxx>

Basic statistics:

- PowerSum<N> ( )
- AbsPowerSum<N> ( )
- Skewness, UnbiasedSkewness
- Kurtosis, UnbiasedKurtosis
- Minimum, Maximum
- FlatScatterMatrix (flattened upper-triangular part of scatter matrix)
- 4 histogram classes (see below)
- StandardQuantiles (0%, 10%, 25%, 50%, 75%, 90%, 100%)
- ArgMinWeight, ArgMaxWeight (store data or coordinate where weight assumes its minimal or maximal value)
- CoordinateSystem (identity matrix of appropriate size)

Modifiers: (S is the statistc to be modified)

- Normalization

  DivideByCount<S> | S/Count
  RootDivideByCount<S> | sqrt( S/Count )
  DivideUnbiased<S> | S/(Count-1)
  RootDivideUnbiased<S> | sqrt( S/(Count-1) )

- Data preparation:

  Central<S> | substract mean before computing S
  Principal<S> | project onto PCA eigenvectors
  Whitened<S> | scale to unit variance after PCA
  Coord<S> | compute S from pixel coordinates rather than from pixel values
  Weighted<S> | compute weighted version of S
  Global<S> | compute S globally rather than per region (per region is default if labels are given)

Aliases for a couple of important features are implemented (mainly as `typedef FullName Alias`). The alias names are equivalent to full names. Here are some examples for supported alias names (these examples also show how to compose statistics from the fundamental statistics and modifiers):

Alias | Full Name
--- | ---
Count | PowerSum<0>
Sum | PowerSum<1>
SumOfSquares | PowerSum<2>
Mean | DivideByCount<PowerSum<1>>
RootMeanSquares | RootDivideByCount<PowerSum<2>>
Moment<N> | DivideByCount<PowerSum<N>>
Variance | DivideByCount<Central<PowerSum<2>>>
StdDev | RootDivideByCount<Central<PowerSum<2>>>
Covariance | DivideByCount<FlatScatterMatrix>
RegionCenter | Coord<Mean>
CenterOfMass | Weighted<Coord<Mean>>

There are a few rules for composing statistics:

- modifiers can be specified in any order, but are internally transformed to standard order: Global<Weighted<Coord<normalization<data preparation<basic statistic
- only one normalization modifier and one data preparation modifier (Central or Principal or Whitened) is permitted
- Count ignores all modifiers except Global and Weighted
- Sum ignores Central and Principal, because sum would be zero
- ArgMinWeight and ArgMaxWeight are automatically Weighted

Here is an example how to use acc::AccumulatorChain to compute statistics. (To use Weighted<> or Coord<> modifiers, see below):

        #include <vigra/multi_array.hxx>
        #include <vigra/impex.hxx>
        #include <vigra/accumulator.hxx>
        using namespace vigra::acc;
        typedef double DataType;
        int size = 1000;
        vigra::MultiArray<2, DataType> data(vigra::Shape2(size, size));

        AccumulatorChain<DataType,
            Select<Variance, Mean, StdDev, Minimum, Maximum, RootMeanSquares, Skewness, Covariance> >
        a;

        std::cout << "passes required: " << a.passesRequired() << std::endl;
        extractFeatures(data.begin(), data.end(), a);

        std::cout << "Mean: " << get<Mean>(a) << std::endl;
        std::cout << "Variance: " << get<Variance>(a) << std::endl;

The acc::AccumulatorChain object contains the selected statistics and their dependencies. Statistics have to be wrapped with acc::Select. The statistics are computed with the acc::extractFeatures function and the statistics can be accessed with acc::get .

Rules and notes:

- order of statistics in Select<> is arbitrary
- up to 20 statistics in Select<>, but Select<> can be nested
- dependencies are automatically inserted
- duplicates are automatically removed
- extractFeatures() does as many passes through the data as necessary
- each accumulator only sees data in the appropriate pass (its "working pass")

The Accumulators can also be used with vector-valued data (vigra::RGBValue, vigra::TinyVector, vigra::MultiArray or vigra::MultiArrayView):

        typedef vigra::RGBValue<double> DataType;
        AccumulatorChain<DataType, Select<...> > a;
        ...

To compute weighted statistics (Weighted<>) or statistics over coordinates (Coord<>), the accumulator chain can be used with 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.

These index specifiers are: (INDEX is of type int)

- DataArg<INDEX>: CoupledHandle holds data at index 'INDEX' (default INDEX=1)
- LabelArg<INDEX>: CoupledHandle holds labels at index 'INDEX' (default INDEX=2)
- WeightArg<INDEX>: CoupledHandle holds weights at index 'INDEX' (default INDEX=outermost index)

Pixel coordinates are always at index 0.

        using namespace vigra::acc;
        vigra::MultiArray<3, double> data(...), weights(...);
        typedef vigra::CoupledIteratorType<3, double, double>::type Iterator; //type of the CoupledScanOrderIterator
        typedef Iterator::value_type Handle; //type of the corresponding CoupledHandle

        AccumulatorChain<Handle,
            Select<DataArg<1>, WeightArg<2>,       //where to look in the Handle (coordinates are always arg 0)
               Mean, Variance,                     //statistics over values
               Coord<Mean>, Coord<Variance>,       //statistics over coordinates,
               Weighted<Mean>, Weighted<Variance>, //weighted values,
               Weighted<Coord<Mean> > > >          //weighted coordinates.
        a;

        Iterator start = createCoupledIterator(data, weights); //coord->index 0, data->index 1, weights->index 2
        Iterator end = start.getEndIterator();

        extractFeatures(start,end,a);

To compute region statistics, use acc::AccumulatorChainArray :

        using namespace vigra::acc;
        vigra::MultiArray<3, double> data(...);
        vigra::MultiArray<3, int> labels(...);
        typedef vigra::CoupledIteratorType<3, double, int>::type Iterator;
        typedef Iterator::value_type Handle;

        AccumulatorChainArray<Handle,
            Select<DataArg<1>, LabelArg<2>,       //where to look in the Handle (coordinates are always arg 0)
               Mean, Variance,                    //per-region statistics over values
               Coord<Mean>, Coord<Variance>,      //per-region statistics over coordinates
               Global<Mean>, Global<Variance> > > //global statistics
        a;

        Iterator start = createCoupledIterator(data, labels);
        Iterator end = start.getEndIterator();

        a.ignoreLabel(0); //statistics will not be computed for region 0 (e.g. background)

        extractFeatures(start,end,a);

        int regionlabel = ...;
        std::cout << get<Mean>(a, regionlabel) << std::endl; //get Mean of region with label 'regionlabel'

In some application it will be known only at run-time which statistics have to be computed. An Accumulator with run-time activation is provided by the 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:

        using namespace vigra::acc;
        vigra::MultiArray<2, double> data(...);
        DynamicAccumulatorChain<double,
            Select<Mean, Minimum, Maximum, Variance, StdDev> > a; // at compile-time
        activate<Mean>(a);      //at run-time
        a.activate("Minimum");  //same as activate<Minimum>(a) (alias names are not recognized)

        extractFeatures(data.begin(), data.end(), a);
        std::cout << "Mean: " << get<Mean>(a) << std::endl;       //ok
        //std::cout << "Maximum: " << get<Maximum>(a) << std::endl; // run-time error because Maximum not activated

Likewise, for run-time activation of region statistics, use acc::DynamicAccumulatorChainArray.

Accumulator merging (e.g. for parallelization or hierarchical segmentation) is possible for many accumulators:

        using namespace vigra::acc;
        vigra::MultiArray<2, double> data(...);
        AccumulatorChain<double, Select<Mean, Variance, Skewness> > a, a1, a2;

        extractFeatures(data.begin(), data.end(), a); //process entire data set at once
        extractFeatures(data.begin(), data.begin()+data.size()/2, a1); //process first half
        extractFeatures(data.begin()+data.size()/2, data.end(), a2); //process second half
        a1 += a2; // merge: a1 now equals a0 (with numerical tolerances)

Not all statistics can be merged (e.g. Principal 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:

        using namespace vigra::acc;
        vigra::MultiArray<2, double> data(...);
        AccumulatorChain<double, Select<Mean, Variance> > a;

        extractFeatures(data.begin(), data.begin()+data.size()/2, a); // this works because
        extractFeatures(data.begin()+data.size()/2, data.end(), a);   // all statistics only work in pass 1

Four kinds of histograms are currently implemented:

IntegerHistogram | Data values are equal to bin indices
UserRangeHistogram | User provides lower and upper bounds for linear range mapping from values to indices.
AutoRangeHistogram | Range mapping bounds are defiend by minimum and maximum of the data (2 passes needed!)
GlobalRangeHistogram | Likewise, but use global min/max rather than region min/max as AutoRangeHistogram will

- 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).
- If UserRangeHistogram is used, the bounds for the linear range mapping from values to indices must be set before seeing data (see below).
- 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.
- Merging is supported if the range mapping of the histograms is the same.
- Histogram accumulators have two members for outliers (left_outliers, right_outliers).

With the StandardQuantiles class, histogram quantiles (0%, 10%, 25%, 50%, 75%, 90%, 100%) are computed from a given histgram using linear interpolation. The return type is TinyVector<double, 7> .

Usage:

        using namespace vigra::acc;
        typedef double DataType;
        vigra::MultiArray<2, DataType> data(...);

        typedef UserRangeHistogram<40> SomeHistogram;   //binCount set at compile time
        typedef UserRangeHistogram<0> SomeHistogram2; // binCount must be set at run-time
        typedef AutoRangeHistogram<0> SomeHistogram3;
        typedef StandardQuantiles<SomeHistogram3> Quantiles3;

        AccumulatorChain<DataType, Select<SomeHistogram, SomeHistogram2, SomeHistogram3, Quantiles3> > a;

        //set options for all histograms in the accumulator chain:
        vigra::HistogramOptions histogram_opt;
        histogram_opt = histogram_opt.setBinCount(50);
        //histogram_opt = histogram_opt.setMinMax(0.1, 0.9); // this would set min/max for all three histograms, but range bounds
                                                             // shall be set automatically by min/max of data for SomeHistogram3
        a.setHistogramOptions(histogram_opt);

        // set options for a specific histogram in the accumulator chain:
        getAccumulator<SomeHistogram>(a).setMinMax(0.1, 0.9); // number of bins must be set before setting min/max
        getAccumulator<SomeHistogram2>(a).setMinMax(0.0, 1.0);

        extractFeatures(data.begin(), data.end(), a);

        vigra::TinyVector<double, 40> hist = get<SomeHistogram>(a);
        vigra::MultiArray<1, double> hist2 = get<SomeHistogram2>(a);
        vigra::TinyVector<double, 7> quant = get<Quantiles3>(a);
        double right_outliers = getAccumulator<SomeHistogram>(a).right_outliers;

© 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)