Trying to speed up a circle machine

Hello All,

I have an algorithm that will draw a circle for each x and y pair, with a given a radius, within an arbitrary image dimension. I am trying to find overlapping pairs with the given radius.

The algorithm completion time increases with the number of particles. I am trying to find a way to decrease the completion time while keeping the number of coordinates pairs high. It would be great if I could have it run ~2000 pairs for ~10 seconds. Currently, it seems to be ~2000 pairs with about ~20 seconds. The bottle neck is in the first for loop as shown below.

The code is below. I am just stuck on how to revamp the code to make it run faster and I don't know of any resources which I can use to help. I am running with these parameters shown below and with the attached particle list.

Matlab '9.2.0.556344 (R2017a)' Intel® Core™ i5-4200M CPU @ 2.50 GHz RAM 8 GB

>> [C, keepers, overlapXY] = circleMachine_help(1024,6,pkList);
Elapsed time is 22.674214 seconds.
Elapsed time is 0.014354 seconds.

function [C, keepers, overlapXY] = circleMachine_help(imageSize,rejectRadius,pkList)

% Create the grid
[xx, yy] = meshgrid(1:imageSize);

% C = false(imageSize,imageSize);
C = zeros(imageSize,imageSize);
tic
% Create the circles on the grid
for n=1:size(pkList,1)
    C = C + ...
        (sqrt((xx-pkList(n,1)).^2+(yy-pkList(n,2)).^2)<=rejectRadius);
end
toc
% list of pks that overlap
overlapXY = zeros(size(pkList,1),2);

keepers = zeros(size(overlapXY,1),2);
tic
% Check to find particles which overlap if not put them on a keepers list
for n=1:size(pkList,1)
    OverLapImage = C(...
        pkList(n,2)-rejectRadius:pkList(n,2)+rejectRadius,...
        pkList(n,1)-rejectRadius:pkList(n,1)+rejectRadius);

      if isempty(OverLapImage(OverLapImage > 1)) == 0

          overlapXY(n,:) = pkList(n,:);

      else

          keepers(n,:) = pkList(n,:);
      end
  end
  toc
  % clean up overlapXY and keeprs
  overlapXY(overlapXY(:,1) == 0 & overlapXY(:,2) == 0,:) = [];
  keepers(keepers(:,1) == 0 & keepers(:,2) == 0,:) = [];

% display image
figure;
imagesc(C);
hold on;
scatter(overlapXY(:,1),overlapXY(:,2),20,'r');
scatter(keepers(:,1),keepers(:,2),8,'g');
title(['Circles of original pkList.' ...
    'Reds are overlaps and green are keepers.']);
hold off;

Answer 1

[Rik]

You can replace the first loop with this strategy: generate a logical the correct image size, set the center coordinates of each circle to true. Now you can create a structuring element and apply a dilate operation. I don't now if it is indeed faster, but I think it is.

A down-side to this is that if your centers are not integer values, this will return a different matrix. You could solve this by first running the first part for integer positions, and then doing the rest with your original loop.

If you encounter problems writing the code, just let me know.

Answer 2

[Guillaume]

I don't have matlab on this machine to test but replacing your circle creation loop with just

