使用分布式粒子拟合图像中自由区域中的最大圆

Question

我正在处理图像,以检测并适合包含分布式粒子的图像的任何自由区域中的最大圆圈:

(能够检测粒子的位置).

一个方向是定义一个圆圈触摸任何三点组合,检查圆圈是否为空,然后在所有空圆圈中找到最大的圆圈.然而,它导致大量组合,即C(n,3),n图像中的粒子总数在哪里.

如果有人能提供我可以探索的任何提示或替代方法,我将不胜感激.

Answer 1

让我的朋友做一些数学,因为数学将永远到最后!

维基百科:

在数学中,Voronoi图是基于到平面的特定子集中的点的距离将平面划分成区域.

例如:

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


voronoi(x,y);

关于这个图的好处是,如果你注意到,那些蓝色区域的所有边/顶点都与它们周围的点的距离相等.因此,如果我们知道顶点的位置,并计算到最近点的距离,那么我们可以选择具有最高距离的顶点作为圆的中心.

有趣的是,Voronoi区域的边缘也被定义为Delaunay三角剖分产生的三角形的外心.

因此,如果我们计算该区域的Delaunay三角剖分及其外心

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

并计算外心和定义每个三角形的任何点之间的距离:

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

然后我们得到所有可能的圆圈的center(cc)和radius(distance),这些圆圈内没有任何点.我们只需要最大的一个!

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

现在让我们绘制

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

注意圆的中心如何在Voronoi图的一个顶点上.

注意:这将在[0-5],[0-5]内找到中心.您可以轻松修改它以更改此约束.您还可以尝试在感兴趣的区域内找到适合其整体的圆圈(而不仅仅是中心).这将需要在最终获得最大值的情况下进行少量添加.

Answer 2

我想提出另一种基于网格搜索和精化的解决方案.它不像Ander那样先进,也不像rahnema1那么短,但它应该很容易理解.而且,它运行得非常快.

该算法包含几个阶段:

1. 我们生成一个均匀间隔的网格.
2. 我们发现网格中点的最小距离为所有提供的点.
3. 我们丢弃距离低于某个百分位数(例如第95位)的所有点.
4. 我们选择包含最大距离的区域(如果我的初始网格足够好,则应该包含正确的中心).
5. 我们在所选区域周围创建一个新的网格网格并再次找到距离(这部分显然是次优的,因为距离计算到所有点,包括远距离和不相关的距离).
6. 我们在区域内迭代细化,同时密切关注前5%值的方差 - >如果它低于某个预设阈值我们就会中断.

几点说明:

• 我假设圆不能超出散点的范围(即散射的边界作为"看不见的墙").
• 适当的百分位数取决于初始网格的精细程度.这也会影响while迭代次数和最佳初始值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');

我得到Ander的示例数据的结果是[x,y,r] = [0.7832, 2.0694, 0.7815](它是相同的).执行时间约为Ander解决方案的一半.

以下是中间图:

从一个点到所有提供点的集合的最大(清晰)距离的轮廓:

在考虑距离边界的距离之后,仅保留距离点的前5%,并且仅考虑包含最大距离的区域(该表面代表保持的值):

最后:

Answer 3

您可以使用Image Processing Toolbox中的bwdist来计算图像的距离变换.这可以被视为创建voronoi图的方法,在@ AnderBiguri的答案中有很好的解释.

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
[x y] = find(dist == radius);
imshow(dist,[]);
hold on
plot(y,x,'ro');

Answer 4

使用"直接搜索"可以解决这个问题的事实(可以在另一个答案中看到)意味着可以将其视为全局优化问题.存在各种方法来解决这些问题,每种方法都适用于某些场景.出于个人的好奇心,我决定使用遗传算法来解决这个问题.

一般来说,这种算法要求我们将解决方案视为一组在某种"适应度函数"下经历"进化"的"基因".碰巧的是,在这个问题中识别基因和适应度函数非常容易:

• 基因:x,y,r.
• 适应度函数:从技术上讲,圆的最大面积,但这相当于最大值r(或最小值-r,因为算法需要一个函数来最小化).
• 特殊约束 - 如果r大于到最近提供点的欧氏距离(即圆圈包含一个点),则生物体"死亡".

下面是这种算法的基本实现(" 基本 ",因为它完全没有优化,并且有很多优化空间,^{没有双关语意图}在这个问题上).

function [x,y,r] = q42806059b(cloudOfPoints)
% Problem setup
if nargin == 0
  rng(1)
  cloudOfPoints = rand(100,2)*5; % equivalent to Ander's 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(S-L)], [], options);
x = g(1); y = g(2); r = g(3);
%{
plot(x,y,'ro'); plot(x,y,'r*'); 
rectangle('Position',[x-r,y-r,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
%}

这是典型运行的47代的"延时"情节:

(蓝点是当前一代,红十字是"insta-killed"生物,绿卦是"非insta-killed"生物,红圈标志着目的地).