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.

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 w/ blob detection, morphological operation (Togeather)

frames = {avi.cdata}; %uses the cdata from the video file

fg = extractForeground(frames); % do foreground extraction
cmap = colormap(gray);

for i = 1:length(fg)
temp0{i} = edge(fg{i}, 'canny', 0.99) + fg{i};
temp2 = temp0{i};
temp2 = cat(3,temp2,temp2,temp2);

fgs = rgb2gray(temp2);
sedisk = strel('square',10);
fgs = imclose(fgs, sedisk);
fgs = imfill(fgs,'holes');
RLL = bwlabel(fgs);

stats = regionprops(RLL,'basic','Centroid');
fig = figure(1),imshow(RLL)
hold on

for n = 1:length(stats)
if(stats(n).Area > 100)
plot(stats(n).Centroid(1), stats(n).Centroid(2),'r*')
end
end
hold off


end;

clear all;

Identify and track the center of the logical object

after morphological operation.
check if the area of blob is greater than a threshold

sedisk = strel('square',15);
fg = imclose(fg, sedisk);
fg = imfill(fg,'holes');
RLL = bwlabel(fg);

stats = regionprops(RLL,'basic','Centroid');
figure(1),imshow(fr_bw)
hold on

for n = 1:length(stats)
if(stats(n).Area > 100)
plot(stats(n).Centroid(1), stats(n).Centroid(2),'r*')
end
end

Problem: can track less 50% of the cars, however there are too many outiners because of the transformation done during stabilization. also hard to determine the area, i did it using trial and error to get the best fit

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)

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

Low complexity background subtraction using frame difference method

Frame differencing, also known as temporal difference, uses the video frame at time t-1 as the background model for the frame at time t. This technique is sensitive to noise and
variations in illumination, and does not consider local consistency
properties of the change mask.
This method also fails to segment the non-background objects if they stop moving. Since it uses only a single previous frame, frame differencing may not be able to identify the interior
pixels of a large, uniformly-colored moving object. This is commonly known as the aperture problem.


a major flaw of this method is that for objects with uniformly distributed intensity values, the pixels are interpreted as part of the background. Another problem is that objects must be continuously moving. If an object stays still for more than a frame period (1/fps), it becomes part of the background.
This method does have two major advantages. One obvious advantage is the modest computational load. Another is that the background model is highly adaptive. Since the background is based solely on the previous frame, it can adapt to changes in the background faster than any other method (at 1/fps to be precise). As we'll see later on, the frame difference method subtracts out extraneous background noise (such as waving trees), much better than the more complex approximate median and mixture of Gaussians methods.
A challenge with this method is determining the threshold value.

The video is a result of stabilization using SIFT features.

Code:

clear all;close all;clc;
source = aviread('stabilized');
thresh = 40;
bg = source(1).cdata; % read in 1st frame as background frame
bg_bw = rgb2gray(bg); % convert background to greyscale
% ----------------------- set frame size variables -----------------------
fr_size = size(bg);
width = fr_size(2);
height = fr_size(1);
fg = zeros(height, width);
% --------------------- process frames -----------------------------------
for i = 2:length(source)
fr = source(i).cdata; % read in frame
fr_bw = rgb2gray(fr); % convert frame to grayscale
fr_diff = abs(double(fr_bw) - double(bg_bw));
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
end
end
bg_bw = fr_bw;
figure(1),subplot(3,1,1),imshow(fr)
subplot(3,1,2),imshow(fr_bw)
subplot(3,1,3),imshow(uint8(fg))
end

Video stabilization - using sift

first video shows Original video with Unwanted camera shake
second video shows result of video without camera shake




Buit in matlab function:
CP2TFORM Infer spatial transformation from control point pairs.
CP2TFORM takes pairs of control points and uses them to infer a
spatial transformation.


CP2TFORM requires a minimum number of control point pairs to infer a
% TFORM structure of each TRANSFORMTYPE:
%
% TRANSFORMTYPE MINIMUM NUMBER OF PAIRS
% ------------- -----------------------
% 'nonreflective similarity' 2
% 'similarity' 3
% 'affine' 3
% 'projective' 4


http://www.mathworks.com/access/helpdesk/help/toolbox/images/cp2tform.html

given a set of points "inp" in the first image
another set of points "outp" in another image
recompute the change in the projected points

Code:

