Optical flow in matlab

Implementation based on Horn’s optical flow algorithm.


function [Vx,Vy] = OpticalFlow(images,alpha,iterations)
%// Calculating optical flow of a sequence of images. 
%// images : 3D array that contains a sequence of images. size of images is (imageHeight, imageWidth, frame number)
%// alpha
%// iterations. 
[height,width,frames]=size(images);

%//initialzation of u and v
Vx = zeros(height,width);
Vy = zeros(height,width);

for k = 1:frames-1
    
    % //initialization of Ex Ey and Et 
    Ex = zeros(height-1,width-1,frames-1);
    Ey = zeros(height-1,width-1,frames-1);
    Et = zeros(height-1,width-1,frames-1);
    
    %//calculating Ex Ey and Et in frame k.
    for x = 2:width-1
        for y = 2:height-1
            Ex(y,x,k) = (images(y+1,x+1,k)-images(y+1,x,k)+images(y,x+1,k)...
                -images(y,x,k)+images(y+1,x+1,k+1)-images(y+1,x,k+1)...
                +images(y,x+1,k+1)-images(y,x,k+1))/4;
            
            Ey(y,x,k) = (images(y,x,k)-images(y+1,x,k)+images(y,x+1,k)...
                -images(y+1,x+1,k)+images(y,x,k+1)-images(y+1,x,k+1)...
                +images(y,x+1,k+1)-images(y+1,x+1,k+1))/4;
            
            Et(y,x,k) = (images(y+1,x,k+1)-images(y+1,x,k)+images(y,x,k+1)...
                -images(y,x,k)+images(y+1,x+1,k+1)-images(y+1,x+1,k)...
                +images(y,x+1,k+1)-images(y,x+1,k))/4;
        end
    end

    for nn = 1:iterations
        for x = 2:width-1
            for y = 2:height-1
                
                Vxbar = (Vx(y-1,x)+Vx(y,x+1)+Vx(y+1,x)+Vx(y,x-1))/6+...
                     (Vx(y-1,x-1)+Vx(y-1,x+1)+Vx(y+1,x+1)+Vx(y+1,x-1))/12;
                
                Vybar = (Vy(y-1,x)+Vy(y,x+1)+Vy(y+1,x)+Vy(y,x-1))/6+...
                    (Vy(y-1,x-1)+Vy(y-1,x+1)+Vy(y+1,x+1)+Vy(y+1,x-1))/12;

                %// chapter 12 of Horn's paper
                temp = (Ex(y,x,k)*Vxbar+Ey(y,x,k)*Vybar+Et(y,x,k))/(alpha^2 + Ex(y,x,k)^2 + Ey(y,x,k)^2);
                %// update u and v 
                Vx(y,x) = Vxbar-Ex(y,x,k)*temp;
                Vy(y,x) = Vybar-Ey(y,x,k)*temp;
            end
        end
    end
    
end

Advertisements

