How can I count the number of circles with particular radius and center which overlap, without a loop? (I need to speed it up)

The section of code below works to count how many times a particular grid point (with grid points much smaller than the circles themselves) is contained within a circle (i.e. how many circles overlap at each grid point) by making a mask for each individual circle and then adding the total to a matrix "cnts". My problem is that I need to do this for potentially hundreds of thousands of circles. Is there a way to speed up the masking/counting process without using the for loop (which I'm assuming is taking most of the time)? Thanks for any suggestions!

---------------------------------------------------------------------

%% loop over each circle to create a mask and count how many circles hit each grid point

cnts = zeros(size(xx));

for i = 1:length(xc)

xc = x_centers(i); % find the center point of circle #i

yc = y_centers(i);

mask = ((xx-xc).^2 + (yy-yc).^2)<(beam_diameter/2)^2; %create a mask for that circle

% add counts to a grid to determine how many circles hit each grid point

cnts = cnts+mask;

end

Answer 1

[Ganesh Regoti]

From your question, you need an alternative which works as for loop and more optimized. The concept of vectorization works well for your case.

You can refer more about vectorization here

https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

[Christine Z]

Thanks for the reply. We tried another method using circshift but initially ended up with more for loops so it was slower (Method 2). Then instead we tried to vectorize the circle itself and still use circshift but couldn't figure out how to get rid of the loop to add the shift of each individual circle. (see code below 'Method 3: circshift method2') so it wasn't significantly faster either. The code comparing everything we have tried so far is below. Let us know if you have any other suggestions!

close all

clear all

clc

% Method 1: Binary mask circle by circle method

disp('circle by circle mask')

tic

x_centers = linspace(10,50, numshifts);

y_centers = linspace(15,17,numshifts);

beam_diameter = 10;

[xx yy] = meshgrid(0:1:100, 0:1:100); % create a grid

cnts = zeros(size(xx));

for i = 1:numshifts

xc = x_centers(i); % find the center point of each pulse/circle

yc = y_centers(i);

mask = ((xx-xc).^2 + (yy-yc).^2)<(beam_diameter/2)^2; %create a mask for that pulse

% add counts to a grid to determine how many pulses are received for

% each grid point

cnts = cnts+mask;

% disp([num2str(i) ' of ' num2str(numshifts)])

end

figure

imagesc(cnts);

toc

% Method 2: Circshift method

disp('Circshift')

S=strel('disk',5,0);

M=zeros(100,100);

ii=1;

jj=1;

numshifts = 100;

tic

for i=1:numshifts

M=circshift(M,[ii 0]);

for j=1:20

M=circshift(M,[0 jj]);

M(1:11,1:11)=M(1:11,1:11)+S.Neighborhood;

end

% disp([num2str(i) ' of ' num2str(numshifts)])

end

figure

imagesc(M);

toc

% Method 3: Circshift method 2

% close all

% clear all

% clc

%I use this to get the indices in “lindex” (you could use a formula):

% Circshift method

disp('Circshift')

S=strel('disk',5,0);

M=zeros(100,100);

Mtrail=M; Mtrail(1,1)=1;

Mtrack=M;

ii=1;

jj=1;

numshifts = 100;

tic

cnt=0

for i=1:numshifts

M=circshift(M,[ii 0]);

Mtrail=circshift(Mtrail,[ii,0]);

for j=1:20

cnt=cnt+1;

M=circshift(M,[0 jj]);

Mtrail=circshift(Mtrail,[0,jj]);

Mtrack(Mtrail==1)=cnt;

lindex(cnt)=find(Mtrail==1);

M(1:11,1:11)=M(1:11,1:11)+S.Neighborhood;

end

end

figure

imagesc(M);

toc

%And then I use this with “lindex” variable. Output is not perfect, but very close. Unfortunately, no speed improvement (comparing a clean version of the above with the below) :-/

disp('Circshift v2')

M=zeros(100*100,1);

S=strel('disk',5,0);

addNbhd=zeros(100,100);

ii=1;

jj=1;

addNbhd(1:11,1:11)=S.Neighborhood;

addNbhd=addNbhd(:);

lindex2=diff(lindex);

lindex2=[lindex(1),lindex2];

tic

for i=1:20*numshifts

M=circshift(M,lindex2(i))+addNbhd;

end

toc

figure

imagesc(reshape(M,[100,100]));