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

3 次查看（过去 30 天）
Christine Z
Christine Z 2019-7-8
评论： Christine Z 2019-7-17
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

请先登录，再回答此问题。

回答（1 个）

Ganesh Regoti
Ganesh Regoti 2019-7-17
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
  1 个评论
Christine Z
Christine Z 2019-7-17
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]));

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Contour Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by