Mathematics - Mean & Median Background subtraction

In Mean Background subtraction or mean filter background subtraction the background is the mean of the previous n frames

( Σ i=0 to n-1 (It- i(x,y)) )/ n

In Median background subtraction the background is the median of the previous n frames

median{ I t-i (x,y) }


i ∈ { 0,1,2,3, ... n-1 }


The mathematics that i used in approximate median background subtraction was:

if It (x,y) > B t(x,y) ⇒ B t(x,y) + 1

if It (x,y) < B t(x,y) ⇒ B t(x,y) - 1



This method will eventually converge where half of the pixels are greater than the background and another half less than the background.

Results can be found here:
Result

Mathematics - Basic Background subtraction

| It (x,y) - Bt(x,y) | > T

Basic background subtraction get an image I at time t which is the observed image and subtracts its pixels (x,y) against a background image B at time t. It then compares the absolute value of it against a threshold value which is determined via trial and error.

The basic formula that i used was to subtract the pixels of an image I at time t against an image I at time t+1 and compare it against a threshold.

| It (x,y) - It+1(x,y) | > T

My results can be seen in the following pages:

Result 1

Result 2

Camshift Tracking algorithm

Camshift stands for "Continuously Adaptive Mean Shift."
It has the basic Mean shift algorithm with the difference of a window that changes in size.

Pro: This method is fast and appears on initial testing to be moderately accurate. It may be possible to improve accuracy by using a different color representation.

Con: There are quite a few parameters: the number of histogram bins, the minimum saturation, minimum and maximum intensity, and the width-to-height ratio for regions. There's also a parameter for enlarging the region while doing Mean Shift to increase the chances of finding the maximum for probability density.



using original code from (Isaac Gerg, Adam Ickes, Jamie McCulloch) "http://www.gergltd.com/cse486/project5/" with some minor changes/refactoring done by me.

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;

Blob Detection using Filters

TO DO:

% ALGORITHM
% 1. Generate a Laplacian of Gaussian filter.
% 2. Build a Laplacian scale space, starting with some initial scale and going for n iterations:
% 1. Filter image with scale-normalized Laplacian at current scale.
% 2. Save square of Laplacian response for current level of scale space.
% 3. Increase scale by a factor k.
% 3. Perform nonmaximum suppression in scale space.
% 4. Display resulting circles at their characteristic scales.

video cropping

Velocity

calculate the middle mass in frame 1

wait X seconds

calculate the middle mass in frame 2

speed = (mm_frame_1 - mm_frame_2) * distance / per_pixel

Lens radial distortion can be modelled by the following equations:

x_actual = xd * (1 + distortion_constant * (xd^2 + yd^2))
y_actual = yd * (1 + distortion_constant * (xd^2 + yd^2))

frame subtraction using adaptive thresholding

code:

clear;close all;
im1=imread('stab1.jpg');
im2=imread('stab2.jpg');
bwim1=adaptivethreshold(im1,15,0.08,0);
bwim2=adaptivethreshold(im2,15,0.03,0);
subplot(3,2,1);
imshow(im1);title('First Frame');
subplot(3,2,2);
imshow(bwim1);title('Adaptive Threshold First Frame');
subplot(3,2,3);
imshow(im2);title('Second Frame');
subplot(3,2,4);
imshow(bwim2);title('Adaptive threshold second Frame');

subt = imabsdiff(im1,im2);
subplot(3,2,5);
imshow(subt);title('subtracted bg');

subO = imabsdiff(bwim1,bwim2);
subplot(3,2,6);
imshow(subO);title('subtracted bg Threshold');

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);