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