Wednesday, December 7, 2016

Activity 11 - Basic Video Processing

After a TIRING semester we're finally at the last activity! It's time to delve into video processing! Hopefully, we'll be able to apply some of the things we've learned throughout this long and arduous image processing journey in some basic video processing! So for the last time (this semester), let's dive in!

BASIC VIDEO PROCESSING

The goal of this final blog post is to use video processing to obtain measure some kinematic quantity. A video is basically a set of images captured and played back fast enough such that the impression of movement is perceived. The difference, however, lies in the fact that there is now a new dimension to consider: time. Because of this, it is instructive to introduce the concepts of frame rate and the Nyquist theorem. The Nyquist sampling states that,

"The sampling frequency should be at least twice the highest frequency contained in the signal." [1]

In this case, our signal is position of some moving object in the video and  our sampling frequency corresponds to the frame rate of our video. For the purposes of this demonstration, we utilize a setup wherein the object of interest has a position that changes periodically and thus the frequency of its movement is fairly constant. Thus, if the object is moving at a frequency of let's say 10 Hz, then the frame rate at which we capture the video should be AT LEAST 30 Hz. Now that's handled, let's go straight into the process. 

I chose to analyze an ideal simple pendulum system and aim to calculate the acceleration due to gravity g. The setup is shown in Figure 1. 

Figure 1. Simple pendulum setup
 As shown above, I constructed a simple pendulum using a thin nylon string (32cm) glued to a tennis ball. This was to ensure that the setup was consistent with the  "massless" and "inextensible" assumptions of the string in a simple pendulum. The pendulum was then suspended over a lamp arm and allowed to swing freely. A red background was placed behind the pendulum to aid in the segmentation of the tennis ball later on. Lastly, a ruler was placed in the vicinity of the scene in order to provide a scale factor for the conversion between pixels and centimeters. Consistent with the small angle approximation, the pendulum was displaced from equilibrium by a horizontal distance of ~4cm (around 1.4 degrees angular displacement) and let go. The oscillation of the pendulum was captured using a Lumix GX1 camera. Three five second trial vidoes were taken to serve as the data set from which g was to be calculated. The videos were taken at 29 fps. The image sequences were extracted from the videos using FFMPEG software with the time interval between frames thus being 1/29=34.48 ms. Using an fps of 29, we note from the Nyquist theorem that the highest frequency that we can measure (before aliasing occurs) is around 14.5 Hz. It is unlikely that the pendulum will oscillate at such a fast rate and thus, we're safe!

Because Scilab was being uncooperative again, I decided to use MATLAB once again! The first step I made was to crop all of the images equally in order to reduce file sizes and help out the segmentation. I did this with the code in Figure 2.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
for j=1:148
    ReadName=['C:\Users\domini8888\Desktop\Activity 11\pics\second\images-',num2str(j),'.jpg'];
    image=imread(ReadName);
    
    croplength=350;
    cropped=image(:,croplength:size(image,2)-croplength,:);
    
    basename='C:\Users\domini8888\Desktop\Activity 11\pics\second_cropped\images-';
    FileName=[basename,num2str(j),'.jpg'];
    imwrite(cropped,FileName);
end
Figure 2. Code used to "batch crop" all the images equally


Just for reference, I created a gif (with GIMP!) using one set of the images. This is shown in Figure 3. Note that I had to drastically reduce the size of the gif just so it wouldn't take too long to load. The time between frames, however, was set to 34.48 ms so what you see in the gif is practically what was filmed.
Figure 3. "Video" taken of the ball in GIF form
Now on to the processing! During the filming process, the scene lighting was made to be consistent. Moreover, it was ensured that the ball was the only object moving in the background. Thus, these two constant factors allow us to apply the same processing sequence on every image in the image sequence taken.

Thus, let's look at one individual photo and try to process it.

Figure 4. A "snapshot" from one of the videos taken
Just for reference, I measured the diameter of the tennis ball using MS Paint to be around 104 pixels. We'll use this value later!

 Now, we must find some way of segmenting the tennis ball out from the background. Of course, the first step would be to see if we can segment the image simply by grayscale thresholding. This requires us to look at the gray scale histogram. 

Figure 5. Gray scaled version of Figure 4 and it's resulting histogram
Taking a look at the histogram in Figure 5, we see that it is highly unlikely that we'll be able to segment the tennis ball using gray scale thresholding cleanly. Thus, when gray scale thresholding fails, it's time for color segmentation (see my blog on color segmentation if you want a more thorough discussion)! Luckily, I made sure that the chromaticities contained in the tennis ball were distinct within the image. Consider the ROI used in Figure 6.

Figure 6. ROI taken from the snapshot in Figure 4. 

To avoid oversegmentation (and be on the safe side), I chose to use parametric color segmentation. The code used to implement this is shown in Figure 7.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
image=imread('C:\Users\domini8888\Desktop\Activity 11\pics\first_cropped\images-1.jpg');
roi=imread('C:\Users\domini8888\Desktop\Activity 11\pics\roi.jpg');
roi=double(roi);
image=double(image);

%convert to NCC
I=roi(:,:,1)+roi(:,:,2)+roi(:,:,3);
I(find(I==0))=1000000;
r=roi(:,:,1)./I;
g=roi(:,:,2)./I;

I_image=image(:,:,1)+image(:,:,2)+image(:,:,3);
I_image(find(I==0))=1000000;
r_image=image(:,:,1)./I_image;
g_image=image(:,:,2)./I_image;

