Motion Tracking/Detection in MATLAB using Kalman Filter
MATAB has great capabilities to process video & pictures. One of these capabilities we have tested is the motion detection in a pre-recorded video, using Kalman Filtering technique.
Here is an example video of Motion Tracking using Kalman Filter:
Here is an example video of Motion Tracking using Kalman Filter:
Video for illustration Purpose Only
Basic background of Kalman Filter:
The Kalman filter, also known as linear quadratic estimation (LQE), is an algorithm that uses a series of measurements observed over time, containing noise
(random variations) and other inaccuracies, and produces estimates of
unknown variables that tend to be more precise than those based on a
single measurement alone. More formally, the Kalman filter operates recursively on streams of noisy input data to produce a statistically optimal estimate of the underlying system state. The filter is named for Rudolf (Rudy) E. Kálmán, one of the primary developers of its theory.
The Kalman filter has numerous applications, e.g, for guidance, navigation & control of vehicles, particularly aircraft and spacecraft.
 |
| Kalman Filter Algorithm (courtsey: colorado.edu) |
MATLAB CODE:
%See the below code update of 11/11/2014 for implementing the code through VideoReader instead of aviread. Now with this update you can solve many of your problems of using your own video for this code.
%See the below code update of 11/11/2014 for implementing the code through VideoReader instead of aviread. Now with this update you can solve many of your problems of using your own video for this code.
%Code in the Bold, & comments in normal font
clear all;
close all;
clc
%% Read video into MATLAB using aviread
video = aviread('rhinos.AVI');
%'n' for calculating the number of frames in the video file
n = length(video);
% Calculate the background image by averaging the first 10 images
temp = zeros(size(video(1).cdata));
[M,N] = size(temp(:,:,1));
for i = 1:10
temp = double(video(i).cdata) + temp;
end
imbkg = temp/10;
% Initialization step for Kalman Filter
centroidx = zeros(n,1);
centroidy = zeros(n,1);
predicted = zeros(n,4);
actual = zeros(n,4);
% % Initialize the Kalman filter parameters
% R - measurement noise,
% H - transform from measure to state
% Q - system noise,
% P - the status covarince matrix
% A - state transform matrix
R=[[0.2845,0.0045]',[0.0045,0.0455]'];
H=[[1,0]',[0,1]',[0,0]',[0,0]'];
Q=0.01*eye(4);
P = 100*eye(4);
dt=1;
A=[[1,0,0,0]',[0,1,0,0]',[dt,0,1,0]',[0,dt,0,1]'];
% loop over all image frames in the video
kfinit = 0;
th = 38;
for i=1:n
imshow(video(i).cdata);
hold on
imcurrent = double(video(i).cdata);
% Calculate the difference image to extract pixels with more than 40(threshold) change
diffimg = zeros(M,N);
diffimg = (abs(imcurrent(:,:,1)-imbkg(:,:,1))>th) ...
| (abs(imcurrent(:,:,2)-imbkg(:,:,2))>th) ...
| (abs(imcurrent(:,:,3)-imbkg(:,:,3))>th);
% Label the image and mark
labelimg = bwlabel(diffimg,4);
markimg = regionprops(labelimg,['basic']);
[MM,NN] = size(markimg);
% Do bubble sort (large to small) on regions in case there are more than 1
% The largest region is the object (1st one)
for nn = 1:MM
if markimg(nn).Area > markimg(1).Area
tmp = markimg(1);
markimg(1)= markimg(nn);
markimg(nn)= tmp;
end
end
% Get the upper-left corner, the measurement centroid and bounding window size
bb = markimg(1).BoundingBox;
xcorner = bb(1);
ycorner = bb(2);
xwidth = bb(3);
ywidth = bb(4);
cc = markimg(1).Centroid;
centroidx(i)= cc(1);
centroidy(i)= cc(2);
% Plot the rectangle of background subtraction algorithm -- blue
hold on
rectangle('Position',[xcorner ycorner xwidth ywidth],'EdgeColor','b');
hold on
plot(centroidx(i),centroidy(i), 'bx');
% Kalman window size
kalmanx = centroidx(i)- xcorner;
kalmany = centroidy(i)- ycorner;
if kfinit == 0
% Initialize the predicted centroid and volocity
predicted =[centroidx(i),centroidy(i),0,0]' ;
else
% Use the former state to predict the new centroid and volocity
predicted = A*actual(i-1,:)';
end
kfinit = 1;
Ppre = A*P*A' + Q;
K = Ppre*H'/(H*Ppre*H'+R);
actual(i,:) = (predicted + K*([centroidx(i),centroidy(i)]' - H*predicted))';
P = (eye(4)-K*H)*Ppre;
% Plot the tracking rectangle after Kalman filtering -- red
hold on
%% Read video into MATLAB using aviread
video = aviread('rhinos.AVI');
%'n' for calculating the number of frames in the video file
n = length(video);
% Calculate the background image by averaging the first 10 images
temp = zeros(size(video(1).cdata));
[M,N] = size(temp(:,:,1));
for i = 1:10
temp = double(video(i).cdata) + temp;
end
imbkg = temp/10;
% Initialization step for Kalman Filter
centroidx = zeros(n,1);
centroidy = zeros(n,1);
predicted = zeros(n,4);
actual = zeros(n,4);
% % Initialize the Kalman filter parameters
% R - measurement noise,
% H - transform from measure to state
% Q - system noise,
% P - the status covarince matrix
% A - state transform matrix
R=[[0.2845,0.0045]',[0.0045,0.0455]'];
H=[[1,0]',[0,1]',[0,0]',[0,0]'];
Q=0.01*eye(4);
P = 100*eye(4);
dt=1;
A=[[1,0,0,0]',[0,1,0,0]',[dt,0,1,0]',[0,dt,0,1]'];
% loop over all image frames in the video
kfinit = 0;
th = 38;
for i=1:n
imshow(video(i).cdata);
hold on
imcurrent = double(video(i).cdata);
% Calculate the difference image to extract pixels with more than 40(threshold) change
diffimg = zeros(M,N);
diffimg = (abs(imcurrent(:,:,1)-imbkg(:,:,1))>th) ...
| (abs(imcurrent(:,:,2)-imbkg(:,:,2))>th) ...
| (abs(imcurrent(:,:,3)-imbkg(:,:,3))>th);
% Label the image and mark
labelimg = bwlabel(diffimg,4);
markimg = regionprops(labelimg,['basic']);
[MM,NN] = size(markimg);
% Do bubble sort (large to small) on regions in case there are more than 1
% The largest region is the object (1st one)
for nn = 1:MM
if markimg(nn).Area > markimg(1).Area
tmp = markimg(1);
markimg(1)= markimg(nn);
markimg(nn)= tmp;
end
end
% Get the upper-left corner, the measurement centroid and bounding window size
bb = markimg(1).BoundingBox;
xcorner = bb(1);
ycorner = bb(2);
xwidth = bb(3);
ywidth = bb(4);
cc = markimg(1).Centroid;
centroidx(i)= cc(1);
centroidy(i)= cc(2);
% Plot the rectangle of background subtraction algorithm -- blue
hold on
rectangle('Position',[xcorner ycorner xwidth ywidth],'EdgeColor','b');
hold on
plot(centroidx(i),centroidy(i), 'bx');
% Kalman window size
kalmanx = centroidx(i)- xcorner;
kalmany = centroidy(i)- ycorner;
if kfinit == 0
% Initialize the predicted centroid and volocity
predicted =[centroidx(i),centroidy(i),0,0]' ;
else
% Use the former state to predict the new centroid and volocity
predicted = A*actual(i-1,:)';
end
kfinit = 1;
Ppre = A*P*A' + Q;
K = Ppre*H'/(H*Ppre*H'+R);
actual(i,:) = (predicted + K*([centroidx(i),centroidy(i)]' - H*predicted))';
P = (eye(4)-K*H)*Ppre;
% Plot the tracking rectangle after Kalman filtering -- red
hold on
rectangle('Position',[(actual(i,1)-kalmanx) (actual(i,2)-kalmany) xwidth ywidth], 'EdgeColor', 'r','LineWidth',1.5);
hold on
plot(actual(i,1),actual(i,2), 'rx','LineWidth',1.5);
drawnow;
end
hold on
plot(actual(i,1),actual(i,2), 'rx','LineWidth',1.5);
drawnow;
end
%end of the code
%Copyright (c) 2012, Xing
%All rights reserved.
%Copyright (c) 2012, Xing
%All rights reserved.
NOTE:
The above code sometimes generate the following error & warning, in the initial step itself (while reading the .avi file using "aviread" function.
Warning: AVIREAD will be removed in a future release. Use VIDEOREADER instead.
> In aviread at 26
In motionnew at 4
Error using aviread (line 80)
Unable to locate decompressor to decompress video stream
Error in motionnew (line 4)
video = aviread('traffic.avi');
> In aviread at 26
In motionnew at 4
Error using aviread (line 80)
Unable to locate decompressor to decompress video stream
Error in motionnew (line 4)
video = aviread('traffic.avi');
# First warning is due to the fact, that aviread function has became obsolete & MATLAB team decided to replace it with a new function called, VideoReader.
# In order to avoid the second error, you need to install the avi video codec, the best thing available in the internet is the K-Lite Mega Codec Pack.
~Team Digital iVision Lab.
# do comment for any code request, help or queries.
Original Unmodified Source can be downloaded here.
If having problem with aviread or mmread function kindly visit the following post.
How to resolve MATLAB Problem With "aviread" or "VideoReader" command, while reading the avi file?
UPDATE : 20/10/2014
Some frequent errors Faced by users, & Solutions.
(1)Cannot find an exact (case-sensitive) match for 'videoReader'
The closest match is: VideoReader in
C:\ProgramFiles\MATLAB\R2013a\toolbox\matlab\audiovideo\@VideoReader\VideoReader.m
Solution: You need to use VideoReader instead of videoReader or VideoReader.
(2)No appropriate method, property, or field cdata for class VideoReader
Solution: Since the above code is implemented for aviread function, it will not be same as for VideoReader class implementation.
(3) Warning: AVIREAD will be removed in a future release. Use VIDEOREADER instead.
(4)Error using aviread (line 140)
Solution: This error is occuring due after implementing the FFMPEG decompression. This is due to the fact that aviread supports only those files only which are 8-bit indexed image or grayscale image. In modern world generally the video are much enhanced & their frames are having more color depth, so this error occurs while using them, even after decompressing step using ffmpeg.
> In aviread at 26
Solution: The warning is due to the fact that aviread command has become obsolete & it will be removed from further version.
The AVI file must be 8-bit Indexed or grayscale images, 16-bit grayscale, or 24-bit TrueColor images
UPDATE : 11/11/2014 (Major Code Update)
Many user are facing problems with the above code, that have been written using "aviread" funciton. So I am re-editing this code to be used using VideoReader class, so that you can easily use it with your own videos.
%Code starts here
clc;
close all;
clear all;
video = VideoReader('rhinos.avi'); %in place of aviread
%nframes = length(video);
nframes=video.NumberOfFrames;
for i=1:nframes
mov(i).cdata=read(video,i)
close all;
clear all;
video = VideoReader('rhinos.avi'); %in place of aviread
%nframes = length(video);
nframes=video.NumberOfFrames;
for i=1:nframes
mov(i).cdata=read(video,i)
%creating '.cdata' field to avoid much changes to previous code
end
temp = zeros(size(mov(1).cdata));
[M,N] = size(temp(:,:,1));
for i = 1:10
temp = double(mov(i).cdata) + temp;
end
imbkg = temp/10;centroidx = zeros(nframes,1);
centroidy = zeros(nframes,1);
predicted = zeros(nframes,4);
actual = zeros(nframes,4);
R=[[0.2845,0.0045]',[0.0045,0.0455]'];
H=[[1,0]',[0,1]',[0,0]',[0,0]'];
Q=0.01*eye(4);
P = 100*eye(4);
dt=1;
A=[[1,0,0,0]',[0,1,0,0]',[dt,0,1,0]',[0,dt,0,1]'];
kfinit = 0;
th = 38;
for i=1:nframes
imshow(mov(i).cdata);
hold on
imcurrent = double(mov(i).cdata);
diffimg = zeros(M,N);
diffimg = (abs(imcurrent(:,:,1)-imbkg(:,:,1))>th) ...
| (abs(imcurrent(:,:,2)-imbkg(:,:,2))>th) ...
| (abs(imcurrent(:,:,3)-imbkg(:,:,3))>th);
labelimg = bwlabel(diffimg,4);
markimg = regionprops(labelimg,['basic']);
[MM,NN] = size(markimg);
for nn = 1:MM
if markimg(nn).Area > markimg(1).Area
tmp = markimg(1);
markimg(1)= markimg(nn);
markimg(nn)= tmp;
end
end
bb = markimg(1).BoundingBox;
xcorner = bb(1);
ycorner = bb(2);
xwidth = bb(3);
ywidth = bb(4);
cc = markimg(1).Centroid;
centroidx(i)= cc(1);
centroidy(i)= cc(2);
hold on
rectangle('Position',[xcorner ycorner xwidth ywidth],'EdgeColor','b');
hold on
plot(centroidx(i),centroidy(i), 'bx');
kalmanx = centroidx(i)- xcorner;
kalmany = centroidy(i)- ycorner;
if kfinit == 0
predicted =[centroidx(i),centroidy(i),0,0]' ;
else
predicted = A*actual(i-1,:)';
end
kfinit = 1;
Ppre = A*P*A' + Q;
K = Ppre*H'/(H*Ppre*H'+R);
actual(i,:) = (predicted + K*([centroidx(i),centroidy(i)]' - H*predicted))';
P = (eye(4)-K*H)*Ppre;
hold on
rectangle('Position',[(actual(i,1)-kalmanx)...
(actual(i,2)-kalmany) xwidth ywidth],'EdgeColor','r','LineWidth',1.5);
hold on
plot(actual(i,1),actual(i,2), 'rx','LineWidth',1.5);
drawnow;
end
end
temp = zeros(size(mov(1).cdata));
[M,N] = size(temp(:,:,1));
for i = 1:10
temp = double(mov(i).cdata) + temp;
end
imbkg = temp/10;centroidx = zeros(nframes,1);
centroidy = zeros(nframes,1);
predicted = zeros(nframes,4);
actual = zeros(nframes,4);
R=[[0.2845,0.0045]',[0.0045,0.0455]'];
H=[[1,0]',[0,1]',[0,0]',[0,0]'];
Q=0.01*eye(4);
P = 100*eye(4);
dt=1;
A=[[1,0,0,0]',[0,1,0,0]',[dt,0,1,0]',[0,dt,0,1]'];
kfinit = 0;
th = 38;
for i=1:nframes
imshow(mov(i).cdata);
hold on
imcurrent = double(mov(i).cdata);
diffimg = zeros(M,N);
diffimg = (abs(imcurrent(:,:,1)-imbkg(:,:,1))>th) ...
| (abs(imcurrent(:,:,2)-imbkg(:,:,2))>th) ...
| (abs(imcurrent(:,:,3)-imbkg(:,:,3))>th);
labelimg = bwlabel(diffimg,4);
markimg = regionprops(labelimg,['basic']);
[MM,NN] = size(markimg);
for nn = 1:MM
if markimg(nn).Area > markimg(1).Area
tmp = markimg(1);
markimg(1)= markimg(nn);
markimg(nn)= tmp;
end
end
bb = markimg(1).BoundingBox;
xcorner = bb(1);
ycorner = bb(2);
xwidth = bb(3);
ywidth = bb(4);
cc = markimg(1).Centroid;
centroidx(i)= cc(1);
centroidy(i)= cc(2);
hold on
rectangle('Position',[xcorner ycorner xwidth ywidth],'EdgeColor','b');
hold on
plot(centroidx(i),centroidy(i), 'bx');
kalmanx = centroidx(i)- xcorner;
kalmany = centroidy(i)- ycorner;
if kfinit == 0
predicted =[centroidx(i),centroidy(i),0,0]' ;
else
predicted = A*actual(i-1,:)';
end
kfinit = 1;
Ppre = A*P*A' + Q;
K = Ppre*H'/(H*Ppre*H'+R);
actual(i,:) = (predicted + K*([centroidx(i),centroidy(i)]' - H*predicted))';
P = (eye(4)-K*H)*Ppre;
hold on
rectangle('Position',[(actual(i,1)-kalmanx)...
(actual(i,2)-kalmany) xwidth ywidth],'EdgeColor','r','LineWidth',1.5);
hold on
plot(actual(i,1),actual(i,2), 'rx','LineWidth',1.5);
drawnow;
end
130 comments: