[an error occurred while processing this directive]
Main Page | Modules | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members | File Members | Related Pages

hugin_base/vigra_ext/ransac.h

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00028 #ifndef _RANSAC_H_
00029 #define _RANSAC_H_
00030 
00031 #include <set>
00032 #include <vector>
00033 #include <stdlib.h>
00034 #include <cstring>
00035 #include <math.h>
00036 #include <ctime>
00037 
00038 #include <boost/random/mersenne_twister.hpp>
00039 #include <boost/random/uniform_int.hpp>
00040 #include <boost/random/variate_generator.hpp>
00041 
00042 //#include "ParameterEsitmator.h"
00043 
00073 class Ransac {
00074 
00075 public:
00092         template<class Estimator, class S, class T>
00093         static int compute(S & parameters, 
00094                            const Estimator & paramEstimator ,
00095                            const std::vector<T> &data, 
00096                            double desiredProbabilityForNoOutliers,
00097                            double maximalOutlierPercentage);
00098 
00099 
00122         template<class Estimator, class S, class T>
00123         static int compute(S &parameters, 
00124                            const Estimator & paramEstimator ,
00125                            const std::vector<T> &data);
00126 
00127 private:
00128 
00132     static unsigned int choose(unsigned int n, unsigned int m);
00133 
00134     template<class Estimator, class T>
00135     static void computeAllChoices(const Estimator & paramEstimator,
00136                                   const std::vector<T> &data,
00137                                   int numForEstimate,
00138                                   short *bestVotes, short *curVotes,
00139                                   int &numVotesForBest, int startIndex,
00140                                   int n, int k, int arrIndex, int *arr);
00141 
00142 
00143     template<class Estimator, class T, class S>
00144     static void estimate(const Estimator & paramEstimator, const std::vector<T> &data,
00145                          int numForEstimate,
00146                          short *bestVotes, short *curVotes,
00147                          int &numVotesForBest, int *arr);
00148 
00149     class SubSetIndexComparator 
00150     {
00151     private:
00152         int m_length;
00153         public:
00154             SubSetIndexComparator(int arrayLength) : m_length(arrayLength)
00155             {}
00156 
00157         bool operator()(const int *arr1, const int *arr2) const 
00158         {
00159             for(int i=0; i<m_length; i++)
00160                 if(arr1[i] < arr2[i])
00161                     return true;
00162                 return false;                   
00163         }
00164     };
00165 
00166 };
00167 
00168 
00169 /*******************************RanSac Implementation*************************/
00170 
00171 template<class Estimator, class S, class T>
00172 int Ransac::compute(S &parameters,
00173                        const Estimator & paramEstimator,
00174                        const std::vector<T> &data,
00175                        double desiredProbabilityForNoOutliers,
00176                        double maximalOutlierPercentage)
00177 {
00178     unsigned int numDataObjects = (int) data.size();
00179     unsigned int numForEstimate = paramEstimator.numForEstimate();
00180         //there are less data objects than the minimum required for an exact fit, or
00181         //all the data is outliers?
00182     if(numDataObjects < numForEstimate || maximalOutlierPercentage>=1.0) 
00183         return 0;
00184 
00185     std::vector<const T *> exactEstimateData;
00186     std::vector<const T *> leastSquaresEstimateData;
00187     S exactEstimateParameters;
00188     int i, j, k, l, numVotesForBest, numVotesForCur, maxIndex, numTries;
00189     short *bestVotes = new short[numDataObjects]; //one if data[i] agrees with the best model, otherwise zero
00190     short *curVotes = new short[numDataObjects];  //one if data[i] agrees with the current model, otherwise zero
00191     short *notChosen = new short[numDataObjects]; //not zero if data[i] is NOT chosen for computing the exact fit, otherwise zero
00192     SubSetIndexComparator subSetIndexComparator(numForEstimate);
00193     std::set<int *, SubSetIndexComparator > chosenSubSets(subSetIndexComparator);
00194     int *curSubSetIndexes;
00195     double outlierPercentage = maximalOutlierPercentage;
00196     double numerator = log(1.0-desiredProbabilityForNoOutliers);
00197     double denominator = log(1- pow((double)(1.0-maximalOutlierPercentage), (double)(numForEstimate)));
00198     int allTries = choose(numDataObjects,numForEstimate);
00199 
00200     //parameters.clear();
00201 
00202 
00203     numVotesForBest = 0; //initalize with 0 so that the first computation which gives any type of fit will be set to best
00204     
00205     // intialize random generator
00206     boost::mt19937 rng;
00207     // start with a different seed every time.
00208     rng.seed(static_cast<unsigned int>(std::time(0)));
00209     // randomly sample points.
00210     maxIndex = numDataObjects-1;
00211     boost::uniform_int<> distribIndex(0, maxIndex);
00212     boost::variate_generator<boost::mt19937&, boost::uniform_int<> >
00213                              randIndex(rng, distribIndex);  // glues randomness with mapping
00214 
00215 //    srand((unsigned)time(NULL)); //seed random number generator
00216     numTries = (int)(numerator/denominator + 0.5);
00217 
00218                     //there are cases when the probablistic number of tries is greater than all possible sub-sets
00219     numTries = numTries<allTries ? numTries : allTries;
00220 
00221     for(i=0; i<numTries; i++) {
00222         //randomly select data for exact model fit ('numForEstimate' objects).
00223         memset(notChosen,'1',numDataObjects*sizeof(short));
00224         curSubSetIndexes = new int[numForEstimate];
00225 
00226         exactEstimateData.clear();
00227 
00228         maxIndex = numDataObjects-1; 
00229         for(l=0; l<(int)numForEstimate; l++) {
00230             //selectedIndex is in [0,maxIndex]
00231             unsigned int selectedIndex = randIndex();
00232 //            unsigned int selectedIndex = (unsigned int)(((float)rand()/(float)RAND_MAX)*maxIndex + 0.5);
00233             for(j=-1,k=0; k<(int)numDataObjects && j<(int)selectedIndex; k++) {
00234                 if(notChosen[k])
00235                     j++;
00236             }
00237             k--;
00238             exactEstimateData.push_back(&(data[k]));
00239             notChosen[k] = 0;
00240             maxIndex--;
00241         }
00242         //get the indexes of the chosen objects so we can check that this sub-set hasn't been
00243         //chosen already
00244         for(l=0, j=0; j<(int)numDataObjects; j++) {
00245             if(!notChosen[j]) {
00246                 curSubSetIndexes[l] = j+1;
00247                 l++;
00248             }
00249         }
00250 
00251         //check that the sub set just chosen is unique
00252         std::pair< std::set<int *, SubSetIndexComparator >::iterator, bool > res = chosenSubSets.insert(curSubSetIndexes);
00253 
00254         if(res.second == true) { //first time we chose this sub set
00255                                  //use the selected data for an exact model parameter fit
00256             if (!paramEstimator.estimate(exactEstimateData,exactEstimateParameters))
00257                 //selected data is a singular configuration (e.g. three colinear points for 
00258                 //a circle fit)
00259                 continue;
00260             //see how many agree on this estimate
00261             numVotesForCur = 0;
00262             memset(curVotes,'\0',numDataObjects*sizeof(short));
00263             for(j=0; j<(int)numDataObjects; j++) {
00264                 if(paramEstimator.agree(exactEstimateParameters, data[j])) {
00265                     curVotes[j] = 1;
00266                     numVotesForCur++;
00267                 }
00268             }
00269             if(numVotesForCur > numVotesForBest) {
00270                 numVotesForBest = numVotesForCur;
00271                 memcpy(bestVotes,curVotes, numDataObjects*sizeof(short));
00272             }
00273             //update the estimate of outliers and the number of iterations we need
00274             outlierPercentage = 1 - (double)numVotesForCur/(double)numDataObjects;
00275             if(outlierPercentage < maximalOutlierPercentage) {
00276                 maximalOutlierPercentage = outlierPercentage;
00277                 denominator = log(1- pow((double)(1.0-maximalOutlierPercentage), (double)(numForEstimate)));
00278                 numTries = (int)(numerator/denominator + 0.5);
00279                 //there are cases when the probablistic number of tries is greater than all possible sub-sets
00280                 numTries = numTries<allTries ? numTries : allTries;
00281             }
00282         }
00283         else {  //this sub set already appeared, don't count this iteration
00284             delete [] curSubSetIndexes;
00285             i--;
00286         }
00287     }
00288 
00289     //release the memory
00290     std::set<int *, SubSetIndexComparator >::iterator it = chosenSubSets.begin();
00291     std::set<int *, SubSetIndexComparator >::iterator chosenSubSetsEnd = chosenSubSets.end();
00292     while(it!=chosenSubSetsEnd) {
00293     delete [] (*it);
00294         it++;
00295     }
00296     chosenSubSets.clear();
00297 
00298     //compute the least squares estimate using the largest sub set
00299     if(numVotesForBest > 0) {
00300         for(j=0; j<(int)numDataObjects; j++) {
00301             if(bestVotes[j])
00302                 leastSquaresEstimateData.push_back(&(data[j]));
00303         }
00304         paramEstimator.leastSquaresEstimate(leastSquaresEstimateData,parameters);
00305     }
00306     delete [] bestVotes;
00307     delete [] curVotes;
00308     delete [] notChosen;
00309 
00310     return numVotesForBest;
00311 }
00312 /*****************************************************************************/
00313 template<class Estimator, class S, class T>
00314 int Ransac::compute(S &parameters,
00315                     const Estimator & paramEstimator,
00316                     const std::vector<T> &data)
00317 {
00318     unsigned int numForEstimate = paramEstimator.numForEstimate();
00319     std::vector<T *> leastSquaresEstimateData;
00320     int numDataObjects = data.size();
00321     int numVotesForBest = 0;
00322     int *arr = new int[numForEstimate];
00323     short *curVotes = new short[numDataObjects];  //one if data[i] agrees with the current model, otherwise zero
00324     short *bestVotes = new short[numDataObjects];  //one if data[i] agrees with the best model, otherwise zero
00325         
00326     //parameters.clear();
00327 
00328     //there are less data objects than the minimum required for an exact fit
00329     if(numDataObjects < numForEstimate) 
00330         return 0;
00331 
00332     computeAllChoices(paramEstimator,data,numForEstimate,
00333                       bestVotes, curVotes, numVotesForBest, 0, data.size(), numForEstimate, 0, arr);
00334 
00335     //compute the least squares estimate using the largest sub set
00336     if(numVotesForBest > 0) {
00337         for(int j=0; j<numDataObjects; j++) {
00338             if(bestVotes[j])
00339                 leastSquaresEstimateData.push_back(&(data[j]));
00340         }
00341         paramEstimator.leastSquaresEstimate(leastSquaresEstimateData,parameters);
00342     }
00343 
00344     delete [] arr;
00345     delete [] bestVotes;
00346     delete [] curVotes; 
00347 
00348     return numVotesForBest;
00349 }
00350 /*****************************************************************************/
00351 template<class Estimator, class T>
00352 void Ransac::computeAllChoices(const Estimator &paramEstimator, const std::vector<T> &data,
00353                                 int numForEstimate,short *bestVotes, short *curVotes,
00354                                 int &numVotesForBest, int startIndex, int n, int k,
00355                                 int arrIndex, int *arr)
00356 {
00357                       //we have a new choice of indexes
00358   if(k==0) {
00359                 estimate(paramEstimator, data, numForEstimate, bestVotes, curVotes, numVotesForBest, arr);
00360     return;
00361   }
00362 
00363                //continue to recursivly generate the choice of indexes
00364   int endIndex = n-k;
00365   for(int i=startIndex; i<=endIndex; i++) {
00366     arr[arrIndex] = i;
00367     computeAllChoices(paramEstimator, data, numForEstimate, bestVotes, curVotes, numVotesForBest,
00368                                         i+1, n, k-1, arrIndex+1, arr);
00369   }
00370 
00371 }
00372 /*****************************************************************************/
00373 template<class Estimator, class T, class S>
00374 void Ransac::estimate(const Estimator & paramEstimator, const std::vector<T> &data,
00375                       int numForEstimate, short *bestVotes, short *curVotes,
00376                       int &numVotesForBest, int *arr)
00377 {
00378         std::vector<T *> exactEstimateData;
00379         std::vector<S> exactEstimateParameters;
00380         int numDataObjects;
00381         int numVotesForCur;//initalize with -1 so that the first computation will be set to best
00382         int j;
00383 
00384         numDataObjects = data.size();
00385         memset(curVotes,'\0',numDataObjects*sizeof(short));
00386         numVotesForCur=0;
00387 
00388         for(j=0; j<numForEstimate; j++)
00389                 exactEstimateData.push_back(&(data[arr[j]]));
00390         paramEstimator.estimate(exactEstimateData,exactEstimateParameters);
00391                              //singular data configuration
00392         if(exactEstimateParameters.size()==0)
00393                 return;
00394 
00395         for(j=0; j<numDataObjects; j++) {
00396                 if(paramEstimator.agree(exactEstimateParameters, data[j])) {
00397                         curVotes[j] = 1;
00398                         numVotesForCur++;
00399                 }
00400         }
00401         if(numVotesForCur > numVotesForBest) {
00402                 numVotesForBest = numVotesForCur;
00403                 memcpy(bestVotes,curVotes, numDataObjects*sizeof(short));
00404         }
00405 }
00406 /*****************************************************************************/
00407 inline unsigned int Ransac::choose(unsigned int n, unsigned int m)
00408 {
00409         unsigned int denominatorEnd, numeratorStart, numerator,denominator; 
00410         if((n-m) > m) {
00411                 numeratorStart = n-m+1;
00412                 denominatorEnd = m;
00413         }
00414         else {
00415                 numeratorStart = m+1;
00416                 denominatorEnd = n-m;
00417         }
00418         
00419         unsigned int i;
00420         for(i=numeratorStart, numerator=1; i<=n; i++)
00421                 numerator*=i;
00422         for(i=1, denominator=1; i<=denominatorEnd; i++)
00423                 denominator*=i;
00424         return numerator/denominator;
00425 
00426 }
00427 
00428 
00429 #endif //_RANSAC_H_

Generated on Mon Sep 20 01:01:27 2010 for Hugintrunk by doxygen 1.3.9.1