% compute the similarity transformation
tt = cp2tform(inp, outp, 'linear conformal');
L = tt.tdata.T;
scale = sqrt(L(1,1)^2 + L(1,2)^2);
L(:,1:2) = L(:,1:2)/scale;
% sum the transforms up
T = T + L;
end

% get the average transformation
T = T/count;

% warping by backprojection.
% inv(T) is the transform from target back to the original image
T = inv(T);

width = size(im, 2);
height = size(im, 1);
% x and y are coordinates on the warped image
[x,y] = meshgrid(1:width, 1:height);
nxy = [x(:), y(:), ones(length(x(:)), 1)] * T;
% newx and newx are corresponding point coordinates in the original image
newx = reshape(nxy(:,1), height, width);
newy = reshape(nxy(:,2), height, width);

Math:

For an ‘affine’ transformation, the parametric motion can be described by the following formulas:

u(x,y) = a1x + a2y + a3
v(x,y) = a4x + a5y + a6

In this case there are six unknowns namely a1 –> a6 . In order to solve for these unknowns we must expand matrices to this form:







Since u1 –> u3 and v1 –> v3 can be easily calculated using the 3 pairs of points from each image and (x1 –> x3 , y1 –> y3) are already known, a1 –> a6 can be solved by using the formula v = (AT*A)-1 * AT b where (AT*A)-1 is defined as the pseudo inverse of A. With these new values imtransform() can modify the first image to correspond with the second image. xdata and ydata are variables describing the offset between the static image and transformed image, while trans is the scaled and sheared transformed image.


For a projective transformation: [up vp wp] = [x y w] T, where

u = up / wp
v = vp / wp.


T is a 3-by-3 matrix, where all nine elements can be different.

T = [ A D G
B E H
C F I]


The above matrix equation is equivalent to these two expressions:

u = (Ax + By + C) / (Gx + Hy + 1)
v = (Dx + Ey + F) / (Gx + Hy + 1)

Summary:

For a projective transformation:

u = (Ax + By + C)/(Gx + Hy + I)
v = (Dx + Ey + F)/(Gx + Hy + I)

Assume I = 1, multiply both equations, by denominator:

u = [x y 1 0 0 0 -ux -uy] * [A B C D E F G H]'
v = [0 0 0 x y 1 -vx -vy] * [A B C D E F G H]'

With 4 or more correspondence points we can combine the u equations and
the v equations for one linear system to solve for [A B C D E F G H]:

[ u1 ] = [ x1 y1 1 0 0 0 -u1*x1 -u1*y1 ] * [A]
[ u2 ] = [ x2 y2 1 0 0 0 -u2*x2 -u2*y2 ] [B]
[ u3 ] = [ x3 y3 1 0 0 0 -u3*x3 -u3*y3 ] [C]
[ u1 ] = [ x4 y4 1 0 0 0 -u4*x4 -u4*y4 ] [D]
[ ... ] [ ... ] [E]
[ un ] = [ xn yn 1 0 0 0 -un*xn -un*yn ] [F]
[ v1 ] = [ 0 0 0 x1 y1 1 -v1*x1 -v1*y1 ] [G]
[ v2 ] = [ 0 0 0 x2 y2 1 -v2*x2 -v2*y2 ] [H]
[ v3 ] = [ 0 0 0 x3 y3 1 -v3*x3 -v3*y3 ]
[ v4 ] = [ 0 0 0 x4 y4 1 -v4*x4 -v4*y4 ]
[ ... ] [ ... ]
[ vn ] = [ 0 0 0 xn yn 1 -vn*xn -vn*yn ]

Or rewriting the above matrix equation:
U = X * Tvec, where Tvec = [A B C D E F G H]'
so Tvec = X\U.

Corner Detection using Harris

Harris and Stephens improved upon Moravec's corner detector by considering the differential of the corner score with respect to direction directly, instead of using shifted patches. (This corner score is often referred to as autocorrelation, since the term is used in the paper in which this detector is described. However, the mathematics in the paper clearly indicate that the sum of squared differences is used.)

C. Harris and M. Stephens (1988). "A combined corner and edge detector" (PDF). Proceedings of the 4th Alvey Vision Conference. pp. 147--151. http://www.csse.uwa.edu.au/~pk/research/matlabfns/Spatial/Docs/Harris/A_Combined_Corner_and_Edge_Detector.pdf.

