Setting the largest circle in a free area with a distributed particle image

Question

I work on images to detect and set the maximum possible circle in any free area of ​​the image containing distributed particles:

(capable of detecting the location of a particle).

One direction is to identify the circle that touches any three-point combination, checking to see if the circle is free, and then I find the largest circle among all the empty circles. However, this leads to a huge number of combinations, i.e. C(n,3) , where n is the total number of particles in the image.

I would appreciate if anyone could provide me with any hint or alternative method that I can explore.

Answer 1

Let's make some maths my friend, since math will always be through!

Wikipedia:

In mathematics, the Voronoi diagram is a partition of a plane into regions based on the distance to points in a specific subset of the plane.

For example:

 rng(1) x=rand(1,100)*5; y=rand(1,100)*5; voronoi(x,y); 

The best part about this diagram is that if you notice, all the edges / vertices of these blue areas will be equal to the distance to the points around them. Thus, if we know the location of the vertices and calculate the distances to the nearest points, then we can select the vertex with the greatest distance as our center of the circle.

It is interesting to note that the edges of Voronoi regions are also defined as the circles of triangles generated by the Delaunay triangulation.

Therefore, if we calculate the triangulation of the Delaney area and their environment

 dt=delaunayTriangulation([x;y].'); cc=circumcenter(dt); %voronoi edges 

And calculate the distances between the centers and any of the points that define each triangle:

 for ii=1:size(cc,1) if cc(ii,1)>0 && cc(ii,1)<5 && cc(ii,2)>0 && cc(ii,2)<5 point=dt.Points(dt.ConnectivityList(ii,1),:); %the first one, or any other (they are the same distance) distance(ii)=sqrt((cc(ii,1)-point(1)).^2+(cc(ii,2)-point(2)).^2); end end 

Then we have the center ( cc ) and radius ( distance ) of all possible circles that do not have a point inside them. We just need the biggest!

 [r,ind]=max(distance); %Tada! 

Now let's schedule

 hold on ang=0:0.01:2*pi; xp=r*cos(ang); yp=r*sin(ang); point=cc(ind,:); voronoi(x,y) triplot(dt,'color','r','linestyle',':') plot(point(1)+xp,point(2)+yp,'k'); plot(point(1),point(2),'g.','markersize',20); 

Notice how the center of the circle is at one vertex of the Voronoi diagram.

---

NOTE : this will detect the center inside [0-5], [0-5]. You can easily change it to change this restriction. You can also try to find a circle that is fully consistent with the area of ​​interest (as opposed to the center). This will require a small addition at the end, where the maximum will be obtained.

Answer 2

I would like to propose another solution based on refinement grid search. It is not as advanced as Ander, or as shorter as rahnema1, but it should be very easy to follow and understand. In addition, it works pretty fast.

The algorithm contains several stages:

• We generate a uniformly distributed grid.
• We find the minimum distances of points in the grid to all provided points.
• We discard all points whose distances are below a certain percentile (for example, the 95th).
• We select the area that contains the largest distance (this should contain the correct center if my initial mesh is thin enough).
• We create a new meshgrid around the selected area and again find the distances (this part is clearly not optimal, since the distances are calculated for all points, including far-reaching and irrelevant ones).
• We repeat the refinement within the region, while monitoring the change in 5% values ​​→ if it falls below a certain threshold that we break.

A few notes:

• I made the assumption that circles cannot go beyond scattered points (that is, the bounding square of the scattering acts like an "invisible wall").
• The corresponding percentile depends on how good the initial grid is. This will also affect the number of while iterations and the optimal starting value for cnt .

 function [xBest,yBest,R] = q42806059 rng(1) x=rand(1,100)*5; y=rand(1,100)*5; %% Find the approximate region(s) where there exists a point farthest from all the rest: xExtent = linspace(min(x),max(x),numel(x)); yExtent = linspace(min(y),max(y),numel(y)).'; % Create a grid: [XX,YY] = meshgrid(xExtent,yExtent); % Compute pairwise distance from grid points to free points: D = reshape(min(pdist2([XX(:),YY(:)],[x(:),y(:)]),[],2),size(XX)); % Intermediate plot: % figure(); plot(x,y,'.k'); hold on; contour(XX,YY,D); axis square; grid on; % Remove irrelevant candidates: D(D<prctile(D(:),95)) = NaN; D(D > xExtent | D > yExtent | D > yExtent(end)-yExtent | D > xExtent(end)-xExtent) = NaN; %% Keep only the region with the largest distance L = bwlabel(~isnan(D)); [~,I] = max(table2array(regionprops('table',L,D,'MaxIntensity'))); D(L~=I) = NaN; % surf(XX,YY,D,'EdgeColor','interp','FaceColor','interp'); %% Iterate until sufficient precision: xExtent = xExtent(~isnan(min(D,[],1,'omitnan'))); yExtent = yExtent(~isnan(min(D,[],2,'omitnan'))); cnt = 1; % increase or decrease according to the nature of the problem while true % Same ideas as above, so no explanations: xExtent = linspace(xExtent(1),xExtent(end),20); yExtent = linspace(yExtent(1),yExtent(end),20).'; [XX,YY] = meshgrid(xExtent,yExtent); D = reshape(min(pdist2([XX(:),YY(:)],[x(:),y(:)]),[],2),size(XX)); D(D<prctile(D(:),95)) = NaN; I = find(D == max(D(:))); xBest = XX(I); yBest = YY(I); if nanvar(D(:)) < 1E-10 || cnt == 10 R = D(I); break end xExtent = (1+[-1 +1]*10^-cnt)*xBest; yExtent = (1+[-1 +1]*10^-cnt)*yBest; cnt = cnt+1; end % Finally: % rectangle('Position',[xBest-R,yBest-R,2*R,2*R],'Curvature',[1 1],'EdgeColor','r'); 

The result that I get for the data from the Ander example is [x,y,r] = [0.7832, 2.0694, 0.7815] (which is the same thing). Runtime is about half the Ander solution.

Here are the intermediate charts:

The contour of the largest (transparent) distance from a point to the set of all provided points:

After considering the distance from the border, keeping only the top 5% of the deleted points and considering only the area that contains the largest distance (a piece of surface represents the stored values):

And finally:

Answer 3

You can use bwdist from the Image Processing Toolbox to calculate the conversion of the image distance. This can be seen as a method for creating a raven diagram, well explained in @AnderBiguri's answer.

 img = imread('AbmxL.jpg'); %convert the image to a binary image points = img(:,:,3)<200; %compute the distance transform of the binary image dist = bwdist(points); %find the circle that has maximum radius radius = max(dist(:)); %find position of the circle [xy] = find(dist == radius); imshow(dist,[]); hold on plot(y,x,'ro'); 

Answer 4

The fact that this problem can be solved using "direct search" (as can be seen in another answer ) means that we can consider this as a global optimization . There are various ways to solve such problems, each of which is suitable for certain scenarios. Out of my personal curiosity, I decided to solve this problem using a genetic algorithm .

Generally speaking, such an algorithm requires that we consider a solution as a set of “genes” prone to “evolution” with a certain “fitness function”. As it happens, it's pretty easy to identify the genes and fitness function in this problem:

• Genes: x , y , r .
• Fitness function: technically the maximum lap area, but this is equivalent to the maximum r (or minimum -r , since the algorithm requires minimizing the function).
• A particular limitation is that if r greater than the Euclidean distance to the nearest of the provided points (that is, the circle contains a point), the body "dies."

The basic implementation of such an algorithm is given below (“basic” because it is completely non-optimized, and there are many possibilities for ^{no pun} optimization in this task).

 function [x,y,r] = q42806059b(cloudOfPoints) % Problem setup if nargin == 0 rng(1) cloudOfPoints = rand(100,2)*5; % equivalent to Ander initialization. end %{ figure(); plot(cloudOfPoints(:,1),cloudOfPoints(:,2),'.w'); hold on; axis square; set(gca,'Color','k'); plot(0.7832,2.0694,'ro'); plot(0.7832,2.0694,'r*'); %} nVariables = 3; options = optimoptions(@ga,'UseVectorized',true,'CreationFcn',@gacreationuniform,... 'PopulationSize',1000); S = max(cloudOfPoints,[],1); L = min(cloudOfPoints,[],1); % Find geometric bounds: % In R2017a: use [S,L] = bounds(cloudOfPoints,1); % Here we also define distance-from-boundary constraints. g = ga(@(g)vectorized_fitness(g,cloudOfPoints,[L;S]), nVariables,... [],[], [],[], [L 0],[S min(SL)], [], options); x = g(1); y = g(2); r = g(3); %{ plot(x,y,'ro'); plot(x,y,'r*'); rectangle('Position',[xr,yr,2*r,2*r],'Curvature',[1 1],'EdgeColor','r'); %} function f = vectorized_fitness(genes,pts,extent) % genes = [x,y,r] % extent = [Xmin Ymin; Xmax Ymax] % f, the fitness, is the largest radius. f = min(pdist2(genes(:,1:2), pts, 'euclidean'), [], 2); % Instant death if circle contains a point: f( f < genes(:,3) ) = Inf; % Instant death if circle is too close to boundary: f( any( genes(:,3) > genes(:,1:2) - extent(1,:) | ... genes(:,3) > extent(2,:) - genes(:,1:2), 2) ) = Inf; % Note: this condition may possibly be specified using the A,b inputs of ga(). f(isfinite(f)) = -genes(isfinite(f),3); %DEBUG: %{ scatter(genes(:,1),genes(:,2),10 ,[0, .447, .741] ,'o'); % All z = ~isfinite(f); scatter(genes(z,1),genes(z,2),30,'r','x'); % Killed z = isfinite(f); scatter(genes(z,1),genes(z,2),30,'g','h'); % Surviving [~,I] = sort(f); scatter(genes(I(1:5),1),genes(I(1:5),2),30,'y','p'); % Elite %} 

And here is a graph of the “temporary passage” of 47 generations of a typical launch:

(Where the blue dots are the current generation, the red crosses are “killed by insta” organisms, the green hexagrams are “non-insta-killed” organisms, and the red circle indicates the destination).

Answer 5

I'm not used to image processing, so this is just an idea:

Implement something like a gaussian filter (blur), which converts each particle (pixels) into a rounded gradient with r = image_size (they all overlap). So you should get an image where the best white pixels should be the best results. Unfortunately, the demonstration in gimp did not succeed because the extreme blur caused the dots to disappear.

Alternatively, you can enlarge all existing pixels by placing all adjacent pixels in the area (example: r = 4), with the left pixels being the same (those with the largest distance to any pixel)