Details behind SIFT Feature detection

First octave blur progressively
second octave: resize to half the image and blur progressively

Getting LoG. Subtract image 1 from image 2. Then subtract image 2 from image 3. then subtract image 3 from image 4
Find the maxima/minima of images from one scale higher and lower.

NOW FOR THE DETAILS:(based on my understanding)
images and details from :
Lowe, David G. “Distinctive Image Features from Scale­ Invariant Keypoints”. International Journal of Computer Vision, 60, 2 (2004)

Scale invariant feature transform is used for detection and extracting local features of an image.
the steps are as follows:

1. Generate a Difference of Gaussian(DoG) or a laplacian pyramid
2. Extrema detection from the DoG pyramid which is the local maxima and minima, the point found is an extrema
3. Eliminate low contrast or poorly localized points, what remains are the keypoints
4. Assign an orientation to the points based on the image properties
5. Compute and generate keypoint descriptors

The details:

The main SIFT implementation consists of 4 major stages described by Lowe . The first step is finding a scale space extrema. This is done via image pyramids. The process involves repeatedly smoothing an image by convolving it with a gaussian operator and subsequently sub-sampling in order to achieve higher levels of the pyramid. For this stage, Lowe suggests 4 octaves and 5 levels of blur. An octave corresponds to doubling the value of σ.



After obtaining a full pyramid, the difference of each octave results in an approximate solution to the Laplacian of Gaussian.



The amount of blur per interval and octave is very important in obtaining keypoints (matlab source screenshot)


this is a manual computation i did in excel


The next step is to iterate through each pixel and compare it with its 8 surounding neighbors and 9 neighbors at one scale higher and lower. If the pixel value is higher or lower amongst its neighbors then it is considered as a keypoint.



After determining the keypoint, tylor expansion is used to generate subpixel values from the pixel data.

Doing so, will generate a lot of keypoints. These keypoints are low contrast points or points along an edge. To eliminate these points, tylor expansion is used to get an intensity value for each subpixel. It is then compared against a threshold. If the value is higher or lower than the threshold, the point is accepted or rejected. Following this step a series of stable keypoints will be left that are scale invariant. To make the keypoints rotation invariant, a weighted histogram of local gradient direction and magnitudes around each keypoint is computed at a selected scale and the most prominent orientation in that region is selected and assigned to the keypoint.


The orientation of these keypoints are broken into 36 bins, each representing 10 degrees from the 360 degree. After creating the histogram any peaks above 80 are converted to another keypoint with a different orientation.
To do this, check the highest peek in the histogram and get the left and right neighbors. you will get 3 points, fit a downward parabola and solve the 3 equations to get a new keypoint. do this for all octaves and all iterations.


The next step is to produce unique ID’s for each keypoint that will differentiate them from each other. To accomplish this, a 16x16 window is created between pixels around the keypoint. It is then split into sixteen 4x4 windows and an 8 bin histogram is generated. The histogram and gradient values are interpolated to produce a 128 dimensional feature vector which is invariant to affine illumination. To obtain illumination invariance, the descriptor is normalized by the square root of the sum of squared components

The keypoints are now invariant to affine transformation and partial illumination changes.
Although SIFT has been used in the computer vision community for quite some time there are several studies that try to improve on it.

Disparity Map

In stereo vision we have 2 camera's at a constant distance from each other. we take a picture of a scene. the image that we get from each camera is slightly different from each other. if we place both of them together we will see that they are a little blurred. some of the objects in the scene will look mis-aligned. between the left and right image. if the alignment of an object in the scene is very BIG, then the object is closer to the camera's if the object is less aligned or clear or lined up better, that means its further away from the camera.













Im using the dataset available here:

DATASET HERE >>>
H. Hirschmüller and D. Scharstein. Evaluation of cost functions for stereo matching.
In IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR 2007), Minneapolis, MN, June 2007.

