images/optFlowHorn.m: true coarse-to-fine implementation that warps I…

…1 in each phase according to current estimate of flow
kkoolerter · Aug 10, 2012 · 72d4c86 · 72d4c86
1 parent 42ab16f
commit 72d4c86
Show file tree

Hide file tree

Showing 2 changed files with 23 additions and 13 deletions.
diff --git a/images/optFlowHorn.m b/images/optFlowHorn.m
@@ -39,13 +39,17 @@
   scale=2^(nScales-s); h1=round(h/scale); w1=round(w/scale);
   if( scale==1 ), I1b=I1; I2b=I2; else
     I1b=imResample(I1,[h1 w1]); I2b=imResample(I2,[h1 w1]); end
-  % smooth images
-  I1b=convTri(I1b,smooth); I2b=convTri(I2b,smooth);
   % initialize Vx,Vy or upsample from previous scale
   if(s==1), Vx=zeros(h1,w1,'single'); Vy=Vx; else
     Vx=imResample(Vx,[h1 w1])*2; Vy=imResample(Vy,[h1 w1])*2; end
+  % transform I1b according to current estimate of Vx and Vy
+  if(s>1), I1b=imtransform2(I1b,[],'pad','none',...
+      'vs',-double(Vx),'us',-double(Vy)); end
+  % smooth images
+  I1b=convTri(I1b,smooth); I2b=convTri(I2b,smooth);
   % run optical flow on current scale
-  [Vx,Vy]=optFlowHorn1(I1b,I2b,Vx,Vy,alpha,nIter);
+  [Vx1,Vy1]=optFlowHorn1(I1b,I2b,alpha,nIter);
+  Vx=Vx+Vx1; Vy=Vy+Vy1;
 end
 
 % show quiver plot on top of reliab
@@ -56,7 +60,7 @@
 
 end
 
-function [Vx,Vy] = optFlowHorn1( I1,I2,Vx,Vy,alpha,nIter )
+function [Vx,Vy] = optFlowHorn1( I1, I2, alpha, nIter )
 % compute derivatives (averaging over 2x2 neighborhoods)
 A00=shift(I1,0,0); A10=shift(I1,1,0);
 A01=shift(I1,0,1); A11=shift(I1,1,1);
@@ -71,10 +75,10 @@
 Z=1./(alpha*alpha + Ex.*Ex + Ey.*Ey);
 % iterate updating Ux and Vx in each iter
 if( 1 )
-  Vx=shift(Vx,0,0); Vy=shift(Vy,0,0);
-  optFlowHornMex(Vx,Vy,Ex,Ey,Et,Z,nIter);
+  [Vx,Vy]=optFlowHornMex(Ex,Ey,Et,Z,nIter);
   Vx=Vx(2:end-1,2:end-1); Vy=Vy(2:end-1,2:end-1);
 else
+  Vx=zeros(size(I1),'single'); Vy=Vx;
   for i = 1:nIter
     Mx=.25*(shift(Vx,-1,0)+shift(Vx,1,0)+shift(Vx,0,-1)+shift(Vx,0,1));
     My=.25*(shift(Vy,-1,0)+shift(Vy,1,0)+shift(Vy,0,-1)+shift(Vy,0,1));

diff --git a/images/private/optFlowHornMex.cpp b/images/private/optFlowHornMex.cpp
@@ -42,21 +42,27 @@ void optFlowHorn( float *Vx, float *Vy, const float *Ex, const float *Ey,
   delete [] Vx0; delete [] Vy0;
 }
 
-// optFlowHornMex(Vx,Vy,Ex,Ey,Et,Z,nIter); operates in place on Vx,Vy
+// [Vx,Vy]=optFlowHornMex(Ex,Ey,Et,Z,nIter); - helper for optFlowHorn
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
-  size_t h, w, nIter; float *Is[6];
+  size_t h, w, nIter; float *Is[4], *Vx, *Vy;
 
   // Error checking on arguments
-  if( nrhs!=7 ) mexErrMsgTxt("Seven inputs expected.");
-  if( nlhs>0 ) mexErrMsgTxt("NO outputs expected.");
+  if( nrhs!=5 ) mexErrMsgTxt("Five inputs expected.");
+  if( nlhs!=2 ) mexErrMsgTxt("Two outputs expected.");
   h = mxGetM(prhs[0]); w = mxGetN(prhs[0]);
-  for( int i=0; i<6; i++ ) {
+  for( int i=0; i<4; i++ ) {
     if(mxGetM(prhs[i])!=h || mxGetN(prhs[i])!=w) mexErrMsgTxt("Invalid dims.");
     if(mxGetClassID(prhs[i])!=mxSINGLE_CLASS) mexErrMsgTxt("Invalid type.");
     Is[i] = (float*) mxGetData(prhs[i]);
   }
-  nIter = (int) mxGetScalar(prhs[6]);
+  nIter = (int) mxGetScalar(prhs[4]);
+
+  // create output matricies
+  plhs[0] = mxCreateNumericMatrix(int(h),int(w),mxSINGLE_CLASS,mxREAL);
+  plhs[1] = mxCreateNumericMatrix(int(h),int(w),mxSINGLE_CLASS,mxREAL);
+  Vx = (float*) mxGetData(plhs[0]);
+  Vy = (float*) mxGetData(plhs[1]);
 
   // run optical flow
-  optFlowHorn(Is[0],Is[1],Is[2],Is[3],Is[4],Is[5],int(h),int(w),nIter);
+  optFlowHorn(Vx,Vy,Is[0],Is[1],Is[2],Is[3],int(h),int(w),nIter);
 }