How does it work?
The first step in Harris is to compute a corner response function. Harris uses a series of filters that are considered the most important steps to find a corner. These kernels are:

• A presmoothing filter pfilt = {0.223755f,0.552490f,0.223755f};
• a gradient filter gfilt = {0.453014f,0.0f,-0.453014f};
• a Blur Filter bfilt = {0.01563f,0.09375f,0.234375f,0.3125f, 4 0.234375f,0.09375f,0.01563f};

The combination of these filters and few parameters are what make Harris successful in finding corners. The parameters are as followed::

• Steering parameter 0.04 - 0.06 (0.05 default)
• Response threshold 10K - 1M (default 25,000)
• neighborhood radius 10 pixels

clear all
close all
clc;

disp('Harris Corner detection Purposes Only....');

fin = 'smallVersion2.avi'; %input File name to be changed
fout = 'movie2.avi'; %output File name to be changed

%get the File info
fileinfo = aviinfo(fin);
%get the number of Frames
nframes = fileinfo.NumFrames;

aviobj = avifile(fout, 'compression', 'none', 'fps',fileinfo.FramesPerSecond);

for i = 1:nframes %You may need to limit the nFrames, may result in a large file size
%Read frames from input video
mov_in = aviread(fin,i);
im_in = frame2im(mov_in);

%Do processing on each frame of the video

im = rgb2gray(im_in);
% Find Harris corners in image1 and image2
if ~isa(im,'double')
im = double(im);
end

subpixel = nargout == 5;

dx = [-1 0 1; -1 0 1; -1 0 1]; % Derivative masks
dy = dx';

Ix = conv2(im, dx, 'same'); % Image derivatives
Iy = conv2(im, dy, 'same');

% Generate Gaussian filter of size 6*sigma (+/- 3sigma) and of
% minimum size 1x1.
sigma =1;
g = fspecial('gaussian',max(1,fix(6*sigma)), sigma);

Ix2 = conv2(Ix.^2, g, 'same'); % Smoothed squared image derivatives
Iy2 = conv2(Iy.^2, g, 'same');
Ixy = conv2(Ix.*Iy, g, 'same');
k = 0.04;
cim = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2; % Original Harris measure.


if nargin > 2 % We should perform nonmaximal suppression and threshold

if subpixel
[r,c,rsubp,csubp] = nonmaxsuppts(cim, 1, 10000, im);
else
[r,c] = nonmaxsuppts(cim, 1, 10000, im);
end
end


show(im,1), hold on, plot(c,r,'r+');
drawnow


%Write frames to output video
%frm = im2frame(im_out);
%aviobj = addframe(aviobj,frm);

%i %Just to display frame number you can ommit
end;

%Don't forget to close output file
aviobj = close(aviobj);
return;




function [r,c, rsubp, csubp] = nonmaxsuppts(cim, radius, thresh, im)

subPixel = nargout == 4; % We want sub-pixel locations
[rows,cols] = size(cim);

% Extract local maxima by performing a grey scale morphological
% dilation and then finding points in the corner strength image that
% match the dilated image and are also greater than the threshold.

sze = 2*radius+1; % Size of dilation mask.
mx = ordfilt2(cim,sze^2,ones(sze)); % Grey-scale dilate.

% Make mask to exclude points within radius of the image boundary.
bordermask = zeros(size(cim));
bordermask(radius+1:end-radius, radius+1:end-radius) = 1;

% Find maxima, threshold, and apply bordermask
cimmx = (cim==mx) & (cim>thresh) & bordermask;

[r,c] = find(cimmx); % Find row,col coords.


if subPixel % Compute local maxima to sub pixel accuracy
if ~isempty(r) % ...if we have some ponts to work with

ind = sub2ind(size(cim),r,c); % 1D indices of feature points
w = 1; % Width that we look out on each side of the feature
% point to fit a local parabola

% Indices of points above, below, left and right of feature point
indrminus1 = max(ind-w,1);
indrplus1 = min(ind+w,rows*cols);
indcminus1 = max(ind-w*rows,1);
indcplus1 = min(ind+w*rows,rows*cols);

% Solve for quadratic down rows
cy = cim(ind);
ay = (cim(indrminus1) + cim(indrplus1))/2 - cy;
by = ay + cy - cim(indrminus1);
rowshift = -w*by./(2*ay); % Maxima of quadradic