p_r=(1 ./(std(r(:))*sqrt(2*pi)))*exp(-((r_image-mean(r(:))) ./ (2*std(r(:)))).^2);
p_g=(1 ./(std(g(:))*sqrt(2*pi)))*exp(-((g_image-mean(g(:))) ./ (2*std(g(:)))).^2);

jointrg=p_r.*p_g;
jointrg=jointrg./max(jointrg(:));
jointrg=mat2gray(jointrg);
%imshow((jointrg));
imwrite(jointrg,'C:\Users\domini8888\Desktop\Activity 11\pics\jointrg.jpg')
Figure 7. Code used to parametrically color segment the tennis ball

Now let's take a look at the output!

Figure 8. Color segmented tennis ball
Pretty good if I must say so! The only problem is that we need to find some way of combining all the regions of the tennis ball into one big blob. We can now go back to activity 8 and apply our morphological cleaning skills! Check out the code used to do this.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
%morphological cleaning + thresholding + cleaning again
SE1=strel('disk',25);
SE2=strel('disk',1);
SE3=strel('disk',9);
edited=imopen(jointrg,SE2);
edited=imclose(edited,SE1);
edited=edited>0.0393;
edited=imopen(edited,SE3);
%imshow(edited);
imwrite(edited,'C:\Users\domini8888\Desktop\Activity 11\pics\test.jpg')
Figure 9. Code used to clean up the image in Figure 8

As shown in the code in Figure 9, an opening operation is first performed with a circular structuring element to eliminate any small indiscernable blobs not part of the tennis ball that are still present in the image. A radius of 25 is well below the radius of the tennis ball (about 52 from what was measured using MS Paint) and thus won't affect the tennis ball much. We then apply a closing operation to fill in the gaps within the region of the tennis ball. The result is then thresholded to binarize the grayscale image and then an opening operation is performed again to clean up artifacts. Figure 9 shows the development of the image as it went through the process. This hopefully helps explain the motivation behind the process implemented.

Figure 10. Development of the image as it went through the process. The variable r indiciates the radius of the circular structuring element used and Thresh indicates the thresholding value used.


Measuring the diameter of the last image in Figure 10, we see that the diameter is abour 94-108 pixels, giving roughly a mean value of 101. This is pretty close to the measured value earlier of 104! Thus, the process for this single image was successful. All that's left is to apply it to all the other images! Shown in Figure 11is the GIF of the segmented tennis ball using the same video used to generate Figure 3.

Figure 11. "Original" video (left) and the segmented video (right)

Now all that's left is to track the movement of the ball. We do this by finding the centroid of the white blob corresponding to the tennis ball. To do this, we use  [x,y]=find(image>0) to find the x and y indices of the white pixels per image. Then, taking the mean of these x and y values, we end up with an estimate for the x and y points for the centroid of the image. Because we're using the small angle approximation, we only really need to look at the x values and plot these against time. Remember that each frame is separated by a time interval of 1/29=34.48ms. The code used to obtain these values per image and plot the resulting oscillation along the x-direction is shown in Figure 12. 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
x_points=[];
y_points=[];
for j=1:150
num=num2str(j);
ReadName=['C:\Users\domini8888\Desktop\Activity 11\pics\first_seg\images-',num2str(j),'.jpg'];
image=imread(ReadName);
[x,y]=find(image>0);
cent_x=mean(x(:));
cent_y=mean(y(:));
x_points=cat(1,y_points,cent_x);
y_points=cat(1,y_points,cent_y);

end

fps=29;
t_int=1./29;
t=0:t_int:(size(x_points,1)-1)*t_int;
plot(t,x_points);

Figure 12. Code used to find the centroid of the tennis ball within each image and then plot the oscillation along x with respect to time.


For the three videos taken, the plots are shown below.
Figure 13. Oscillation along the x-axis of the pendulum bob as a function of time as measured from the three videos taken
From this data, we can find the period of oscillation by finding the time differences between the occurences of the maxima. This is done using the code in Figure 14.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
[pks,locs]=findpeaks(x_points);

pkloc=[];
for i=1:size(pks,1)
    if pks(i)>300
        pkloc=cat(1,pkloc,locs(i));
    end
end

disp(pkloc);

period=[];
for i=1:size(pkloc)-1
    diff=t(pkloc(i+1))-t(pkloc(i));
    period=cat(1,period,diff);
end

disp(mean(period));
Figure 14. Code used to calculate the period of oscillation

Thus, for the three cases, we obtain periods of 1.2184 s, 1.2184s, and 0.9138s. This gives us a mean of 1.1169s and a standard deviation of 0.1759s. Using the formula for the period of a simple pendulum [2], and ImageJ to get a string length of 32.3cm (scale value of 12.66pixels/cm using the ruler in the image), we can calculate the g to be L/(T/2pi)^2=10.2215 m/s^2. Comparing this with the textbook value of 9.81m/s^2, we obtain a percent deviation of 4.19%! Pretty good!!!!


I really feel I put in a lot of effort into this activity. It was pretty summative of a lot of things and lessons I've learned throughout the semester in Applied Physics 186. I also believe that I went above and beyond in my explanations of the steps made. For all of this, I believe I deserve a 12/10 for my effort.

Acknowledgements:
I'd like to acknowledge Roland Romero for helping me generate GIF files using GIMP, and pointing out to me how simple the centroid calculation could be.

References:
[1](n.d.). Retrieved from http://redwood.berkeley.edu/bruno/npb261/aliasing.pdf
[2] The Simple Pendulum - Boundless Open Textbook. (n.d.). Retrieved December 07, 2016, from https://www.boundless.com/physics/textbooks/boundless-physics-textbook/waves-and-vibrations-15/periodic-motion-123/the-simple-pendulum-431-8324/
Posted by at

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)