Hoagies&Stogies: Eastern Orthodox

Thanks to Rob Hickok for organizing a loooong overdue H&S (and sorry from me for forgetting for so long to post about the audio)!

Way back on Oct 24, we gathered at St Gregory of Nyssa EO church, some of us attended their vespers service, we shared a meal with their members, and had a H&S discussion after with father Simeon Corona, and Ben Rochester of Pilgrim Presbyterian.

Audio was uploaded shortly after the event, and can be obtained from here (again, sorry about the delay)!

PROSAC for OpenCV (2)

I made some changes to my code from last time. I won’t repost the whole thing, just snippets.

I removed typedef uint int and changed all uint to int, since the OpenCV coding standard seems not to be finicky about using unsigned ints.

In preparation for integrating into OpenCV (3.0), while retaining the ability to test standalone, I defined a macro for random integers:

// for testing outside opencv
#include <stdlib.h> // rand()
#include <math.h> // ceil()
#define rand_uint_less_than(limit) rand()%limit
// for use inside opencv
#define rand_uint_less_than(limit) cvRandInt(&rng) % limit

In order to easily switch back and forth between PROSAC and RANSAC for comparative testing, I moved the contents of the constructor to an initPROSAC() method, and created a separate initRANSAC() method. All the initRANSAC() method does is set poolmax=N-1; which causes the preexisting randomSample() method to behave in a RANSAC way (i.e. sample m from the entire pool). The new constructor is empty. Perhaps a better pattern would be to have a SacHelper base class with the randomSample() method, and derive RansacHelper and ProsacHelper each with their specific constructors. Maybe I’ll do that at some point.

My test driver was updated to exercise both capabilities, like this:

ProsacHelper pro_or_ran;
#define EM 4
uint sample[EM];
for (int which=0; which<2; ++which) {
  switch(which) {
   case 0: pro_or_ran.initRANSAC(EM,100); break;
   case 1: pro_or_ran.initPROSAC(EM,100); break; }
  for (int it=0; it<100; ++it) {
    for (int mi=0; mi<EM; ++mi)
      sample[mi] = pro_or_ran.randomSample();
      cout << it << ":";
      for (int mi=0; mi<EM; ++mi)
        cout << " " << sample[mi];
      cout << endl;

And output looks something like this:

 0: 83 86 77 15
 1: 93 35 86 92
 2: 49 21 62 27
 3: 90 59 63 26
 4: 40 26 72 36
 5: 11 68 67 29
 0: 0 1 2 3
 1: 3 1 0 4
 2: 0 4 2 5
 3: 0 2 3 6
 4: 2 0 5 6
 5: 5 4 6 7
 6: 2 6 3 7
 7: 5 4 6 8
 8: 5 0 4 8
 9: 7 4 5 8
 10: 0 4 3 9

I have cloned the git repository of OpenCV (3.0) and the opencv_extra with testdata, “cmake . ; make” was sufficient to build perfectly first time, and the unit test suites opencv_test_features2d and opencv_test_calib3d are passing. Next post I hope to show first results of slapping this class into the calib3d source code, into RANSACPointSetRegistrator::getSubset() in function modules/calib3d/src/ptsetreg.cpp. This class is used to conduct RANSAC for cv::findHomography() and cv::findFundamentalMatrix(). Hopefully I can slap my ProsacHelper in and not break the test suite.


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;

typedef unsigned int uint;

#ifdef DEBUG
#define DPRINT(x) cout << "\t" << x << endl;
#define DPRINT(x)

class ProsacHelper {
 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)
   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

     if (poolit >= poolits && // time to expand the poolmax
         poolmax < N) { // ...unless we can't
       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);
       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
         for (uint j=0; j<i; ++j) {
           if (sample[i] == sample[j]) {
       } 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;

 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.

Eulogy for Dave

My friend Dave died this week, a victim of an unexpected stroke — tragically young (early 50s?)

Although I hadn’t seen Dave much in the last two years since I moved, there was a period of maybe 5 years that we met regularly for coffee every Monday morning, with a few other guys, before work. During that time, we all suffered with him through his divorce, and rejoiced with him when he found, wooed, and married Susan. And we got to hear a lot about progress in the lives of his two daughters, of whom he was so proud.

Dave was a manly man; he knew about manly things like cars and preventative maintenance. He also had sophisticated tastes in wine, food, and jazz. On the other hand, he also had a corny sense of humor, epitomized by his love for Ernest movies (“you know what ah mean, Vern?” Ernest Goes to Camp, Ernest Goes to Jail, …).

Dave loved my kids. He taught two of them through 1st-2nd grade sunday school, and he was always telling me or Tina some story about something they said or did in class that amazed him or made him laugh. I know he told other parents about things their kids did in class, but I like to think he loved my boys the most.

I don’t feel any compunction to say only nice things about Dave; he surely had his warts like the rest of us, but I didn’t see much evidence of them. When he was going through his divorce, he told me what he thought his faults were in the marriage: being distant and curt when laid up with back injuries, getting impatient and angry. (He often encouraged me not to take my wife for granted.) He might have had a bit of a temper, but I only ever saw him calmly deal with frustrating issues that would have enraged me.

But I don’t need to have seen his faults to know that Dave was a sinner; the bible tells me so, and he confessed it himself. Dave might have been less of a sinner than most of us, but he was still enough of a sinner to earn and deserve the eternal wrath of God in hell. Thankfully, God gave him the grace to be able to claim Christ’s sacrifice, bearing that punishment on his behalf. That’s the only reason Dave’s untimely death is not a total tragedy; he will be resurrected on the last day. Meanwhile, a lot of people will be missing him.

Memorial services at New Life PCA, la Mesa, Sat Apr 11, 2pm.

Mini-Reunion III

Another two years, another Mini-Reunion. Had a great time with the Gibsons again. They offered a great dinner of pork (but no mashed potatoes), various family members rushing in and out for babysitting and other meetings, family music, games, breakfast of al dente rice pudding before going off to my two-hour meeting and flying back home!

Come and Eat

Sorry, this is not an insightful meditation on John 6, or Rev 3:20, or Isaiah 55; rather, I want you to come, buy and eat, with money! I’m a Calvinist Cadet Corps Counselor (and webmaster for the Southern California Council). We’re raising money in support of the 2010 (Inter)National Counselor’s Convention, which our council is hosting here in San Diego this summer. If you’re in the area, please come out next Sat. 24th and come have pancakes with us; or (and!) come out Mon. 26th and eat at Souplantation for us.

Here’s a flyer for the pancake breakfast, if you want more information.

Here’s a flyer for Souplantationwhich you must bring to Souplantation if we are to redeem our 20%. (Note: this is only at the Fletcher Parkway Souplantation, from 5-8pm on Mon Apr 26th)

Thanks for your support!

Sir #1

Yesterday afternoon there were some SCA-types in a local park, so I brought #1 round, and he had a blast learning about (and trying on!) a wide variety of helmets, gauntlets, chain mail, shields, coats, swords, etc. They were super-nice and glad to show off their toys, and even took a few pictures for us! Click through to flickr, and look at the next couple of pics in the photostream…