this is because i dont have 2 cameras to take pictures.
they also keep a list of algorithms on who has the best implementation: HERE>>



data is the similarity of the right image from the left, the code uses left to right. i made a mistake in correcting the image dataset. i thought the image left side is from the left camera while its from the other infact.

comparing the wrong images will give unwanted results. the left image and the right image should clearly be defined in order for the disparity map to be correct.
another important factor is the disparity range, i used 0-16 with a window size of 9x9. you can see the results are not so clear.

for this disparity map image i used the 9x9 window with a disparity range of 0-16. the similarity measurement which was used is the sum of square difference.

also the higher the disparity range the longer it will take to process. so im guessing each image should have a different disparity range just like different threshold in background subtraction. you can see from the image below that the disparity range gave a different result.


Its clear that a higher disparity range will give a better result as you can see from the image below i used a disparity range of 0-4 with the same window size of 9x9

CODE TO FOLLOW

Meanshift Tracking algorithm

Mean shift

The main function of this algorithm is histogram estimation. Since moving objects can be identified by their color histogram. Mean-shift tracking algorithm is an iterative scheme based on comparing the histogram of the original object in the current image frame and histogram of candidate regions in the next image frame. The aim is to maximize the correlation between two histograms.

Code:


Mv = aviread('StabilizationResult.avi')
rmin = 0; %min row value for search window
rmax = 0; %max row value for search window
cmin = 0; %min col value for search window
cmax = 0; %max col value for search window
numofframes = 0;
frameNo = 10; %starting frame
vidSize=[128 160]; %video size
centerold = [0 0];
centernew = [0 0];

M =Mv;

% get number of frames
numberofframes = length(M);

Frame1 = M(frameNo);
Image1 = imresize(Frame1.cdata,vidSize,'bilinear');
disp('Click and drag mouse for an initial box window');
% get search window for first frame
[ cmin, cmax, rmin, rmax ] = select( Image1 );
cmin = round(cmin);
cmax = round(cmax);
rmin = round(rmin);
rmax = round(rmax);
wsize(1) = abs(rmax - rmin);
wsize(2) = abs(cmax - cmin);

hue=Image1(:,:,1);
histogram = zeros(256);

for i=rmin:rmax
for j=cmin:cmax
index = uint8(hue(i,j)+1);
%count number of each pixel
histogram(index) = histogram(index) + 1;
end
end

% for each frame
for i = frameNo:numberofframes
Frame = M(i);
I = imresize(Frame.cdata,vidSize,'bilinear');

hue= I(:,:,1);
[rows cols] = size(hue);
probmap = zeros(rows, cols);
for r=1:rows
for c=1:cols
if(hue(r,c) ~= 0)
probmap(r,c)= histogram(hue(r,c));
end
end
end
probmap = probmap/max(max(probmap));
probmap = probmap*255;

count = 0;

rowcenter = 0; % any number just so it runs through at least twice
colcenter = 0;
rowcenterold = 30;
colcenterold = 30;
while (((abs(rowcenter - rowcenterold) > 2) && (abs(colcenter - colcenterold) > 2)) || (count < 15) )
rowcenterold = rowcenter;
colcenterold = colcenter;

[ rowcenter colcenter M00 ] = meanshift(rmin, rmax, cmin,...
cmax, probmap);

rmin = round(rowcenter - wsize(1)/2);
if rmin<1 br=""> rmin=1;
end
rmax = round(rowcenter + wsize(1)/2);
if rmax<1 br=""> rmax=1;
end
cmin = round(colcenter - wsize(2)/2);
if cmin<1 br=""> cmin=1;
end
cmax = round(colcenter + wsize(2)/2);
if cmax<1 br=""> cmax=1;
end
wsize(1) = abs(rmax - rmin);
wsize(2) = abs(cmax - cmin);

count = count + 1;
end

G = I;
trackim=G;

