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 comment »