ivote3.m 4.79 KB
``````function [px, py] = ivote3(image, rmax, phi, iter, t0, sigma, t1)

%	This function performs iterative voting on an image in order to identify the centers of blobs

%	image 	=	a matrix of values representing a 2D image
%	rmax	=	estimated radius of a blob
%	phi		=	angular range of the voting area (in radians)
%	iter	=	number of iterations for voting
%	sigma	=	standard deviation of the weighting filter used for voting
%	t0		=	initial gradient magnitude threshold (limits the # of voting pixels)
%	t1		=	final threshold for valid points

dphi = phi / (iter);

%compute the gradient and gradient magnitude
[Igrad_x, Igrad_y] = gradient(image);
Igrad_mag = sqrt(Igrad_x.^2 + Igrad_y.^2);

%represent the gradient direction in polar coordinates
Igrad_theta = atan2(Igrad_x,Igrad_y) + pi;

%compute pixels that will be processed

% calculate a mask representing the voting pixels
S = Igrad_mag > t0;

Sx = size(S, 1);
Sy = size(S, 2);
S(1:rmax, :) = 0;
S(Sx - rmax:Sx, :) = 0;
S(:, 1:rmax) = 0;
S(:, Sy - rmax:Sy) = 0;

%calculate the coordinates of all of the voting pixels

[sx, sy] = find(S);     %this provides a coordinate for each voter

%eliminate pixels close to the image boundary

nV = nnz(S);            %calculate the number of valid pixels (voters)

% create a meshgrid describing coordinates relative to the voter
% position
range = -rmax:rmax;
[Vx, Vy] = meshgrid(range, range);
%create a mask for the voting area
M_dist =  Vx.^2 + Vy.^2 < rmax^2;
M_theta = atan2(Vx, Vy);

pxPerRow = size(image,1);

g_v_prime = zeros(nV, ceil(rmax^2*phi/3));
validPoints = zeros(nV,1);

for i = 1: iter

t_iter = tic;

%set the vote image to zeros
Ivote = zeros(size(image));

for v = 1:nV

%determine the position of the current voter
vx = sx(v);
vy = sy(v);

%determine the orientation (theta) of the current voter
vtheta = Igrad_theta(vx, vy);

% find the angular distance between M_theta and V_theta at each
ang_diff = abs(M_theta - vtheta);
M_diff = min(2*pi - ang_diff,ang_diff)<phi;

%compute the final mask
M = (M_dist .* M_diff);

%get the coordinates of each pixel in the final mask
[vx_prime, vy_prime] = find(M);

%transform the local coordinates of the pixels in the voting region
%to the global coordinates of the image

npts =  numel(vx_prime);
validPoints(v) = npts;

g_v_prime(v,1:npts) = vx + (vx_prime - (rmax + 1)) + (vy + (vy_prime - (rmax + 1))-1).*pxPerRow;

% compute the Gaussian kernel for the current voter
%c_x = (max(vx_prime) + min(vx_prime))/2
%c_y = (max(vy_prime) + min(vy_prime))/2
c_x = rmax/2 * cos(vtheta)+rmax;
c_y = rmax/2 * sin(vtheta)+rmax;
Gauss_size_x = vx_prime -c_x;
Gauss_size_y = vy_prime -c_y;
gauss_p1 = 1/(2*pi*sigma^2);
gauss_p2 = -1/(2*sigma^2);
Gauss_p3 = (Gauss_size_x.^2 + Gauss_size_y.^2);
gaussian = gauss_p1*exp(gauss_p2*Gauss_p3);
% *******************************************
%determine the gradient magnitude at the current voter location
vmag = Igrad_mag(vx, vy);

%add the voter's current magnitude weighted by gaussian function to the vote image
for n=1: npts
%             Ivote( g_v_prime(v,n)) = Ivote( g_v_prime(v,n)) + gaussian(n)*vmag;
Ivote( g_v_prime(v,n)) = Ivote( g_v_prime(v,n)) + vmag;
end

end

%update the voting direction for each voter
for v = 1:nV

%determine the position of the current voter
vx = sx(v);
vy = sy(v);

%get the local value of the voting image
local_Ivote = Ivote(g_v_prime(v,1:validPoints(v)));

%find the index of the maximum value
[~, local_max_idx] = max(local_Ivote);

%convert this into a global subscript
gx = rem(g_v_prime(v,local_max_idx)-1, pxPerRow) + 1;
gy = (g_v_prime(v,local_max_idx) - gx)/pxPerRow + 1;

%compute the vector from the voter position to this position
new_vx = gx - vx;
new_vy = gy - vy;

%convert this to a new theta value
new_vtheta = atan2(new_vy, new_vx);

%update the gradient image with the new voter direction
Igrad_theta(vx, vy) = new_vtheta;

end

%reduce phi for the next iteration
phi = phi - dphi;
toc(t_iter)

end

% add threshold and find the local maxima
cell_center = Ivote;
cell_center (cell_center <t1) = 0;
cell_center = imregionalmax(cell_center);
imshow(cell_center);
%  Create an Overlay
output = mat2gray(image);
output(:, :, 2) = mat2gray(cell_center);
output(:, :, 3) = mat2gray(image);
figure; imagesc(output);
imwrite(output, 'output image.bmp');
[px,py] = find(cell_center);``````