Time for action – implementing a spatiotemporal averaging filter

Now that we know the steps of the process, we can write a function that performs it. The function will have some similarities to the spatiotemporal deinterlacing function implemented in the previous example. It will take as input a grayscale video in a matrix variable (if it is not grayscale it will convert each frame to grayscale) and the filtering neighborhood, and will produce as output a matrix containing the filtered video:

function [vid] = SpatioTemporalAveraging(vid,nhood)
 
% Function for de-interlacing a video using Field Merging
% Inputs:
%        vid     – Input video matrix (we assume grayscale video) 
%        nhood   - Dimensions of filtering neighborhood
% Output:   
%        vidOut     - Output video matrix (filtered)
 
vid = single(vid);  % Convert matrix to single to perform averaging
  vidOut = zeros(size(vid)); % Create a zero-valued matrix to hold the result
  % For all frames except the ones that fall off limits
  for frame = ceil(nhood(3)/2):size(vid,3)-ceil(nhood(3)/2)
    % For all rows except the ones that fall off limits
    for row = ceil(nhood(1)/2):size(vid,1)-ceil(nhood(1)/2)
      % For all columns except the ones tha fall off limits
      for column = ceil(nhood(2)/2):size(vid,2)-ceil(nhood(2)/2)
        % Crop the neighborhood area           
        neighborhood = vid(row-round(nhood(1)/2)+1:row+round(nhood(1)/2),...
        column-round(nhood(2)/2)+1:column+round(nhood(2)/2), 
        frame-round(nhood(3)/2)+1:frame+round(nhood(3)/2));
        % Calculate its mean and assign it to the central pixel    
        vidOut(row,column,frame) = mean(neighborhood(:));
      end
    end
  end
vidOut = uint8(vidOut);   % Convert matrix to type uint8

Now, it is time to test our function. Beware that this is a time-consuming function, because of the three nested loops. It's going to take a while for each frame, so don't use it on a very large video:

  1. We'll use our driving video comprising 28 frames:
    >> obj = VideoReader('inter.avi'),
    >> vid = read(obj);
  2. Now, we will convert our color video matrix to grayscale, frame-by-frame, right after we initialize a new uint8 type matrix with zeros, to store our result in:
    >> grayVid = uint8(zeros(size(vid,1), size(vid,2), size(vid,4)));
    >> for i = 1:size(vid,4),
    grayVid(:,:,i) = rgb2gray(vid(:,:,:,i));
    end
  3. Let's test our function now, for a cubical neighborhood of 3 rows, 3 columns, and 3 frames:
    >> filteredVid = SpatioTemporalAveraging(grayVid,[3 3 3]);
  4. After a while (the time needed depends on your CPU), you will see in your workspace a variable called filteredVid. Let's display a frame of this matrix next to the same frame of our original matrix, to see what the results of the filtering process were:
    >> subplot(1,2,1),imshow(grayVid(:,:,10)),title('Original frame')
    >> subplot(1,2,2),imshow(filteredVid(:,:,10)),title('Filtered frame')
    Time for action – implementing a spatiotemporal averaging filter

What just happened?

This example has demonstrated how you can program a spatiotemporal blurring process for grayscale video streams in MATLAB. The function began with a type conversion of the input video matrix, so that averaging was feasible followed by the initialization of a new matrix that stored the resulting video.

Then, we moved on to three nested for loops, by making sure to set the limits of each of the loops in a manner that ensured our process won't fall off the borders of the image. This choice resulted in an image with a small black border in all its four edges. The size of the border will depend on the size of the filtering neighborhood.

In each step of the triple for loop, we cropped the cubical neighborhood centered at the pixel under examination (row, column, and frame):

neighborhood = vid(row-round(nhood(1)/2)+1:row+round(nhood(1)/2),...column-round(nhood(2)/2)+1:column+round(nhood(2)/2), ...frame-round(nhood(3)/2)+1:frame+round(nhood(3)/2));

Then, the average value of the neighborhood was calculated and used to replace the value of the central pixel of the neighborhood:

vidOut(row,column,frame) = mean(neighborhood(:));

The steps 1 through 4 evaluated the usage of the filtering process. First the video was loaded and stored to a matrix variable. Then, in step 2, the input video was converted to grayscale, frame-by-frame. Each of the frames is stored in a type uint8 matrix initialized with zero values. The step 3 applied the filter to our grayscale video and finally, in step 4 we displayed one frame from the resulting video next to the same frame from the original video. As you can see, the result of such a spatiotemporal filtering process was adding a shaky effect to our video. The more the motion in the area of the frame was, the more intense the effect will be.

Have a go hero – creating a spatiotemporal median filter

Before we move on, you should try to alter the filtering function we made in the previous example, to give a spatiotemporal median filtering result instead of the averaging result we have programmed. The only change will be in the line where the averaging calculation was performed. If you want, you can also add a third input to the previous function to let the users choose which of the two methods they prefer to use. A third addition to the code would be to add padding, to eliminate the thin black border from the result.

Using convolution for spatiotemporal averaging

The previous implementation of the spatiotemporal averaging filter was interesting, but had a real big disadvantage; it was really slow. The alternative way to perform such a process is to take advantage of the multidimensional version of convolution. In Chapter 5, 2-Dimensional Image Filtering, we discussed the method of convolution and described how it can be used for image filtering. Now, it is time to expand the convolution operation to more dimensions and perform spatiotemporal (three-dimensional) filtering to our video frames.

In order to implement the process of spatiotemporal averaging by means of convolution, we have to use the function convn provided by MATLAB. This function is called exactly like conv2, with the only difference that the inputs must be n-dimensional (in our case three-dimensional). Let's see how this works.

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.
Reset