logo标签列表

在具有分布式粒子的图像的空闲区域内拟合最大圆形

matlabmathimage-processinggeometryparticles
52

我正在处理图片，以便在包含分布颗粒的图像的任何自由区域中检测和适配最大可能的圆形：1

(能够检测出颗粒的位置).

一种方法是定义一个接触三个点组合的圆形，检查圆形是否为空，然后在所有空圆形中找到最大的圆。但是，这会导致组合数量巨大，即 C(n,3)，其中n是图像中颗粒的总数。

如果有人能够给我提供任何提示或替代方法，我将不胜感激。

2
@AnderBiguri 我相信他的意思是最大的圆 没有 任何粒子在里面 - Daniel Alder
3
我认为他的意思是，最大的圆可以适应这个空间，而不包含任何点... - Wolfie
@Wolfie 呃...我误解了“3点”的说法。 - Ander Biguri
3
@AnderBiguri，我认为所提出的算法是：1）获取所有点的三元组。2）使用这些三点集定义圆（这将是唯一的）。3）丢弃包含点的圆。4）找到剩余圆中最大的圆。问题在于有许多三元组，因此该算法速度较慢。 - Wolfie
1
@AnderBiguri 非常感谢您提供的解决方案。它对我的问题起到了很好的作用。 - T50740
显示剩余6条评论
5个回答
91

让我们做些数学吧，因为数学总会到达终点！

维基百科：

在数学中，沃罗诺伊图是根据平面上特定子集中的点到距离进行划分的区域。

例如：

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


voronoi(x,y);

输入图像描述

这张图的好处在于，如果你注意到，所有蓝色区域的边缘/顶点都与它们周围的点等距。因此，如果我们知道顶点的位置，并计算到最近点的距离，那么我们可以选择具有最大距离的顶点作为圆的中心。

有趣的是，沃罗诺伊区域的边缘也被定义为由德劳内三角剖分生成的三角形的外心。

因此，如果我们计算出该区域的德劳内三角剖分及其外心，则...

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

然后我们有了所有可能没有任何点在其中的圆的中心 (cc) 和半径 (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]区间内找到圆心。您可以轻松修改它以更改此限制。您也可以尝试找到完全适合所需区域内的圆（而不仅仅是中心）。这需要在最后添加一个小额外操作来获取最大值。

1
@GuyRT 的确。我认为这种方法的优点在于可以在连续空间中工作。结果不需要是整数像素索引。 - Ander Biguri
1
@fyrepenguin 如果你仔细阅读，这个答案的结尾有一个关于此事的免责声明。圆的中心不能越界（有一个明确的 if 条件），但是圆可以部分越界。如果需要，很容易修复。 - Ander Biguri
1
@AnderBiguri 我猜我没有仔细阅读，我的错。然而，我正在测试它，并且最终得到了一个以x轴上大约-11为中心的圆。很不幸，我现在身边没有带有MATLAB的电脑，所以无法向您展示我所说的内容的图像。 - fyrepenguin
1
@fyrepenguin 这种情况之前确实发生过，所以我才加了 if。也许我漏掉了什么。如果可以的话，明天再试一次。 - Ander Biguri
1
@AnderBiguri 我现在看到了那个if语句，但我不确定为什么它对我来说表现得有些奇怪，我明天也可以再检查一下。这确实是一个非常聪明的问题和答案。 - fyrepenguin
显示剩余10条评论
25
我想提出另一种基于网格搜索和细化的解决方案。它不像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%的点，并仅考虑包含最大距离的区域（表示保留值的曲面片段）：

保留最大区域后

最后：

展示找到的圆

2
嗯，你的图比我的好看！ - Ander Biguri
1
@AnderBiguri 不用担心...只需要看编辑之前的版本 :) - Dev-iL
3
除了实际应用之外，情节看起来很棒！ - moonwave99
14
你可以使用图像处理工具箱中的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');

enter image description here

2
你应该使用和我相同的数据，这样我们才能更好地进行比较！bwdist的结果看起来非常奇怪。 - Ander Biguri
4
你知道在 imread 中可以使用 URL 吗？例如这样可以正常工作：img = imread('https://i.stack.imgur.com/AbmxL.jpg');。 - Dev-iL
1
@AnderBiguri 我使用了问题中的数据。由于任务是图像处理，我认为你应该将你的方法应用在图像上。 - rahnema1
@Dev-iL 我无法查看或直接下载附加的图片！ - rahnema1
@rahnema1 对不起，我没有意识到。非常好的答案！ - Ander Biguri
14

事实上，这个问题可以使用“直接搜索”方法来解决（正如另一个答案所示），这意味着我们可以将其视为一个全局优化问题。有各种方法可以解决这样的问题，每种方法适用于特定情况。出于个人好奇心，我决定使用遗传算法来解决此问题。

一般来说，这样的算法要求我们将解决方案看作是一组“基因”，在某种“适应性函数”的作用下进行“进化”。恰巧，在这个问题中很容易识别基因和适应性函数：

  • 基因：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代的“时间流逝”图：

时间流逝

(蓝色点代表当前的一代，红叉代表被“即刻杀死”的生物，绿色六芒星代表“非即刻杀死”的生物，红色圆圈标记着目的地)。

3
虽然这种方法对读者来说似乎只是一种练习（意思是可能最好使用其他方法，包括你的方法），但天哪，这太酷了！那个 gif :P - Ander Biguri
1

我不擅长图像处理，这只是一个想法：

实现类似高斯滤波器（模糊）的东西，将每个粒子（像素）转换为带有r = image_size的圆形渐变（所有重叠）。这样，您应该可以得到一张图片，其中最白的像素应该是最好的结果。不幸的是，Gimp中的演示失败了，因为极端模糊使点消失了。

或者，您可以增量地通过标记区域内所有相邻像素（例如：r = 4）来扩展所有现有像素，剩下的像素将是相同的结果（与任何像素的距离最大的像素）。

1
这是一个不错的开始，但最终图像处理会受到像素大小的限制。始终会存在1*像素大小的最小误差。 - Ander Biguri
2
无论如何，这种方法与@AnderBiguri的回答相比是不可比较的 :-P - Daniel Alder
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接