% Solve for quadratic across columns
cx = cim(ind);
ax = (cim(indcminus1) + cim(indcplus1))/2 - cx;
bx = ax + cx - cim(indcminus1);
colshift = -w*bx./(2*ax); % Maxima of quadradic

rsubp = r+rowshift; % Add subpixel corrections to original row
csubp = c+colshift; % and column coords.
else
rsubp = []; csubp = [];
end
end

if nargin==4 & ~isempty(r) % Overlay corners on supplied image.
figure(1), imshow(im,[]), hold on
if subPixel
plot(csubp,rsubp,'r+'), title('corners detected');
else
plot(c,r,'r+'), title('corners detected');
end
end

Optical Flow using Lucas Kanade

Optical flow is the pattern of apparent motion of objects, surfaces, and edges in a visual scene caused by the relative motion between an observer (an eye or a camera) and the scene.

The optical flow methods try to calculate the motion between two image frames which are taken at times t and t + 1 at every voxel position.

One method for determining optical flow is the Differential methods of estimation, based on partial derivatives of the image signal and/or the sought flow field and higher-order partial derivatives. the method ill be using is the Lucas–Kanade Optical Flow Method –

Lucas Kanade with Pyramids Algorithm:

Problem to solve: Video stabilization

TO DO: MATLAB CODE

Video - Hough Transform

The hough function implements the Standard Hough Transform (SHT). The Hough transform is designed to detect lines, using the parametric representation of a line:

rho = x*cos(theta) + y*sin(theta)

The variable rho is the distance from the origin to the line along a vector perpendicular to the line. theta is the angle between the x-axis and this vector.
The hough function generates a parameter space matrix whose rows and columns correspond to these rho and theta values, respectively. The houghpeaks function finds peak values in this space, which represent potential lines in the input image.
The houghlines function finds the endpoints of the line segments corresponding to peaks in the Hough transform and it automatically fills in small gaps.

The code:

clear all
close all
clc;
disp('Testing Purposes Only....');
%fin = 'sampleVideo.avi';
fin = 'smallVersion.avi';
fout = 'test2.avi';
avi = aviread(fin);

% Convert to RGB to GRAY SCALE image.
avi = aviread(fin);
pixels = double(cat(4,avi(1:2:end).cdata))/255; %get all pixels (normalize)

nFrames = size(pixels,4); %get number of frames
for f = 1:nFrames
pixel(:,:,f) = (rgb2gray(pixels(:,:,:,f))); %convert images to gray scale
end
rows=128;
cols=160;
nrames=f;
for l = 2:nrames

%edge detection
edgeD(:,:,l) = edge(pixel(:,:,l),'canny');
g(:,:,l) = double(edgeD(:,:,l));

%subtract background
d(:,:,l)=(abs(pixel(:,:,l)-pixel(:,:,l-1))); %subtract current pixel from previous

%convert to binary image
k=d(:,:,l);
bw(:,:,l) = im2bw(k, .2);
bw1=bwlabel(bw(:,:,l));

[H,theta,rho] = hough(edgeD(:,:,l));

P = houghpeaks(H,5,'threshold',ceil(0.1*max(H(:))));
x = theta(P(:,2));
y = rho(P(:,1));
lines = houghlines(edgeD(:,:,l),theta,rho,P,'FillGap',5,'MinLength',7);
imshow(pixel(:,:,l)); hold on;
max_len = 0;
for k = 1:length(lines)

xy = [lines(k).point1; lines(k).point2];
plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');
% Plot beginnings and ends of lines
plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');
len = norm(lines(k).point1 - lines(k).point2);

if ( len > max_len)
max_len = len;
xy_long = xy;
end
end

% highlight the longest line
plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color','cyan');
plot(xy_long(:,1),xy_long(:,2),'x','LineWidth',2,'Color','yellow');
plot(xy_long(:,1),xy_long(:,2),'x','LineWidth',2,'Color','red');
drawnow;
hold off

end

results:

The same process on an Image:

%# load image, process it, find edges

