[2D] Optical flow

This is a function for optical flow in Matlab, using Lucan-kanande/Horn-Schunck optical flow estimation. It includes a sample code and the function itself.

Assumptions behind the Lucas-Kanade and Horn-Schunck Optical Flow Estimation methods

Both the Horn Schunck and Lucas-Kanade method assumes that brightness constancy does not
change over time. Also, the Horn Schunck method assumes that the flowfield is globally smooth (or
that neighboring velocities are nearly identical). The Lucas-Kanade method assumes that the velocity
is locally constant, and neighboring points belong to the same patch that have similar motion.

Differences between Lucas-Kanade and Horn-Schunck methods

The Lucas-Kanade method is a sparse/local method , while the Horn-Schunck method is a
dense/global one. The Lucas-Kanade would obtain less noise compared to the dense method as
Horn-Schunck. The Horn-Schunck is a dense one, which means it will have denser flow fields. The
pixel displacement of the Lucas-Kanade is constant within a small vicinity, whilst that of Horn-
Schunck is smooth over the time domain.

The threshold to constrain Lucas-Kanade Optical Flow estimation

This is done to ensure there is a numerical stability. We obtain a useful measure of reliability of
motion for each computation. Fields that have a fairly undefined structure are excluded from
potential error because of this (flow while there seems to be no flow).

Tradeoffs associated with a smaller or larger window size

The function takes quite a while, mostly due to the fact that I did not put much effort in making it fast
and efficient. Because of that, a larger window size makes the computation take even longer because
with a larger window size it means it has more values to process.
A larger window will also produce a less ‘sharp’ directional image. Because the windows are larger,
more ‘windows’ will include the parts that are actually moving. This means that the area with
direction will get bigger. On the other hand, a larger window will include the parts that might be
missed by smaller windows in terms of brightness constancy. For instance, a small part (for instance
the back of the dog) is quite large and black, and for small windows this would not seem to move. Big
windows on the other hand will include these potential problem area’s as well, so that these areas
are also moving as they are. So there is a tradeoff between local accuracy and robustness when
choosing this window size. On the follow page I have shown window sizes of 2,5,10,20,30 to illustrate
the differences. Also see that i use the “quiver” function below.

Different used window-sizes


The left and richt image (see the dog moving?)

Example code

  1. %Tim Zaman, TU Delft ME, 12-12-2010
  2. clc
  3. clear all
  4. close all
  5. windowSize=5;
  6. im1=imread('image1.bmp'); %Read image 1
  7. im1=double(rgb2gray(im1))/255; %Convert to gray & range 0-1
  8. im2=imread('image2.bmp'); %Read image 2
  9. im2=double(rgb2gray(im2))/255; %Convert to gray & range 0-1
  10. [u,v]=optical_flow(im1,im2,windowSize);
  11. quiver_uv(u,v)

