I recently read about the simple but brilliantly clever improvement to RANSAC, called PROSAC. Maybe it’s not really all that brilliant, because as far as I can tell, it is not picked up and built upon in the authors’ subsequent RANSAC-related papers. But I think it’s cool, and easy enough that I want to play around with adding it to the OpenCV baseline, maybe even eventually get it solid enough to request it to be added officially.

I don’t want to try to explain or reiterate the entire paper, the general idea is, if your prospective matches are usefully ordered (say, small descriptor distances up front, or larger 1st/2nd-best-dist-ratio up front), then it makes no sense for RANSAC to give every element the same chance to be in a sample that defines a model. Why not give priority to the (probably) better elements at the front of the list? So authors Chum&Matas derive a meaningful schedule for drawing samples with max element m, then m+1, then m+2, etc. So instead of simply RANSAC=RANdom SAmple Consensus, we have PROgressive SAmple Consensus. Chum&Matas also provide more useful stopping criteria than RANSAC alone, which I haven’t gotten into yet.

So without further ado, here’s my first stab at code for refactoring a RANSAC implementation over to PROSAC (no thanks to wordpress for not maintaining indentation of my paste!):

#include <stdlib.h> // rand() #include <math.h> // ceil() #ifdef DEBUG #include <iostream> using std::cout; using std::endl; #endif typedef unsigned int uint; #ifdef DEBUG #define DPRINT(x) cout << "\t" << x << endl; #else #define DPRINT(x) #endif class ProsacHelper { public: ProsacHelper(uint m_, uint N_, uint TN_=200000) : m(m_) // sample size , N(N_) // full size , TN(TN_) // max iterations , poolmax(m-1) // initial poolmax is m-1 , poolit(0) // initial iteration is 0 , si(0) // be ready to serve first element of first sample { // if (m<2 || m>100) throw... DPRINT("m = " << m); DPRINT("N = " << N); DPRINT("TN = " << TN); DPRINT("poolmax = " << poolmax); // initialize current T(n=m) by direct computation T = TN; for (uint i=0; i<=m-1; ++i) { T *= m-i; T /= N-i; } DPRINT("T = " << T); // spend only 1 it in first pool (sample=0...m-1) poolits=1; for (uint i=0; i<m; ++i) sample[i] = i; DPRINT("3:0 0 1 2 3"); } uint randomSample() { uint ret = sample[si]; // this is what will be returned si++; // advance within current sample if (si==m) { // we have exhausted current sample of m values for iteration. Time // to advance to the next iteration, and perhaps the next poolmax, and // draw the sample for the next iteration poolit++; if (poolit >= poolits && // time to expand the poolmax poolmax < N) { // ...unless we can't poolmax++; DPRINT("poolmax = " << poolmax); // eq (3) n=poolmax, n+1=poolmax+1 DPRINT("curT = " << T); double nxtT = T * (poolmax+1) / (poolmax+1-m); DPRINT("nxtT = " << nxtT); double dpoolits = nxtT - T; poolits = ceil(dpoolits); DPRINT("its with max= "<<poolmax<<": "<< dpoolits<<" --> "<<poolits); poolit=0; T = nxtT; } // now that we know what the poolmax is, // create the next sample uint n_unique_elts = (poolmax >= N ? m : m-1); // loop to populate elements 0..m-2 from 0..n-2 for (uint i=0; i<n_unique_elts; ++i) { bool unique=false; // draw random r until different from all previous do { sample[i]=rand()%poolmax; // 0..poolmax-1 unique=true; for (uint j=0; j<i; ++j) { if (sample[i] == sample[j]) { unique=false; break; } } } while (!unique); } if (poolmax < N) { // last element must be end of pool sample[m-1] = poolmax; } si=0; // reset sample counter for next call DPRINT(poolmax<<":"<<poolit<<": " << sample[0] << " " << sample[1] << " " << sample[2] << " " << sample[3]); } return ret; } protected: uint m; // sample size: homographies=4, fundamental matrices=7 or 8 uint N; // full number of items uint poolmax; // current pool max (pool = 0..poolmax) uint poolit; // iteration within current pool size uint TN; // max iterations double T; // T(n)=avg # regular RANSAC seeds out of TN with max<=n-1 uint poolits; // # iterations to spend in current pool; i=0..poolits-1 uint sample[100]; // cache holding the current sample of m indices uint si; // current position within the current sample };

So the whole point is that you construct a ProsacHelper with parameters relevant to your case (m=4 for fitting a homography, or 7 or 8 for fitting a fundamental matrix; and N=number of elements; you can probably leave TN default), and then call randomSample() repeatedly. It remembers its state within each sample of m (with member si), and (as long as you don’t call it to use the random indices for other purposes) it will give you m-samples according to the schedule laid out by PROSAC.

In the OpenCV calib3d module, source file modelest.cpp, method CvModelEstimator2::getSubset(), the line

idx[i] = idx_i = cvRandInt(&rng)%count;

would be replaced by

idx[i] = idx_i = prosac_helper.randomSample();

–that is, assuming a ProsacHelper object could be set up beforehand and persisted across calls to getSubset(). I hope to show how this works out in a next post.

I did a few things a little differently from the paper. Member ‘poolsize’ is n, the size of the current pool of samples (0..n-1, or more specifically, n-1 and m-1 from 0..n-2). Rather than maintain t=global iteration number and T’n=iteration to bump up n=poolsize, I maintain only T'(n) – T'(n-1) = ceil(T(n) – T(n-1)) = ‘poolits’, and ‘poolit’ is my counter from 0..poolits-1.

I did spend time in there to guarantee that each m-sample does not repeat itself. I figured that this time is probably overshadowed by the time required to make sure that random samples are not degenerate (i.e. 3/4 collinear when fitting a homography).

Test code like this:

ProsacHelper prosac(4, 100); uint sample[4]; for (int it=0; it<100; ++it) { for (int mi=0; mi<4; ++mi) sample[mi] = prosac.randomSample(); cout << it << ": " << sample[0] << " " << sample[1] << " " << sample[2] << " " << sample[3] << endl; }

creates output like this:

0: 0 1 2 3 1: 3 2 1 4 2: 0 3 1 5 3: 0 3 1 6 4: 2 1 5 6 5: 0 3 1 7 6: 2 3 0 7 7: 7 5 6 8 8: 2 6 3 8 9: 3 7 1 8 ...

You can see that, since T(n) increases by less than one the first few times, poolmax=3,4,5 each only get one iteration, but then as T(n) grows more quickly, poolmax=6,7 get 2 iterations each, poolmax=8 gets 3 iterations, and it grows from there.

Filed under: General |

## Leave a Reply