Video stabilization using Lucas-Kanade method

Question

I got an assignment in a video processing course - to stabilize a video using the Lucas-Kanade method. Assuming the matlab code I wrote for performing LK on 2 images works (i.e, I can get the optical flow matrices u, v that allow me to warp one image towards another) - how can I use it to perform a video stabilization?

As far as I understand, the main idea is to calculate the OF for each 2 adjacent frames, and to warp all the images towards the first one. that means that I need to keep the total warp that I accumulate during the whole process- and practically: saving total warp matrices u_tot, v_tot and add to them the u, v matrices I got from 2 adjacent frames, to perform a warp of the latter frame according to u_tot, v_tot and so on (notice the matlab code below).

Am I going in the right direction? I'd like to get an advise... thanks

---

matlab code:

EDIT: changes the code + added the WarpImage imp. (23.4.17)

function StabilizedVid =LucasKanadeVideoStabilization(InputVid,WindowSize,MaxIter,NumLevels,NumOfFrames)
    % InputVid is a VideoFileReader object
    % WindowSize - window size for the LK phase
    % MaxIter - maximal number of iterations for the LK phase
    % NumLevels - number of levels in the pyramid for the LK phase
    % NumOfFrames - number of frames to process

    reset(InputVid);
    StabilizedVid = vision.VideoFileWriter('StabilizedVid.avi');
    cnt = 1;
    if (~isDone(InputVid))
        frame1 = im2double(rgb2gray(step(InputVid)));
        first_frame = true;
    else
        error('no frames in the video');
    end

    while ~isDone(InputVid)
        cnt = cnt+1;
        frame2 = im2double(rgb2gray(step(InputVid)));
        [frame1, frame2] = sizeCheck(frame1, frame2, NumLevels);
        if (first_frame)
            u_tot = zeros(size(frame1));
            v_tot = zeros(size(frame1));
            step(StabilizedVid,frame1);
            first_frame = false;
        end  
        disp(['processing frames no. ', num2str(cnt-1), ', ', num2str(cnt)]);
        tic
        [u ,v] = LucasKanadeOpticalFlow(frame1,frame2,WindowSize,MaxIter,NumLevels); 
        disp(['time: ', num2str(toc)]);
        u_tot = u_tot + u;
        v_tot = v_tot + v;
        frame2_warp = WarpImage(frame2, u_tot, v_tot);
        step(StabilizedVid,frame2_warp);
        if cnt >= NumOfFrames
            break
        end
        frame1 = frame2_warp;
    end
    release(StabilizedVid);
end

%%%%%%%%%%%%%%%%%%%%

function I_warp = WarpImage(I,u,v)
    % I - image to warp
    % u,v - the optical flow parameters gotten from the LK algorithm        

    [x, y] = meshgrid(1:size(I,2),1:size(I,1));
    % interpulate
    I_warp = interp2(I, x+u, y+v, 'cubic');
    % in case of NaN values - put the original values of I instead
    I_warp(isnan(I_warp)) = I(isnan(I_warp));

end

Answer 1

For video stabilization, you need to estimate motion vector with respect to a reference frame $F_{reference}$ (which let's assume it to be the first frame). Then your motion model should be determined, for simplicity let's assume it to be Translational motion. Therefore, we need to estimate motion in $x$ direction and motion in $y$ direction for each frame. Now you have two options to estimate the change. If you are interested to estimate it through optical flow you can do either these:

1. Estimate regular conventional optical flow from frame to frame
2. Estimate Optical flow of first frame to the $i$ $th$ frame.

I accumulate during the whole process- and practically...

If you've chose the first, then Yes, you need to accumulate the optical flow, but I think it would be simpler to go for the second method.

As a hint, I suggest track a token or compute optical flow of a small patch of the image.

Edit:

I feel the problem is with your Warp function, the code below is part of video stabilization assignment in my class, try something like this sample. The translational parameters are estimated using token based tracking.

function [ ou_I ] = frame_translation_func( in_I,Tx,Ty,ou_RoWs,ou_CuLoMns,ou_layer)
ou_I = uint8(zeros(ou_RoWs,ou_CuLoMns,ou_layer));
for i2=1:ou_RoWs
    for j2=1:ou_CuLoMns
        if (i2 - Tx<1) || (j2 - Ty<1) || (i2 - Tx>ou_RoWs) || (j2 - Ty>ou_CuLoMns)
            ou_I(i2,j2,:)=0;
        else
            ou_I(i2,j2,:) = in_I(i2 - Tx,j2 - Ty,:);
        end
    end
end
end