%make box of current search window on saved image
for r= rmin:rmax
trackim(r, cmin:cmin+2) = 0;
trackim(r, cmax-2:cmax) = 0;
end
for c= cmin:cmax
trackim(rmin:rmin+2, c) = 0;
trackim(rmax-2:rmax, c) = 0;
end

windowsize = 100 * (M00/256)^.5;
sidelength = sqrt(windowsize);
rmin = round(rowcenter-sidelength/2);
if rmin<1 br=""> rmin=1;
end
rmax = round(rowcenter+sidelength/2);
if rmax<1 br=""> rmax=1;
end
cmin = round(colcenter-sidelength/2);
if cmin<1 br=""> cmin=1;
end
cmax = round(colcenter+sidelength/2);
if cmax<1 br=""> cmax=1;
end
wsize(1) = abs(rmax - rmin);
wsize(2) = abs(cmax - cmin);


outname = sprintf('./answer/%d.jpg', i);
imwrite(trackim, outname);

figure(1),imshow(trackim); hold on
end
hold off;

Tracking using Optical Flow

Optical flow with Hierarchical Lucas Kanade (using pyramids)

results are very strange!














Horn’s optical flow algorithm. Results are unclear.




Code:

function [u,v,cert] = HierarchicalLK(im1, im2, numLevels, windowSize, iterations, display)
%HIERARCHICALLK Hierarchical Lucas Kanade (using pyramids)
% [u,v]=HierarchicalLK(im1, im2, numLevels, windowSize, iterations, display)

if (size(im1,1)~=size(im2,1)) | (size(im1,2)~=size(im2,2))
error('images are not same size');
end;

if (size(im1,3) ~= 1) | (size(im2, 3) ~= 1)
error('input should be gray level images');
end;


% check image sizes and crop if not divisible
if (rem(size(im1,1), 2^(numLevels - 1)) ~= 0)
warning('image will be cropped in height, size of output will be smaller than input!');
im1 = im1(1:(size(im1,1) - rem(size(im1,1), 2^(numLevels - 1))), :);
im2 = im2(1:(size(im1,1) - rem(size(im1,1), 2^(numLevels - 1))), :);
end;

if (rem(size(im1,2), 2^(numLevels - 1)) ~= 0)
warning('image will be cropped in width, size of output will be smaller than input!');
im1 = im1(:, 1:(size(im1,2) - rem(size(im1,2), 2^(numLevels - 1))));
im2 = im2(:, 1:(size(im1,2) - rem(size(im1,2), 2^(numLevels - 1))));
end;

%Build Pyramids
pyramid1 = im1;
pyramid2 = im2;

for i=2:numLevels
im1 = reduce(im1);
im2 = reduce(im2);
pyramid1(1:size(im1,1), 1:size(im1,2), i) = im1;
pyramid2(1:size(im2,1), 1:size(im2,2), i) = im2;
end;

% base level computation
disp('Computing Level 1');
baseIm1 = pyramid1(1:(size(pyramid1,1)/(2^(numLevels-1))), 1:(size(pyramid1,2)/(2^(numLevels-1))), numLevels);
baseIm2 = pyramid2(1:(size(pyramid2,1)/(2^(numLevels-1))), 1:(size(pyramid2,2)/(2^(numLevels-1))), numLevels);
[u,v] = LucasKanade(baseIm1, baseIm2, windowSize);

for r = 1:iterations
[u, v] = LucasKanadeRefined(u, v, baseIm1, baseIm2);
end

%propagating flow 2 higher levels
for i = 2:numLevels
disp(['Computing Level ', num2str(i)]);
uEx = 2 * imresize(u,size(u)*2); % use appropriate expand function (gaussian, bilinear, cubic, etc).
vEx = 2 * imresize(v,size(v)*2);

