MeanShift.h

Go to the documentation of this file.

#ifndef GRT_MEAN_SHIFT_HEADER
#define GRT_MEAN_SHIFT_HEADER

#include "../../CoreModules/MLBase.h"

GRT_BEGIN_NAMESPACE

class MeanShift : public MLBase{
public:
 MeanShift() : MLBase("MeanShift") {
 }
 
 virtual ~MeanShift(){
 
 }
 
 bool search( const VectorFloat &meanStart, const Vector< VectorFloat > &points, const Float searchRadius, const Float sigma = 20.0 ){
 
 //clear the results from any previous search
 clear();
 
 const unsigned int numDimensions = (unsigned int)meanStart.size();
 const unsigned int numPoints = (unsigned int)points.size();
 const Float gamma = 1.0 / (2 * SQR(sigma) );
 unsigned int iteration = 0;
 VectorFloat numer(2,0);
 VectorFloat denom(2,0);
 VectorFloat kernelDist(2,0);
 Float pointsWithinSearchRadius = 0;
 
 mean = meanStart;
 VectorFloat lastMean = mean;
 
 //Start the search loop
 while( true ){
 
 //Reset the counters
 pointsWithinSearchRadius = 0;
 std::fill(numer.begin(),numer.end(),0);
 std::fill(denom.begin(),denom.end(),0);
 std::fill(kernelDist.begin(),kernelDist.end(),0);
 
 //Update the numerator and denominator for points that are with the search radius
 for(unsigned int i=0; i<numPoints; i++){
 
 //Compute the distance of the current point to the mean
 Float distToMean = euclideanDist( mean, points[i] );
 
 //If the point is within the search radius then update numer and denom
 if( distToMean < searchRadius ){
 
 for(unsigned int j=0; j<numDimensions; j++){
 kernelDist[j] = gaussKernel( points[i][j], mean[j], gamma );
 numer[j] += kernelDist[j] * points[i][j];
 denom[j] += kernelDist[j];
 }
 
 pointsWithinSearchRadius++;
 }
 }
 
 //Update the mean
 Float change = 0;
 for(unsigned int j=0; j<numDimensions; j++){
 
 mean[j] = numer[j] / denom[j];
 
 change += grt_sqr( mean[j] - lastMean[j] );
 
 lastMean[j] = mean[j];
 }
 change = grt_sqrt( change );
 
 trainingLog << "iteration: " << iteration;
 trainingLog << " mean: ";
 for(unsigned int j=0; j<numDimensions; j++){
 trainingLog << mean[j] << " ";
 }
 trainingLog << " change: " << change << std::endl;
 
 if( change < minChange ){
 trainingLog << "min changed limit reached - stopping search" << std::endl;
 break;
 }
 
 if( ++iteration >= maxNumEpochs ){
 trainingLog << "max number of iterations reached - stopping search." << std::endl;
 break;
 }
 
 }
 numTrainingIterationsToConverge = iteration;
 trained = true;
 
 return true;
 }
 
 VectorFloat getMean() const {
 return mean;
 }
 
 Float gaussKernel( const Float &x, const Float &mu, const Float gamma ){
 return exp( gamma * grt_sqr(x-mu) );
 }
 
 Float gaussKernel( const VectorFloat &x, const VectorFloat &mu, const Float gamma ){
 
 Float y = 0;
 const UINT N = x.getSize();
 for(UINT i=0; i<N; i++){
 y += grt_sqr(x[i]-mu[i]);
 }
 return exp( gamma * y );
 }
 
 Float euclideanDist( const VectorFloat &x, const VectorFloat &y ){
 
 Float z = 0;
 const UINT N = x.getSize();
 for(UINT i=0; i<N; i++){
 z += grt_sqr(x[i]-y[i]);
 }
 return sqrt( z );
 }
 
protected:

 VectorFloat mean;
 
};

GRT_END_NAMESPACE

#endif //GRT_MEAN_SHIFT_HEADER