: In matlab, implement the mean shift clustering algorithm as discussed in lectu

Question

: In matlab, implement the mean shift clustering algorithm as discussed in lectures and papers. Hand-in your code for this question. To do this, use a “flat” kernel function (you will need to specify the value for the radius parameter, ). You can choose to implement as either a “blurring” or “non-blurring” process. To test your algorithm, create some 2-D datasets using matlab’s Gaussian random number generator randn: o a = randn(200,2); o b = a + 4; o c = a; o c(:,1) = 3*c(:,1); o c = c - 4; o d = [a; b]; o e = [a; b; c]; o plot(a(:,1),a(:,2),’+’); o hold on o plot(b(:,1),b(:,2),’o’); o plot(c(:,1),c(:,2),’*’);

Explanation / Answer

Solution :

function [y, MS] = meanShiftPixCluster(x,hs,hr,th,plot_On)
% FUNCTION: meanShiftPixCluster implements classic mean shift pixel
% clustering algorithm introduced in Cmaniciu etal.'s PAMI paper
% "Mean shift: a robust apporach toward feature space analysis", 2002.
% -------------------------------------------------------------------------
% Input:
% x = an input image (either gray or rgb, please expect long time processing if image size is large)
% hs = the bandwidth of spatial kernel (see Eq.(35) in the cited paper)
% hr = the bandwidth of feature kernel (see Eq.(35) in the cited paper)
% th = the threshod of the convergence crIterion (default = .25)
% plot_On = switch on/off the image display of intermediate results (default = 1)
%
% Output:
% y = the output pixel clustered image
% MS = the output of averaged mean shift

% -------------------------------------------------------------------------

%% Checking of arguments
if narg_in<3
error('Type help for the function syntax')
elseif narg_in == 3
th = 1/100; plot_On = 1;
elseif narg_in == 4
if th<0 || th >255
error('threshold range should be in [0,255]')
else
plot_On = 1;
end
elseif narg_in == 5
if sum(ismember(plot_On,[0,1])) == 0
error('plot_On option is to be either 0 or 1')
end
elseif narg_in>5
error('Too many input_arguments')
end
%% initialize
x = double(x);
[height,width,depth] = size(x);
y = x;
Done = 0;
Iter = 0;
if plot_On
figure(randi(1000)+1000);
end
% padding image is to deal with the pixels on borders
xPad = padarray(x,[height,width,0],'symmetric');
% build up the look up table to boost the computation speed
Weightmap = exp( -(0:255^2)/hr^2 );
MS = [];
%% main loop
while ~Done
Weight_accum = 0;
yAccum = 0;
% only 99.75% area of the entire non-zero Gaussian kernel is to be considered
for i = -hs:hs
for j = -hs:hs
if ( i~=0 || j~=0 )
  
spatialKernel = 1;
% uncomment the following line to activate the Gausian kernel
% spatialKernel = exp(-(i^2+j^2)/(hs/3)^2/2);
xThis = xPad(height+i:2*height+i-1, width+j:2*width+j-1, 1:depth);
xDiffSq = (y-xThis).^2;
  
intensityKernel = repmat( prod( reshape( Weightmap( xDiffSq+1 ), height, width, depth) , 3 ), [1,1, depth]);
  
weightThis = spatialKernel.*intensityKernel;
  
Weight_accum = Weight_accum+ weightThis;
  
yAccum = yAccum+xThis.*weightThis;
end
end
end
  
yThis = yAccum./(Weight_accum+eps);

yMS = mean(abs(round(yThis(:))-round(y(:))));
y = round(yThis);
MS(Iter+1) = yMS;
% xPad = padarray(y,[height,width,0],'symmetric');
if plot_On
subplot(121), imshow(uint8(y)),axis image, title(['Iteration times = ' num2str(Iter) '; averaged mean-shift = ' num2str(yMS)]);
subplot(122), plot(0:Iter, MS ), xlabel('Iteration #'), ylabel('averaged mean shift');axis square
drawnow
end
if yMS <= th % exit if converge
Done = 1;
else % otherwise,update the estimated y and repeat the mean shift
Iter = Iter+1;
end
end

// ** Thank You ** //

Related Questions

Navigate

• Browse (All)
• Browse #
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.