curIm1 = pyramid1(1:(size(pyramid1,1)/(2^(numLevels - i))), 1:(size(pyramid1,2)/(2^(numLevels - i))), (numLevels - i)+1);
curIm2 = pyramid2(1:(size(pyramid2,1)/(2^(numLevels - i))), 1:(size(pyramid2,2)/(2^(numLevels - i))), (numLevels - i)+1);

[u, v] = LucasKanadeRefined(uEx, vEx, curIm1, curIm2);

for r = 1:iterations
[u, v, cert] = LucasKanadeRefined(u, v, curIm1, curIm2);
end
end;

if (display==1)
figure, quiver(reduce((reduce(medfilt2(flipud(u),[5 5])))), -reduce((reduce(medfilt2(flipud(v),[5 5])))), 0), axis equal
end


Code:

%Implementation based on Horn’s optical flow algorithm.
function [u, v] = LucasKanade(im1, im2, windowSize);
%LucasKanade lucas kanade algorithm, without pyramids (only 1 level);

%REVISION: NaN vals are replaced by zeros

[fx, fy, ft] = ComputeDerivatives(im1, im2);

u = zeros(size(im1));
v = zeros(size(im2));

halfWindow = floor(windowSize/2);
for i = halfWindow+1:size(fx,1)-halfWindow
for j = halfWindow+1:size(fx,2)-halfWindow
curFx = fx(i-halfWindow:i+halfWindow, j-halfWindow:j+halfWindow);
curFy = fy(i-halfWindow:i+halfWindow, j-halfWindow:j+halfWindow);
curFt = ft(i-halfWindow:i+halfWindow, j-halfWindow:j+halfWindow);

curFx = curFx';
curFy = curFy';
curFt = curFt';


curFx = curFx(:);
curFy = curFy(:);
curFt = -curFt(:);

A = [curFx curFy];