I = rgb2gray( imread('pillsetc.png') );
I = imcrop(I, [30 30 450 350]);
J = imfilter(I, fspecial('gaussian', [17 17], 5), 'symmetric');
BW = edge(J, 'canny');
%# Hough Transform and show matrix
[H T R] = hough(BW);
imshow(imadjust(mat2gray(H)), [], 'XData',T, 'YData',R, 'InitialMagnification','fit')
xlabel('\theta (degrees)'), ylabel('\rho')
axis on, axis normal, hold on
colormap(hot), colorbar
%# detect peaks
P = houghpeaks(H, 4);
plot(T(P(:,2)), R(P(:,1)), 'gs', 'LineWidth',2);
%# detect lines and overlay on top of image
lines = houghlines(BW, T, R, P);
figure, imshow(I), hold on

for k = 1:length(lines)
xy = [lines(k).point1; lines(k).point2];
plot(xy(:,1), xy(:,2), 'g.-', 'LineWidth',2);
end
hold off

Video - Boundary Tracing

The MATLAB toolbox includes two functions you can use to find the boundaries of objects in a binary image:

• bwtraceboundary
• bwboundaries

The bwtraceboundary function returns the row and column coordinates of all the pixels on the border of an object in an image. You must specify the location of a border pixel on the object as the starting point for the trace.
The bwboundaries function returns the row and column coordinates of border pixels of all the objects in an image. For both functions, the nonzero pixels in the binary image belong to an object and pixels with the value 0 (zero) constitute the background.

The code:

clear all
close all
clc;
disp('Testing Purposes Only....');
%fin = 'sampleVideo.avi';
fin = 'smallVersion2.avi';
fout = 'test2.avi';
avi = aviread(fin);


% Convert to RGB to GRAY SCALE image.
avi = aviread(fin);
pixels = double(cat(4,avi(1:2:end).cdata))/255; %get all pixels (normalize)

nFrames = size(pixels,4); %get number of frames
for f = 1:nFrames
pixel(:,:,f) = (rgb2gray(pixels(:,:,:,f))); %convert images to gray scale
end
rows=128;
cols=160;
nrames=f;
for l = 2:nrames

%edge detection
edgeD(:,:,l) = edge(pixel(:,:,l),'canny');
g(:,:,l) = double(edgeD(:,:,l));

%subtract background
d(:,:,l)=(abs(pixel(:,:,l)-pixel(:,:,l-1))); %subtract current pixel from previous

%convert to binary image
k=d(:,:,l);
bw(:,:,l) = im2bw(k, .2);
bw1=bwlabel(bw(:,:,l));

%By default, bwboundaries finds the boundaries of all objects in an image
BW_filled = imfill(bw(:,:,l),'holes');
boundaries = bwboundaries(BW_filled);

imshow(pixel(:,:,l));

hold on
for k=1:10
b = boundaries{k};
plot(b(:,2),b(:,1),'g','LineWidth',3);
end
drawnow;
hold off

end

Results:

Merge a sequence of video images to another (red Channel)

Combining 2 different sequence of images into one via a red channel for distinction.

The results:


The Code:


clear all
close all
clc;
disp('Testing Purposes Only....');
fin = 'smallVersion.avi';
fout = 'bgSubtract1.avi';
avi = aviread(fin);

% Convert to RGB to GRAY SCALE image.
avi = aviread(fin);
pixels = double(cat(4,avi(1:2:end).cdata))/255; %get all pixels (normalize)
nFrames = size(pixels,4); %get number of frames

for f = 1:nFrames
pixel(:,:,f) = (rgb2gray(pixels(:,:,:,f))); %convert images to gray scale
end

nrames=f;
for l = 2:nrames
d(:,:,l)=(abs(pixel(:,:,l)-pixel(:,:,1))); %subtract current pixel from background
z(:,:,l)=(abs(pixel(:,:,l)-pixel(:,:,l-1))); %subtract current pixel from previous


%-------seperate channel
I =z(:,:,l);
rz = cat(3,I,I,I);
for i=1:128 %video size
for j=1:160 %video size
if rz(i,j,1) > 0.2 %theshold
rz(i,j,1) = 1; %red channel
else
rz(i,j,1) = 0;
end
rz(i,j,2) = 0; %green channel
rz(i,j,3) = 0; %Blue Channel
end
end

%------merge
fg= rz;
bg= cat(3,pixels(:,:,l),pixels(:,:,l),pixels(:,:,l)); %array of 3
coef = 0.6;
dif = fg-bg;

out = bg + coef.*dif;

imshow(out);
%-----
drawnow;
%hold off

end

Merge a sequence of video images to another

In motion detection you would want to combine your video to the original video. to do this, you combine your original video with the processed video.

the code:

clear all
close all
clc;
disp('Testing Purposes Only....');
fin = 'sampleVideo.avi';
avi = aviread(fin);
% Convert to RGB to GRAY SCALE image.
avi = aviread(fin);
pixels = double(cat(4,avi(1:2:end).cdata))/255; %get all pixels (normalize)
nFrames = size(pixels,4);
%get number of frames

for f = 1:nFrames
pixel(:,:,f) = (rgb2gray(pixels(:,:,:,f))); %convert images to gray scale
end

nrames=f;
for l = 2:nrames
d(:,:,l)=(abs(pixel(:,:,l)-pixel(:,:,1))); %subtract current pixel from background
z(:,:,l)=(abs(pixel(:,:,l)-pixel(:,:,l-1))); %subtract current pixel from previous

%------merge
fg= pixel(:,:,l); %foreground
bg= z(:,:,l); %background
alpha= 0.1; %alpha
dif = fg-bg;
out = bg + alpha.*dif;
imshow(out);
%--------
drawnow;
end

Results:

Motion Detection - background image subtraction

One of the most common approaches is to compare the current frame with the previous one. It's useful in video compression when you need to estimate changes and to write only the changes, not the whole frame. But it is not the best one for motion detection applications.

here the process i used:

1. get each frames pixel and normalize to 255
2. Convert it to gray scale
3. subtract pixels from previous pixels (previous post subtracted images)
4. show each image in a sequence

heres the code:

clear all
close all
clc;
disp('Testing Purposes Only....');
fin = 'smallVersion.avi';
avi = aviread(fin);

% Convert to RGB to GRAY SCALE image.
avi = aviread(fin);
pixels = double(cat(4,avi(1:2:end).cdata))/255; %get all pixels (normalize)
nFrames = size(pixels,4); %get number of frames

for f = 1:nFrames
pixel(:,:,f) = (rgb2gray(pixels(:,:,:,f))); %convert images to gray scale
end

nrames=f;
for l = 2:nrames
%d(:,:,l)=(abs(pixel(:,:,l)-pixel(:,:,1))); %subtract current pixel from background
z(:,:,l)=(abs(pixel(:,:,l)-pixel(:,:,l-1))); %subtract current pixel from previous pixel

%imshow(d(:,:,l));
imshow(z(:,:,l));

hold on
drawnow;
hold off

end

Results:

Image 1: Frame subtracted from Background or first frame (overlaps)

Image 2: Frame subtracted from previous Frame

Video - Motion Detection

There are many approaches for motion detection in a continuous video stream. All of them are based on comparing of the current video frame with one from the previous frames or with something that we'll call background.

as you can see from my implementation, the results are not very promising, because of the colors in the video stream. the next approach is using image and converting to gray scale.

here the code for background subtraction using video frames:

clear all
close all
clc;
disp('Testing Purposes Only....');
fin = 'sampleVideo.avi';
fout = 'test2.avi';
fileinfo = aviinfo(fin);
nframes = fileinfo.NumFrames;
aviobj = avifile(fout, 'compression', 'none', 'fps',fileinfo.FramesPerSecond);
for i = 2:nframes %1
%Read frames from input video
mov_in = aviread(fin,i);
im_in = frame2im(mov_in);
%Do processing on each frame of the video
%----------------------------------------------------------------------

mov_in2 = aviread(fin,i-1); %or 1 for the first frame
im_in2 = frame2im(mov_in2);
background = imopen(im_in2,strel('disk',15)); %get the background image
I2= imsubtract(im_in,background);
im_out =I2;
%----------------------------------------------------------------------
%Write frames to output video
frm = im2frame(im_out);
aviobj = addframe(aviobj,frm);
%i %Just to display frame number
end;
%Don't forget to close output file
aviobj = close(aviobj);
msgbox('DONE');
return;

results:

Video Editting using Image Frames

one easy way for image processing and adding extra images is:

1. read each frame of the video
2. perform image analysis
3. show the image and hold
4. add additional and stop hold
5. continue drawing the image

this way youll see the images drawn on screen in a sequence that represents a video. the next stop is exporting it back to video. here is the structure im using:

clear all
close all
clc;
disp('Testing Purposes Only....');
fin = 'sampleVideo.avi';
fout = 'test2.avi';
avi = aviread(fin);

%view output frame by frame in an image
video = {avi.cdata};
for a = 1:length(video)

%---------------------------------
%Do image processing

newImage = rgb2gray(video{a});
%---------------------------------

imshow(newImage); %or use imagesc
axis image off
hold on

%---------------------------------
%add tracing to image
rectangle('Position',[50 50 70 100],'EdgeColor','r');
%---------------------------------

drawnow;
hold off
end;

Video Spatial Transformation - Rotation

To rotate an image in Matlab, use the imrotate function. imrotate accepts two primary arguments:

• The image to be rotated
• The rotation angle

You specify the rotation angle in degrees. If you specify a positive value, imrotate rotates the image counterclockwise; if you specify a negative value, imrotate rotates the image clockwise.

This example rotates the image I 35 degrees in the counterclockwise direction.

J = imrotate(I,35);

As optional arguments to imrotate, you can also specify

• The interpolation method
• The size of the output image

In Video:

fin = 'rawVideo.avi';
fout = 'test.avi';
fileinfo = aviinfo(fin);
nframes = fileinfo.NumFrames;
aviobj = avifile(fout, 'compression', 'none', 'fps',fileinfo.FramesPerSecond);
for i = 1:20
%Read frames from input video
mov_in = aviread(fin,i);
im_in = frame2im(mov_in);
%Do processing on each frame of the video
%----------------------------------------------------------------------
%In this example - Image Rotation
im_out = imrotate(im_in,35);
%----------------------------------------------------------------------
%Write frames to output video
frm = im2frame(im_out);
aviobj = addframe(aviobj,frm);
%i %Just to display frame number
end;
%Don't forget to close output file
aviobj = close(aviobj);

msgbox('DONE');
return;

The Result:

Video Spatial Transformation - Interpolation

Problems arise with the resizing of video. By default, imresize uses nearest-neighbor interpolation to determine the values of pixels in the output image, but you can specify other interpolation methods.

Interpolation is the process used to estimate an image value at a location in between image pixels.

For example, if you resize an image so it contains more pixels than it did originally, the toolbox uses interpolation to determine the values for the additional pixels. The imresize and imrotate geometric functions use two-dimensional interpolation as part of the operations they perform. The improfile image analysis function also uses interpolation.

The interpolation methods all work in a fundamentally similar way. In each case, to determine the value for an interpolated pixel, they find the point in the input image that the output pixel corresponds to. They then assign a value to the output pixel by computing a weighted average of some set of pixels in the vicinity of the point. The weightings are based on the distance each pixel is from the point.

The methods differ in the set of pixels that are considered:

• For nearest-neighbor interpolation, the output pixel is assigned the value of the pixel that the point falls within. No other pixels are considered.
• For bilinear interpolation, the output pixel value is a weighted average of pixels in the nearest 2-by-2 neighborhood.
• For bicubic interpolation, the output pixel value is a weighted average of pixels in the nearest 4-by-4 neighborhood.

The number of pixels considered affects the complexity of the computation. Therefore the bilinear method takes longer than nearest-neighbor interpolation, and the bicubic method takes longer than bilinear. However, the greater the number of pixels considered, the more accurate the effect is, so there is a tradeoff between processing time and quality.

Heres the list of argument values you can use:

1. 'nearest' Nearest-neighbor interpolation (the default)
2. 'bilinear' Bilinear interpolation
3. 'bicubic' Bicubic interpolation

Example for Images:

I = imread('circuit.tif');
J = imresize(I ,[100 150],'bilinear');
imshow(I)
figure, imshow(J)

In Videos:

fin = 'rawVideo.avi'; fout = 'test.avi'; fileinfo = aviinfo(fin); nframes = fileinfo.NumFrames; aviobj = avifile(fout, 'compression', 'none', 'fps',fileinfo.FramesPerSecond); for i = 1:20 %Read frames from input video mov_in = aviread(fin,i); im_in = frame2im(mov_in); %Do processing on each frame of the video
%----------------------------------------------------------------------
%In this example - Image Resize using interpolation
im_out = imresize(im_in,1.5,'nearest'); %----------------------------------------------------------------------
%Write frames to output video frm = im2frame(im_out); aviobj = addframe(aviobj,frm); %i %Just to display frame number end; %Don't forget to close output file aviobj = close(aviobj); msgbox('DONE'); return;

the Results:

Image 1: Nearest


Image 2: Bilinear


Image 3: Bicubic