Mean.h

Go to the documentation of this file.

#ifndef __GfxTL_MEAN_HEADER__
#define __GfxTL_MEAN_HEADER__
#include <GfxTL/Covariance.h>
#include <GfxTL/Jacobi.h>
#include <GfxTL/MathHelper.h>
#include <GfxTL/MatrixXX.h>
#include <GfxTL/VectorXD.h>
#include <GfxTL/WeightFunc.h>
 
namespace GfxTL
{
    template <class PointT, class PointsForwardIt, class WeightsForwardIt>
    void
    Mean(PointsForwardIt begin, PointsForwardIt end, WeightsForwardIt weights, PointT* mean)
    {
        mean->Zero();
        typename PointT::ScalarType totalWeight = 0;
        for (; begin != end; ++begin, ++weights)
        {
            *mean += typename PointT::ScalarType(*weights) * ((const PointT)(*begin));
            totalWeight += *weights;
        }
        if (totalWeight)
        {
            *mean /= totalWeight;
        }
    }
 
    template <class PointsForwardIt, class WeightsForwardIt, class ScalarT>
    void
    Mean(PointsForwardIt begin,
         PointsForwardIt end,
         WeightsForwardIt weights,
         VectorXD<3, ScalarT>* mean)
    {
        mean->Zero();
        ScalarT totalWeight = 0, w;
        for (; begin != end; ++begin, ++weights)
        {
            w = (ScalarT)*weights;
            (*mean)[0] += w * (*begin)[0];
            (*mean)[1] += w * (*begin)[1];
            (*mean)[2] += w * (*begin)[2];
            totalWeight += w;
        }
        if (totalWeight)
        {
            *mean /= totalWeight;
        }
    }
 
    template <class PointT, class PointsForwardIt>
    void
    Mean(PointsForwardIt begin, PointsForwardIt end, PointT* mean)
    {
        Mean(begin, end, UnitWeightIterator(), mean);
    }
 
    // This computes an average of unoriented normals
    template <class NormalsItT, class WeightsItT, class MeanT>
    bool
    MeanOfNormals(NormalsItT begin, NormalsItT end, WeightsItT weights, MeanT* mean)
    {
        typedef typename MeanT::ScalarType ScalarType;
 
        enum
        {
            Dim = MeanT::Dim
        };
 
        GfxTL::MatrixXX<Dim, Dim, ScalarType> cov, eigenVectors;
        GfxTL::VectorXD<Dim, ScalarType> center, eigenValues;
        center.Zero();
        CovarianceMatrix(center, begin, end, weights, &cov);
        if (!Jacobi(cov, &eigenValues, &eigenVectors))
        {
            mean->Zero();
            return false;
        }
        // find the maximal eigenvalue and corresponding vector
        ScalarType maxEigVal = eigenValues[0];
        unsigned int maxEigIdx = 0;
        for (unsigned int i = 1; i < Dim; ++i)
            if (eigenValues[i] > maxEigVal)
            {
                maxEigVal = eigenValues[i];
                maxEigIdx = i;
            }
        *mean = MeanT(eigenVectors[maxEigIdx]);
        return true;
    }
 
    template <class NormalsItT, class MeanT>
    bool
    MeanOfNormals(NormalsItT begin, NormalsItT end, MeanT* mean)
    {
        return MeanOfNormals(begin, end, UnitWeightIterator(), mean);
    }
}; // namespace GfxTL
 
#endif