U = pinv(A'*A)*A'*curFt;

u(i,j)=U(1);
v(i,j)=U(2);
end;
end;

u(isnan(u))=0;
v(isnan(v))=0;

%u=u(2:size(u,1), 2:size(u,2));
%v=v(2:size(v,1), 2:size(v,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fx, fy, ft] = ComputeDerivatives(im1, im2);
%ComputeDerivatives Compute horizontal, vertical and time derivative
% between two gray-level images.

if (size(im1,1) ~= size(im2,1)) | (size(im1,2) ~= size(im2,2))
error('input images are not the same size');
end;

if (size(im1,3)~=1) | (size(im2,3)~=1)
error('method only works for gray-level images');
end;


fx = conv2(double(im1),double(0.25* [-1 1; -1 1])) + conv2(double(im2), double(0.25*[-1 1; -1 1]));
fy = conv2(double(im1), double(0.25*[-1 -1; 1 1])) + conv2(double(im2), double(0.25*[-1 -1; 1 1]));
ft = conv2(double(im1), double(0.25*ones(2))) + conv2(double(im2),double( -0.25*ones(2)));

% make same size as input
fx=fx(1:size(fx,1)-1, 1:size(fx,2)-1);
fy=fy(1:size(fy,1)-1, 1:size(fy,2)-1);
ft=ft(1:size(ft,1)-1, 1:size(ft,2)-1);



TO DO: More experiments




5 Frame difference optical flow using hierarchical lucas kanade on stabilized image sequence.
[u,v,w] = HierarchicalLK(ai1,ai2,3,1,1,1);
results show noise.




2 consecutive frames, stablized vs. Original
stabilized has more inconsistent movement! (Results are inconclusive)

Adaptive Thesholding

An adaptive thresholding algorithm that seperates the foreground from the background with nonuniform illumination.

Thresholding is used to segment an image by setting all pixels whose intensity values are above a threshold to a foreground value and all the remaining pixels to a background value.

Whereas the conventional operator uses a global threshold for all pixels, adaptive thresholding changes the threshold dynamically over the image. This more sophisticated version of thresholding can accommodate changing lighting conditions in the image, e.g. those occurring as a result of a strong illumination gradient or shadows.

Code:
for i=1:1:100 %(number of image)

file_name=strcat('',int2str(i),'.jpg');
im1=imread(file_name);

image=adaptivethreshold(im1,5,0.4,1); %additional adaptive threshold


SE = strel('square', 15);
biner_dilate=imdilate(image,SE); %pixel dilatation

biner_filter= bwareaopen(biner_dilate,150);

SE = strel('square', 20);
biner_dilate_2=imdilate(biner_filter,SE); %pixel dilatation


region=regionprops(biner_dilate_2); %regionprops

koordinat_x=region(1,1).BoundingBox(1,1);
koordinat_y=region(1,1).BoundingBox(1,2);
width_x=region(1,1).BoundingBox(1,3);
width_y=region(1,1).BoundingBox(1,4);

center_x=region(1,1).Centroid(1,1);
center_y=region(1,1).Centroid(1,2);


figure(1)

frame_count=strcat('frame ke :',int2str(i));
subplot(2,2,1); imshow(im1); title(frame_count);

subplot(2,2,2); imshow(image); title('adaptive threshold');

subplot(2,2,3); imshow(biner_dilate_2); title('test dilate 2');

subplot(2,2,4); imshow(biner_filter); title('filter');

subplot(2,2,1),rectangle('Position',[koordinat_x,koordinat_y,width_x,width_y], 'Edge', 'g', 'LineWidth',2);
hold on
subplot(2,2,1),plot(center_x,center_y,'g+');
hold off

end

function bw=adaptivethreshold(IM,ws,C,tm)
%ADAPTIVETHRESHOLD An adaptive thresholding algorithm that seperates the
%foreground from the background with nonuniform illumination.
% bw=adaptivethreshold(IM,ws,C) outputs a binary image bw with the local
% threshold mean-C or median-C to the image IM.
% ws is the local window size.
% tm is 0 or 1, a switch between mean and median. tm=0 mean(default); tm=1 median.
%
% Contributed by Guanglei Xiong (xgl99@mails.tsinghua.edu.cn)
% at Tsinghua University, Beijing, China.
%
% For more information, please see
% http://homepages.inf.ed.ac.uk/rbf/HIPR2/adpthrsh.htm

if (nargin<3)
error('You must provide the image IM, the window size ws, and C.');
elseif (nargin==3)
tm=0;
elseif (tm~=0 && tm~=1)
error('tm must be 0 or 1.');
end

IM=mat2gray(IM);

if tm==0
mIM=imfilter(IM,fspecial('average',ws),'replicate');
else
mIM=medfilt2(IM,[ws ws]);
end
sIM=mIM-IM-C;
bw=im2bw(sIM,0);
bw=imcomplement(bw);

Object Tracking in a LIVE VIDEO STREAM

% %run first in command
% ONCE PER MATLAB SESSION
% vid = videoinput('winvideo', '1', 'YUY2_160x120');
% set(vid,'ReturnedColorSpace','rgb');
% set(vid,'TriggerRepeat',Inf);
% vid.FrameGrabInterval = 5;
% start(vid);
%------

figure;
while(vid.FramesAcquired<=1000) % Stop after 1000 frames data = getdata(vid,2); diff_im = imabsdiff(data(:,:,:,1),data(:,:,:,2)); %background subtraction
diff = rgb2gray(diff_im);
diff_bw = im2bw(diff,0.2);
bw2 = imfill(diff_bw,'holes');
s = regionprops(diff_im, 'centroid');
cd = s.Centroid;
centroids = cat(1, s.Centroid);
imshow(data(:,:,:,2));
hold(imgca,'on');
plot(imgca,centroids(:,1),centroids(:,2),'g*');

hold on;
rectangle('Position',[cd(:,1) cd(:,2) 20 20],'LineWidth',2,'EdgeColor','b');
hold(imgca,'off');

end

stop(vid)

Matlab morphological operation

IM2 = imdilate(IM,SE) dilates the grayscale, binary, or packed binary image IM, returning the dilated image, IM2.
IM2 = imerode(IM,SE) erodes the grayscale, binary, or packed binary image IM, returning the eroded image IM2.


The argument SE is a structuring element object or array of structuring element objects returned by the strel function.


The morphological close operation is a dilation followed by an erosion, using the same structuring element for both operations. The morphological open operation is an erosion followed by a dilation, using the same structuring element for both operations.

these operations are performed on a stabilized video sequence done by SIFT. and basic background subtraction.

se1 = strel('square',11) % 11-by-11 square
se2 = strel('line',10,45) % line, length 10, angle 45 degrees
se3 = strel('disk',15) % disk, radius 15
se4 = strel('ball',15,5) % ball, radius 15, height 5


sedisk = strel('square',10);
fg = imclose(fg, sedisk);


sedisk = strel('line',10,90);
fg = imclose(fg, sedisk);

sedisk = strel('line',15,5);
fg = imclose(fg, sedisk);


sedisk = strel('disk',10);
fg = imclose(fg, sedisk);


Deciding on the morphological structure is much harder than deciding on the threshold of background subtraction

High complexity background subtraction using mixture of guassian

The mean u of each Gaussian function, can be thought of as an educated guess of the pixel value in the next frame—we assume here that pixels are usually background. The weight and standard deviations of each component are measures of our confidence in that guess (higher weight & lower σ = higher confidence). There are typically 3-5 Gaussian components per pixel—the number typically depending on memory limitations.

In the image, the left images show the startup where all pixels are assumed to be background.through out the course of the video the background will be constructed. in contrast to median filter where it estimates and clears up the background.

Code:

% ----------------------- frame size variables -----------------------

fr = source(1).cdata; % read in 1st frame as background frame
fr_bw = rgb2gray(fr); % convert background to greyscale
fr_size = size(fr);
width = fr_size(2);
height = fr_size(1);
fg = zeros(height, width);
bg_bw = zeros(height, width);

% --------------------- mog variables -----------------------------------

C = 3; % number of gaussian components (typically 3-5)
M = 3; % number of background components
D = 2.5; % positive deviation threshold
alpha = 0.01; % learning rate (between 0 and 1) (from paper 0.01)
thresh = 0.25; % foreground threshold (0.25 or 0.75 in paper)
sd_init = 6; % initial standard deviation (for new components) var = 36 in paper
w = zeros(height,width,C); % initialize weights array
mean = zeros(height,width,C); % pixel means
sd = zeros(height,width,C); % pixel standard deviations
u_diff = zeros(height,width,C); % difference of each pixel from mean
p = alpha/(1/C); % initial p variable (used to update mean and sd)
rank = zeros(1,C); % rank of components (w/sd)


% --------------------- initialize component means and weights -----------

pixel_depth = 8; % 8-bit resolution
pixel_range = 2^pixel_depth -1; % pixel range (# of possible values)

for i=1:height
for j=1:width
for k=1:C

mean(i,j,k) = rand*pixel_range; % means random (0-255)
w(i,j,k) = 1/C; % weights uniformly dist
sd(i,j,k) = sd_init; % initialize to sd_init

end
end
end

%--------------------- process frames -----------------------------------

for n = 1:length(source)

fr = source(n).cdata; % read in frame
fr_bw = rgb2gray(fr); % convert frame to grayscale

% calculate difference of pixel values from mean
for m=1:C
u_diff(:,:,m) = abs(double(fr_bw) - double(mean(:,:,m)));
end

% update gaussian components for each pixel
for i=1:height
for j=1:width

match = 0;
for k=1:C
if (abs(u_diff(i,j,k)) <= D*sd(i,j,k)) % pixel matches component

match = 1; % variable to signal component match

% update weights, mean, sd, p
w(i,j,k) = (1-alpha)*w(i,j,k) + alpha;
p = alpha/w(i,j,k);
mean(i,j,k) = (1-p)*mean(i,j,k) + p*double(fr_bw(i,j));
sd(i,j,k) = sqrt((1-p)*(sd(i,j,k)^2) + p*((double(fr_bw(i,j)) - mean(i,j,k)))^2);
else % pixel doesn't match component
w(i,j,k) = (1-alpha)*w(i,j,k); % weight slighly decreases

end
end

w(i,j,:) = w(i,j,:)./sum(w(i,j,:));

bg_bw(i,j)=0;
for k=1:C
bg_bw(i,j) = bg_bw(i,j)+ mean(i,j,k)*w(i,j,k);
end

% if no components match, create new component
if (match == 0)
[min_w, min_w_index] = min(w(i,j,:));
mean(i,j,min_w_index) = double(fr_bw(i,j));
sd(i,j,min_w_index) = sd_init;
end

rank = w(i,j,:)./sd(i,j,:); % calculate component rank
rank_ind = [1:1:C];

% sort rank values
for k=2:C
for m=1:(k-1)

if (rank(:,:,k) > rank(:,:,m))
% swap max values
rank_temp = rank(:,:,m);
rank(:,:,m) = rank(:,:,k);
rank(:,:,k) = rank_temp;

% swap max index values
rank_ind_temp = rank_ind(m);
rank_ind(m) = rank_ind(k);
rank_ind(k) = rank_ind_temp;

end
end
end

% calculate foreground
match = 0;
k=1;

fg(i,j) = 0;
while ((match == 0)&&(k<=M))

if (w(i,j,rank_ind(k)) >= thresh)
if (abs(u_diff(i,j,rank_ind(k))) <= D*sd(i,j,rank_ind(k)))
fg(i,j) = 0;
match = 1;
else
fg(i,j) = fr_bw(i,j);
end
end
k = k+1;
end
end
end

figure(1),subplot(3,1,1),imshow(fr)
subplot(3,1,2),imshow(uint8(bg_bw))
subplot(3,1,3),imshow(uint8(fg))
end

medium complexity background subtraction using approximate median

In this method the previous N frames of video are buffered, and the background is calculated as the median of buffered frames. The problem is the background is cleared of all objects after few frames to show a cleared background.

Median filtering has been shown to be very robust and to have performance comparable to higher complexity methods. However, storing and processing many frames of video (as is often required to track slower moving objects) requires an often prohibitively large amount of memory. This can be alleviated somewhat by storing and processing frames at a rate lower than the frame rate— thereby lowering storage and computation requirements at the expense of a slower adapting background.

The approximate median method works as such: if a pixel in the current frame has a value larger than the corresponding background pixel, the background pixel is incremented by 1. Likewise, if the current pixel is less than the background pixel, the background is decremented by one. In this way, the background eventually converges to an estimate where half the input pixels are greater than the background, and half are less than the background—approximately the median (convergence time will vary based on frame rate and amount movement in the scene.)




for i = 2:length(source)

fr = source(i).cdata;
fr_bw = rgb2gray(fr); % convert frame to grayscale

fr_diff = abs(double(fr_bw) - double(bg_bw)); % avoid negative overflow

for j=1:width
for k=1:height

if ((fr_diff(k,j) > thresh))
fg(k,j) = fr_bw(k,j);
else
fg(k,j) = 0;
end

if (fr_bw(k,j) > bg_bw(k,j))
bg_bw(k,j) = bg_bw(k,j) + 1;
elseif (fr_bw(k,j) <>
bg_bw(k,j) = bg_bw(k,j) - 1;
end

end
end

figure(1),subplot(3,1,1),imshow(fr)
subplot(3,1,2),imshow(uint8(bg_bw))
subplot(3,1,3),imshow(uint8(fg))
end