Optical flow function Matlab

  1. function [u,v]=optical_flow(im1,im2,windowSize)
  2. %Compute differences
  3. x_c=im1(1:end-1,2:end)-im2(1:end-1,1:end-1);
  4. y_c=im2(2:end,1:end-1)-im1(1:end-1,1:end-1);
  5. t_c=im2(1:end-1,1:end-1)-im1(1:end-1,1:end-1);
  6. %initialize for speed
  7. u = zeros(size(x_c));
  8. v = u;
  9. for x = 1:size(x_c,1)-windowSize %fpr all x
  10. for y = 1:size(x_c,2)-windowSize
  11. %Get the windows of the dimensions
  12. win_x=imcrop(x_c,[x y windowSize windowSize]);
  13. win_y=imcrop(y_c,[x y windowSize windowSize]);
  14. win_t=imcrop(t_c,[x y windowSize windowSize]);
  15. %Convert windows to vectors to produce A for solving
  16. later
  17. A = [win_x(:) win_y(:)];
  18. %Compute threshold t (smallest eigenvalue of A'A)
  19. th=min(eig(A'*A));
  20. %Optical flow is only valid in regions with t<0.01 and
  21. rank =2
  22. %if true, then it should not be computed
  23. if th<0.01 || rank(A'*A)~=2
  24. unk=[0 0];
  25. else
  26. unk = pinv(A'*A)*A'*win_t(:); %compute unkown or in
  27. lecture 'x'
  28. end
  29. %put in u and v
  30. u(x,y)=unk(1);
  31. end
  32. end

Quiver

  1. function quiver_uv (u, v)
  2. % Resize u and v so we can actually see something in the quiver plot
  3. scalefactor = 50/size(u,2);
  4. u_ = scalefactor*imresize(u,scalefactor,'bilinear');
  5. v_ = scalefactor*imresize(v,scalefactor,'bilinear');
  6. % Run quiver taking into account matlab coordinate system quirkiness
  7. % and scaling the magnitude of (u,v) by 2 so it is more visible.
  8. quiver(u_(end:-1:1,:),-v_(end:-1:1,:),2);
  9. axis('tight');

Tim Zaman

MSc Biorobotics. Specialization in computer vision and deep learning. Works at NVIDIA.

17 Responses

  1. meda says:

    Hi,
    dankjewel for your post!
    You showed how to compute, given two images, the velocity components (u,v). These components are not integers.
    Here’s my question: how can I register the second image according to the optical flow?
    I have seen that warping techniques (via bilinear, cubic, etc interpolation) could come in handy but are there other techniques/approaches? I would like to correct movement using optical flow in order to obtain a steady scene in each frame.
    Thank you in advance (and compliments for your research work)

  2. Nits says:

    I do not understand, how to mark a region and see it moving in another image below, please help me.

  3. Nits says:

    hi,

    i wanna help

  4. kahi says:

    salut , je voulais m’envoyer le programme de l’implimentation de la methode de horn & schunck et celle de lucas-kanade de flot optique sur matlab pour l’estimation et la detection de mouvement si vous pouvez slvp

  5. wajid says:

    that cleared my point of view about optical flow, since I am new to matlab I am getting problem at line
    if th<0.01 || rank(A'*A)~=2

    Matlab gives error "Subscript indices must either be real positive integers or logicals."

    can anyone guide me please

  6. prashanth says:

    hi,
    this code helped me a lot…
    i am doing background subtraction using this optical flow method..
    can u please help me with this

  7. Ben Katz says:

    I have the same problem as Victor. For some unknown reason Tim is NOT bothering to reply to Victor!

    • Tim Zaman says:

      Show some respect please. Instead of whining that it doesn’t work *for you*, try being more academic and see what the problem is and try to solve it. Don’t be lazy.
      And its not my code that’s the problem, it’s you input. If you know how the code works you can see that your image should be an integer multiple of the ‘windowsize’. In that case the problem should not occur. Don’t try to copy my code 1:1, its just for reference, if you would like to have me make a readymade solution for you, try donating —>

  8. Victor says:

    Hi guys,

    I can’t make work the optical flow code because the ‘end’ statement is missing for both ‘for’ statements.

    Also this line is causing a problem:

    th=min(eig(A’*A));

    MTIMES is not fully supported for integer classes. At least one input must be scalar.

    since the matrices are of different sizes changing the operator for .* won’t work neither.

    I’ll appreciate if those who got it working could give me a hand.

    Many thanks

  9. oo says:

    very good result.

    • kahi says:

      slt, est ce que vous pouvez m’envoyer se programme sur ma boite slvp j’ai besoin de lui merci

  10. Runlada says:

    Thank you for your code. It ‘s very good. My project is clear.

  11. ting says:

    Your code is greet and i am able to run it.
    However, is there have any documentation or pdf. file or link which can make me understand your code and the idea you wrote this code? thanks

  12. Sasan says:

    Hi, why can’t I use the optical flow method?
    Fo example when I enter Your code, in command line10 it gives me this error:
    .??? Undefined function or method ‘optical_flow’ for arguments of type ‘double’.

    Is there an installation I have to do before using optical flow?
    Plz help me it’s urgent

    • Tim says:

      dont you see the matlab function called optical flow above here?? use that as a function (= place that text in an m-file)
      dont you know how to make functions in matlab??

      • sasan says:

        do u mean File/new/function???? if yes, I paste your code there. but it keeps giving me the same error.

        Is this what you mean?

      • Tim says:

        are you aure the function itself is defined well and has no errors? make sure its in your workspace and open it in matlab and check for red errors..