Optflow C++

Here is a slimmed down mainline C++ code that uses Seppo Pulkkinen's optflow library. This library uses CImg internally.

#include "CImg_config.h"
#include <CImg.h>
#include <sstream>
#include <string>

#include "DenseVectorFieldIO.h"
#include "DualDenseMotionExtractor.h"
#include "PyramidalLucasKanade.h"
#include "SparseVectorFieldIO.h"
#include "VectorFieldIllustrator.h"

using namespace cimg_library;

int main() {

  CImg< unsigned char > I1("../examples/test1.png");
  CImg< unsigned char > I2("../examples/test2.png");

  const int W = I1.dimx();
  const int H = I1.dimy();
  CImg< unsigned char > I1_smoothed;
  CImg< unsigned char > I2_smoothed;
  CImg< unsigned char > motionImageF(W, H, 1, 3);
  CImg< double > VF, VB;

  I1_smoothed = I1.get_channel(0);
  I2_smoothed = I2.get_channel(0);

  motionImageF.get_shared_channel(0) = I1_smoothed * 0.75;
  motionImageF.get_shared_channel(1) = I1_smoothed * 0.75;
  motionImageF.get_shared_channel(2) = I1_smoothed * 0.75;

  I1_smoothed.blur(3.0, 3.0, 3.0);
  I2_smoothed.blur(3.0, 3.0, 3.0);

  DenseMotionExtractor* e = new PyramidalLucasKanade(8,3,0.0025,0.0,4,true);
  e->compute(I1_smoothed, I2_smoothed, VF, VB);
  printf("%f\n",VF[100,100]);


  return 0;
}

To compile drop this file under lib, run make, create the so, then compile as

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:.
/usr/bin/c++ -L. -Doptflow_EXPORTS -fPIC -I. -Wall -O2 -frounding-math \
-loptflow -o main main.cpp

Mumford on Math

"Mathematicians believe in this Platonic universe in that, there is a pre-existing bunch of facts which are true and you never invent anything, you are discovering".

Optical Flow, Lucas Kanade in Python

Following is the Lucas Kanade optical flow algorithm in Python. We used it successfully on two png images, as well as through OpenCV to follow a point in successive frames. More details are at Github.

import numpy as np
import scipy.signal as si
from PIL import Image

def gauss_kern():
   h1 = 15
   h2 = 15
   x, y = np.mgrid[0:h2, 0:h1]
   x = x-h2/2
   y = y-h1/2
   sigma = 1.5
   g = np.exp( -( x**2 + y**2 ) / (2*sigma**2) );
   return g / g.sum()

def deriv(im1, im2):
   g = gauss_kern()
   Img_smooth = si.convolve(im1,g,mode='same')
   fx,fy=np.gradient(Img_smooth)  
   ft = si.convolve2d(im1, 0.25 * np.ones((2,2))) + \
       si.convolve2d(im2, -0.25 * np.ones((2,2)))
                 
   fx = fx[0:fx.shape[0]-1, 0:fx.shape[1]-1]  
   fy = fy[0:fy.shape[0]-1, 0:fy.shape[1]-1];
   ft = ft[0:ft.shape[0]-1, 0:ft.shape[1]-1];
   return fx, fy, ft

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as si
from PIL import Image
import deriv
import numpy.linalg as lin

def lk(im1, im2, i, j, window_size) :
   fx, fy, ft = deriv.deriv(im1, im2)
   halfWindow = np.floor(window_size/2)
   curFx = fx[i-halfWindow-1:i+halfWindow,
              j-halfWindow-1:j+halfWindow]
   curFy = fy[i-halfWindow-1:i+halfWindow,
              j-halfWindow-1:j+halfWindow]
   curFt = ft[i-halfWindow-1:i+halfWindow,
              j-halfWindow-1:j+halfWindow]
   curFx = curFx.T
   curFy = curFy.T
   curFt = curFt.T

   curFx = curFx.flatten(order='F')
   curFy = curFy.flatten(order='F')
   curFt = -curFt.flatten(order='F')
 
   A = np.vstack((curFx, curFy)).T
   U = np.dot(np.dot(lin.pinv(np.dot(A.T,A)),A.T),curFt)
   return U[0], U[1]