pkList2 = permute(pkList, [3 2 1]);
C = sum(hypot(xx - pkList(1, 1, :), yy - pkList(1, 2, :) <= rejectRadius, 3);

would probably be faster, at the expense of memory usage. The above requires a temporary matrix of size imageSize x ImageSize x size(pkList, 1).

Note that I'm using hypot instead of hand calculating the radial distance. It probably does not make a difference to speed (but is much easier on the eyes). For a probably infinitesimal speed boost, you could remove the sqrt calculation and instead compare to the square of rejectRadius.

Another option for massively boosting your loop would be to just do the comparison over

xx(max(pkList(n, 2)-rejectRadius, 1) : min(pkList(n, 2)+rejectRadius, imageSize),
   max(pkList(n, 1)-rejectRadius, 1) : min(pkList(n, 1)-rejectRadius, imageSize))
%same for yy

Instead of testing imageSize x imageSize points at each step of the loop, only test rejectRadius^2 x rejectRadius^2 points around your centre. The others are guaranteed to be outside the circle.

Answer 3

[bdatko]

Hello Everyone,

Thank you for the help. I have been working through both @Guillaume and @"Rik Wisselink" suggestions. I have also been talking with my peers and with all the help I have reworked the entire algorithm with huge speed up. I think the first version of the code still has merits since visually it makes sense at the cost of speed. Instead of drawing the image I am looping through the list of particles and calculating the distance of every particle from every other particle. The algorithm will find overlapping particles within a chosen radius and inspect which particle to keep by their brightness.

I did notice that the code takes several iterations to converge to the answer. The first iteration is shown below. The first iteration fails to find many overlapping particles. I am not sure why this is the case but maybe more experienced eyes can recognize the source of the issue.

Anyway, below is the code and attached is the data set. You will notice that the code still draws the picture at the end. This is only to verify that search is working. In between the tic toc is the algorithm. Again, the run below is only the first iteration.

LB = CircleMachine_help_v2(6,pkListMaxInt,1024);
Elapsed time is 0.238588 seconds.

function keepers = ...
    CircleMachine_help_v2(rejectRadius,pkListMaxInt,imageSize)

tic

% Calculate the distance of every particle to each other particle,
% including itself
distancePkList = zeros(size(pkListMaxInt,1),size(pkListMaxInt,1));
for jj=1:size(pkListMaxInt,1)
    tmp = sqrt((pkListMaxInt(jj,1)-pkListMaxInt(:,1)).^2+...
        (pkListMaxInt(jj,2)-pkListMaxInt(:,2)).^2);
    distancePkList(jj,:) = tmp;
end
clear tmp;

% find particles that are not overlapping
keepers = zeros(size(pkListMaxInt,1),3);
keepersIndex = 1;
for n=1:size(distancePkList,1)

      if isempty(find(distancePkList(n,:) < (2*rejectRadius + 1) &...
              distancePkList(n,:) ~= 0,1)) == 1
          keepers(keepersIndex,:) = pkListMaxInt(n,:);
          keepersIndex = keepersIndex + 1;
      end

end

% clean up zeros from the storage
keepers(keepers(:,1) == 0 & keepers(:,2) == 0,:) = [];

% remove the lower half of the distancePkList
% This will ensure you don't inspect every particle pair twice
distancePkList=triu(distancePkList);
[rr, cc] = find(distancePkList < (2*rejectRadius + 1) & ...
    distancePkList ~= 0);

RowCol = [rr, cc];
RowCol = sort(RowCol);

% create the Leader Board
LB = zeros(size(pkListMaxInt,1),3);
LBindex = 1;
for n=1:size(RowCol,1)
    players = RowCol(n,:);
    % compete
    if pkListMaxInt(players(1),3) > pkListMaxInt(players(2),3)

          [LB, LBindex ] = CheckList(...
              LB,pkListMaxInt(players(1),:),pkListMaxInt(players(2),:),LBindex);
      else
          [LB, LBindex ] = CheckList(...
              LB,pkListMaxInt(players(2),:),pkListMaxInt(players(1),:),LBindex);

      end
  end

      function [LeaderBoard, index ] = CheckList(...
              LeaderBoard,winner,loser,index)

          if isempty(find(LeaderBoard(:,1) == winner(1) & ...
                  LeaderBoard(:,2) == winner(2),1)) == 0
              return
          end

          if isempty(find(LeaderBoard(:,1) == loser(1) & ...
                  LeaderBoard(:,2) == loser(2),1))
              % Loser wasn't on the list put winner on the list
              LeaderBoard(index,:) = winner;
              index = index + 1;
          else
              % Loser was on list replace with winner
              rows = find(LeaderBoard(:,1) == loser(1) & ...
                  LeaderBoard(:,2) == loser(2));
              LeaderBoard(rows,:) = winner;
          end
      end

LB(LB(:,1) == 0 & LB(:,2) == 0,:) = [];

keepers = [keepers; LB];
toc
% plot
% Create the grid
[xx, yy] = meshgrid(1:imageSize);

% C = false(imageSize,imageSize);
C = zeros(imageSize,imageSize);

% Create the circles on the grid
for n=1:size(pkListMaxInt,1)

      C = C + ...
          (sqrt((xx-pkListMaxInt(n,1)).^2+(yy-pkListMaxInt(n,2)).^2)<=rejectRadius);
  end

% display image
figure;
imshow(C);
hold on;
scatter(keepers(:,1),keepers(:,2),20,'g');
title('White Circles are every particle. Green are the keepers centeres.')
hold off;

end

[Image Analyst]

How do we call that code? You didn't give us input values for rejectRadius,pkListMaxInt,imageSize so we don't know what to pass to the function. The mat file has this

pkListMaxInt: [1945×3 single]

but what about the other inputs???

[bdatko]

I apologize. I call the function shown below. Where 6 is the radius of the particle in pixels and 1024 is the image size in pixels (assuming a square image).

LB = CircleMachine_help_v2(6,pkListMaxInt,1024);