41 Responses to Optical flow in matlab

  1. matt says:

    Thanks for the code!

    can you please explain what you ment by “images”?
    if i have 2 images that i want to check the flow between them how do i create one “images” var from them?

    cheers!
    Matt

  2. aruna mastani says:

    i want the program for my research

  3. chi3x10 says:

    Matt,
    Yes, create an 3-D array in the format images(imageHeight, imageWidth, numberOfFrames)

  4. Bin Liu says:

    I apply this implementation to the images sequence “office” from http://of-eval.sourceforge.net/.
    The results are quite inconsistent. I am wondering if I am using this code correctly.
    The code I wrote:
    %read the image data
    X(:,:,1)=imread(‘anim.20.pgm’); % X data type unit8
    X(:,:,2)=imread(‘anim.21.pgm’);
    %using your implementation
    [vxx,vyy]=OpticalFlow(Y,0.1,5);

    Is there any special requirements about the input data X?

    • Vt says:

      Hey Bin Liu, this is a detailed implementation of Horn paper, it’s really nice. Have you been able to fix inconsistent results with the one appearing on http://of-eval.sourceforge.net/?

      Jason thanks for the shared knowledge; did you able to get expected results on synthetic data?

  5. matt says:

    Hey,
    I’m sorry, i’m pretty new to Matlab and i still don’t understand how to use your algorithm.
    can you please explain with the next example? :
    lets say i have 2 RGB images that i change to grayscale and get the matrixes ‘E1’ and ‘E2’ with size = [240 352].
    how do i create the var ‘images’ from the 2 matrixes i’ve got?
    and then, do you think that the parameters: the created ‘images’,’alpha’ = 1, and ‘iterations’ = 16, are reasonable?

    Thank you so much!
    and i apologize for driving you nuts here 🙂
    Matt

  6. matt says:

    Bin Liu answered my question.
    i will check if i have the same problem as he stated above and let you know.

    thanks 🙂
    Matt

  7. arpit says:

    Hi

    m doing a project on Video Segmentation [ Intelligent Video Monitoring System ] n currently working on optical flow.

    Plz tell how to get SEQUENCE OF IMAGES..??

    Ben, m new to Matlab.. can u Plz tell what is Y in ur code snippet.?

  8. blanka says:

    For the image issue asked before: use the follwing snippet.

    cat(3,Image0,Image1);

  9. Berkan says:

    Hi, I believe in order for me to visualize the results I have two options:

    1-Apply a velocity threshold and convert Vx and Vy into binary values and then use regionprops to apply a bounding box around the objects that are being tracked.
    2-I can plot the optical flow lines on the (I dont know ho exactly) original image.

    Has anyone tried this algorithm and visualized the results? I am not quite sure where to feed the Vx and Vy to be able to detect moving objects.

    thanks,
    berkan

  10. samira says:

    Dear Sir,

    Thank you very much for your code. I have a single motion blurred image. Is it possible to use this code to know the optical flow of a single motion blurred image? How can I do that based on your code since your code is intended for sequence of images?

    I would like to thank for your kind help.

  11. chi3x10 says:

    Berkan,

    You can use matlab function “quiver” to plot vectors.

  12. chi3x10 says:

    Samira,

    My code is intended for sequence of images. By the way, how do you get the motion vectors from a single image?

  13. Berkan says:

    Hello Samira,

    What you are looking for is not optical flow. In order to deal with blurring (deblur), you can try “deconvblind” function at the VIP toolbox.

    hth.

  14. Berkan says:

    OK, I have tried this algortihm on a 31-frame intensity video (height*width*31).
    One interesting thing is that the output Vx and Vy matrices are of form (height*width), whereas I expected them to be (height*width*30), i.e. giving a different set of Vx-Vy vector for each pixel, comparing each consecutive frame.
    You have to put the code in a for loop if you want to achieve a frame by frame optical flow output.
    Still, the results (at least for the intentions of object tracking) are not very impressive I must say…

    • Diana says:

      Berkan,

      I noticed the same thing, expected a 30 frame output, but it’s just 1 frame. From the code, it seems it’s only returning the optic flow for the last 2 frames, even though it’s calctulating all of it. So I added a ‘k’ in the Vx, Vy assignment at the end of the for loop. What do you think?

      Vx(y,x,k) = Vxbar-Ex(y,x,k)*temp;
      Vy(y,x,k) = Vybar-Ey(y,x,k)*temp;

      D.

    • chi3x10 says:

      The algorithm is to calculate the overall optical flow of all the frames. To get the result you expected, you will have to use two consecutive frames each time and combine the results. Or you can try Diana’s method.

      Hope this helps.

  15. Saurabh says:

    hey,i did’t understand how can i get input arguement images ….m new to matlab

  16. Saurabh says:

    how to plot final result & how to choose alpha & iterations also use this [height,width,frames]=size(images); statement

  17. Dadou says:

    how can i view the optical flow .I try :
    quiver(Vx,Vy) but no resultat
    please help me .Thank u a lot.

  18. Anonymous says:

    %% to run the code above

    a1=imread(‘ur image 1’);
    a2=imread(‘ur image 2’);

    ai1 = rgb2gray(a1);
    ai2 = rgb2gray(a2);

    X(:,:,1)=ai1;
    X(:,:,2)=ai2;

    [Vx,Vy] = OpticalFlow1(X,0.1,5);

    %% to view result
    quiver(Vx,Vy);

    %% for Avi to frames or seq …

    a=aviread(‘ur video.avi’);

    s=size(a);

    for u=s(1):s(2)

    f=a(u).cdata;
    ai = rgb2gray(f);
    X(:,:,u)=ai;

    end

    hope it will be helpful !!!

  19. sathya says:

    hai guys u r doing a great work i’m not new to matlab but i’m unable to run the code which you have provided.
    i need to know how to run it at same time is any image input need to be given and can any 1 specify the software that are used to test optical flow and on wat platform, requirements etc.

    dont consider it as spam plzzzzzzz

    i desperately need it

  20. Sully says:

    Hello! I used this implementation in an application regarding obstacle avoidance. I want to ask you what is the measurement unit for the Vx and Vy components. I need this information to obtain the unit for the time-to-contact. Thanks!

  21. ami voot says:

    hi all ,I need the code for calculating average angular error for optical flow measurement .Please help me and it’s very urgent .I’m looking forward for it.

  22. sean says:

    Hi I’m new to matlab and I’m trying calculate the velocity of a falling object.
    Why do I get this error when I use optical flow?

    ??? Undefined function or method ‘OpticalFlow’ for input arguments of
    type ‘uint8’.

    Please help me

  23. Eng.Dodoo says:

    i’m sorry, but this Optical flow is not robust…

  24. usman says:

    those referring to the sourceforge website and finding inconsistencies is because you need to smooth the image first. I used smoothImg() function from Horn Schunk code at mathworks for smoothening before applying this code and I got the results shown at the site. The direction of optical flow is opposite to what is shown though.

  25. Eng.Dodoo says:

    Ok i will try it again, Thanks alot

  26. monika says:

    what do you mean by term images? if we have a video and we have to compute optical flow between the consecutive frames of that video then how should we input it and then what is the meaning of the argument images which you hav specified in your code.

  27. Moumita Das says:

    im1 = double(imread(‘C:\Users\User\Desktop\Measure_Optical Flow\1.bmp’))/256;
    im2 = double(imread(‘C:\Users\User\Desktop\Measure_Optical Flow\2.bmp’))/256;
    figure; imshow(im1);
    figure; imshow(im2);

    a1=imread(‘C:\Users\User\Desktop\Measure_Optical Flow\1.bmp’);
    a2=imread(‘C:\Users\User\Desktop\Measure_Optical Flow\2.bmp’);

    ai1 = rgb2gray(a1);
    ai2 = rgb2gray(a2);

    X(:,:,1)=ai1;
    X(:,:,2)=ai2;

    [Vx,Vy] = OpticalFlow1(X,0.1,5);

    %% to view result
    quiver(Vx,Vy);

    %% for Avi to frames or seq …

    %a=aviread();

    s=size(a);

    for u=s(1):s(2)

    f=a(u).cdata;
    ai = rgb2gray(f);
    X(:,:,u)=ai;

    end
    function[Vx,Vy] = OpticalFlow(images,alpha,iterations)
    %// Calculating optical flow of a sequence of images.
    %// images : 3D array that contains a sequence of images. size of images is (imageHeight, imageWidth, frame number)
    %// alpha
    %// iterations.
    [height,width,frames]=size(images);

    %//initialzation of u and v
    Vx = zeros(height,width);
    Vy = zeros(height,width);

    for k = 1:frames-1

    % //initialization of Ex Ey and Et
    Ex = zeros(height-1,width-1,frames-1);
    Ey = zeros(height-1,width-1,frames-1);
    Et = zeros(height-1,width-1,frames-1);

    %//calculating Ex Ey and Et in frame k.
    for x = 2:width-1
    for y = 2:height-1
    Ex(y,x,k) = (images(y+1,x+1,k)-images(y+1,x,k)+images(y,x+1,k)…
    -images(y,x,k)+images(y+1,x+1,k+1)-images(y+1,x,k+1)…
    +images(y,x+1,k+1)-images(y,x,k+1))/4;

    Ey(y,x,k) = (images(y,x,k)-images(y+1,x,k)+images(y,x+1,k)…
    -images(y+1,x+1,k)+images(y,x,k+1)-images(y+1,x,k+1)…
    +images(y,x+1,k+1)-images(y+1,x+1,k+1))/4;

    Et(y,x,k) = (images(y+1,x,k+1)-images(y+1,x,k)+images(y,x,k+1)…
    -images(y,x,k)+images(y+1,x+1,k+1)-images(y+1,x+1,k)…
    +images(y,x+1,k+1)-images(y,x+1,k))/4;
    end
    end

    for nn = 1:iterations
    for x = 2:width-1
    for y = 2:height-1

    Vxbar = (Vx(y-1,x)+Vx(y,x+1)+Vx(y+1,x)+Vx(y,x-1))/6+…
    (Vx(y-1,x-1)+Vx(y-1,x+1)+Vx(y+1,x+1)+Vx(y+1,x-1))/12;

    Vybar = (Vy(y-1,x)+Vy(y,x+1)+Vy(y+1,x)+Vy(y,x-1))/6+…
    (Vy(y-1,x-1)+Vy(y-1,x+1)+Vy(y+1,x+1)+Vy(y+1,x-1))/12;

    %// chapter 12 of Horn’s paper
    temp = (Ex(y,x,k)*Vxbar+Ey(y,x,k)*Vybar+Et(y,x,k))/(alpha^2 + Ex(y,x,k)^2 + Ey(y,x,k)^2);
    %// update u and v
    Vx(y,x) = Vxbar-Ex(y,x,k)*temp;
    Vy(y,x) = Vybar-Ey(y,x,k)*temp;
    end
    end
    end

    end

    Hi I am not able to run this your code .. there have some error in the function .Can u plz help me out

  28. Moumita Das says:

    HI,i did’t understand how can i get input argument images …. I am new to mat lab… Plz help me…

  29. I am getting NaN values for Vx and Vy. What should i do,i want to identify stationary and moving pixels in the images. Please help….

  30. pavan says:

    Try this code it works…..
    kindly zoom the figure window to observe the optical flow….

    a1=imread(‘TW1.jpg’);
    a2=imread(‘TW2.jpg’);

    b1 = rgb2gray(a1);
    b2 = rgb2gray(a2);

    images(:,:,1)=b1;
    images(:,:,2)=b2;

    [height,width,frames]=size(images);

    %//initialzation of u and v
    Vx = zeros(height,width);
    Vy = zeros(height,width);

    for k = 1:frames-1

    % //initialization of Ex Ey and Et
    Ex = zeros(height-1,width-1,frames-1);
    Ey = zeros(height-1,width-1,frames-1);
    Et = zeros(height-1,width-1,frames-1);

    %//calculating Ex Ey and Et in frame k.
    for x = 2:width-1
    for y = 2:height-1
    Ex(y,x,k) = (images(y+1,x+1,k)-images(y+1,x,k)+images(y,x+1,k)-images(y,x,k)+images(y+1,x+1,k+1)-images(y+1,x,k+1)+images(y,x+1,k+1)-images(y,x,k+1))/4;

    Ey(y,x,k) = (images(y,x,k)-images(y+1,x,k)+images(y,x+1,k)-images(y+1,x+1,k)+images(y,x,k+1)-images(y+1,x,k+1)+images(y,x+1,k+1)-images(y+1,x+1,k+1))/4;

    Et(y,x,k) = (images(y+1,x,k+1)-images(y+1,x,k)+images(y,x,k+1)-images(y,x,k)+images(y+1,x+1,k+1)-images(y+1,x+1,k)+images(y,x+1,k+1)-images(y,x+1,k))/4;
    end
    end

    for nn = 1:5
    for x = 2:width-1
    for y = 2:height-1

    Vxbar = (Vx(y-1,x)+Vx(y,x+1)+Vx(y+1,x)+Vx(y,x-1))/6+(Vx(y-1,x-1)+Vx(y-1,x+1)+Vx(y+1,x+1)+Vx(y+1,x-1))/12;

    Vybar = (Vy(y-1,x)+Vy(y,x+1)+Vy(y+1,x)+Vy(y,x-1))/6+(Vy(y-1,x-1)+Vy(y-1,x+1)+Vy(y+1,x+1)+Vy(y+1,x-1))/12;

    %// chapter 12 of Horn’s paper
    temp = (Ex(y,x,k)*Vxbar+Ey(y,x,k)*Vybar+Et(y,x,k))/(0.5^2 + Ex(y,x,k)^2 + Ey(y,x,k)^2);
    %update u and v
    Vx(y,x) = Vxbar-Ex(y,x,k)*temp;
    Vy(y,x) = Vybar-Ey(y,x,k)*temp;
    end
    end
    end

    end

    quiver(Vx,Vy);

  31. heloo guys..i am doing some work to detect motion of a vehicle in opposite direction using matlab..can any body help me regarding this!

  32. Thanks extremely helpful. Will certainly share site with my pals

  33. With havin so much content and articles do you ever run into
    any problems of plagorism or copyright infringement? My site
    has a lot of completely unique content I’ve either authored myself
    or outsourced but it looks like a lot of it is popping it up all over the web without
    my permission. Do you know any methods to help stop content from being ripped off?
    I’d really appreciate it.

  34. Usually I do not read article on blogs, but I would like to say that this write-up very pressured me to take a look at and do it! Your writing taste has been surprised me. Thanks, quite nice post.|

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: