Mean shift - Toy Example

meanshift.m

clc; close all; clear all;

% Generate data
org = zeros(200,200,3);
org(100:198,100:198,:) = 255;
I = org(:,:,1);
imshow(I)

tic
%           %
% Constants %
%           %
% Threshold of convergence (in pixels for each deg. of freedom)
T = 1;
% Number of pixels to expand search window by.
P = 1;

%            %
% Mean Shift %
%            %

% Initial search window size 
% and location of search window.
% [xmin ymin width height]
sBox = [2, 2, 98, 98]; % Object location from the previous frame/location

% Get the center of the search window
cY = sBox(1) + sBox(3)/2;
cX = sBox(2) + sBox(4)/2;

% To draw the bounding (search) box
shapeInserter = vision.ShapeInserter('BorderColor','Custom','CustomBorderColor',uint8([255 0 0]));

% Initialization
oldBox = sBox(1:2);
MeanConverging = true;

filename = 'meanshift.gif';

% Create search window box on image.
R = step(shapeInserter, org, int32(sBox));
R = insertMarker(R,[cY cX]);

imshow(R); pause(1);

frame = getframe(1);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);
imwrite(imind,cm,filename,'gif', 'Loopcount',inf);

   
MAXITER = 100;
for i = 1:MAXITER

    % Compute centroid of search window
    roiM = [sBox(1)-P, sBox(2)-P, sBox(3)+2*P, sBox(4)+2*P];

    % Moment - sum up the 1's within the roi
    M00 = sum(sum(double(imcrop(I,roiM))));

    % Moment along rows
    % Multiply the pixels by row index then sum up within the roi
    iMultiple = roiM(2):roiM(2)+roiM(4);
    M10 = sum(sum(bsxfun(@times,iMultiple',double(imcrop(I,roiM)))));

    % Moment along cols
    % Multiply the pixels by column index then sum up within the roi
    jMultiple = roiM(1):roiM(1)+roiM(3);
    M01 = sum(sum(bsxfun(@times,jMultiple',double(imcrop(I,roiM))')));
    
    % Get the center of the search window
    cX = round(M10 / M00);
    cY = round(M01 / M00);

    oldBox = sBox(1:2);
    % Update the search box location
    sBox(1:2) = [floor(cY - (sBox(3)/2)), floor(cX - (sBox(4)/2))];

    % Check threshold
    % If the rect location does not change its converged
    if( abs(oldBox(1)-sBox(1)) < T || abs(oldBox(2)-sBox(2)) < T )
        break;
    end
    
    % Debug
    disp(['Iteration : ',num2str(i)])
    disp(['Zeroth moment : ',num2str(M00)]);
    disp(['First moments : ',num2str(M10),',' +num2str(M01)]);
    disp(['Center : ',num2str(cY),',',num2str(cX)])
    disp(['Location : ',num2str(sBox(1)),',',num2str(sBox(2)),',',num2str(sBox(3)),',',num2str(sBox(4))]);
    
    % Create search window box on image.
    R = step(shapeInserter, org, int32([sBox(1) sBox(2) sBox(3)+1 sBox(4)+1]));
    R = insertMarker(R,[cY cX]);
    
    imshow(R); pause(1);
    
    frame = getframe(1);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    imwrite(imind,cm,filename,'gif','WriteMode','append');
   
end