用分布式粒子拟合图像中自由区域的最大圆
Fitting largest circle in free area in image with distributed particle
我正在处理图像以检测并拟合包含分布式粒子的图像的任何自由区域中的最大可能圆：
（能够检测粒子的位置）。
一个方向是定义一个与任意三点组合相接的圆，检查圆是否为空，然后在所有空圆中找出最大的圆。然而，它会导致大量的组合，即 C(n,3)，其中 n 是图像中的粒子总数。
如果有人能提供任何提示或我可以探索的替代方法，我将不胜感激。
我不习惯图像处理，所以这只是一个想法：
实现类似于高斯滤波器（模糊）的东西，它将每个粒子（像素）转换为具有 r=image_size 的圆形渐变（它们全部重叠）。这样，您应该得到一张图片，其中最白的像素应该是最好的结果。不幸的是，gimp中的演示失败了，因为极度模糊使点消失了。
或者，您可以通过标记一个区域中的所有相邻像素来逐步扩展所有现有像素（例如：r=4），剩下的像素将是相同的结果（与任何像素距离最大的像素）
让我们做一些数学，我的朋友，因为数学总会走到尽头！
维基百科：
In mathematics, a Voronoi diagram is a partitioning of a plane into
regions based on distance to points in a specific subset of the plane.
例如：
rng(1)
x=rand(1,100)*5;
y=rand(1,100)*5;
voronoi(x,y);
这张图的好处是，如果你注意到的话，那些蓝色区域的所有 edges/vertices 都与它们周围的点距离相等。因此，如果我们知道顶点的位置，并计算到最近点的距离，那么我们可以选择距离最大的顶点作为我们的圆心。
有趣的是，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
然后我们得到所有内部没有点的可能圆的圆心 (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] 内的中心。您可以轻松修改它以更改此约束。您还可以尝试在感兴趣的区域内找到完全适合的圆圈（而不是仅仅位于中心）。这将需要在获得最大值的末尾进行少量添加。
您可以使用图像处理工具箱中的 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');
我想提出另一种解决方案，该解决方案基于经过优化的网格搜索。它不像 Ander 的那样先进，也不像 rahnema1 的那样简短，但应该非常容易理解和理解。而且，它运行得相当快。
算法包含几个阶段：
- 我们生成一个均匀分布的网格。
- 我们找到网格中的点到所有提供的点的最小距离。
- 我们丢弃所有距离低于某个百分位数（例如第 95 个）的点。
- 我们选择包含最大距离的区域（如果我的初始网格足够好，这应该包含正确的中心）。
- 我们在所选区域周围创建一个新的网格并再次查找距离（这部分显然不是最优的，因为计算了所有点的距离，包括远点和不相关的点）。
- 我们在区域内迭代细化，同时关注前 5% 的值的方差 -> 如果它低于某个预设阈值，我们就会中断。
几个注意事项：
- 我假设圆圈不能超出散点的范围（即散点的边界正方形充当"invisible wall"）。
- 合适的百分位数取决于初始网格的精细程度。这也会影响
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%的距离点，只考虑包含最大距离的区域（那块面代表保留值）：
最后：
事实上这个问题可以用 "direct search" 来解决（在 ) means one can look at this as a global optimization 问题中可以看出。有多种方法可以解决此类问题，每种方法都适用于特定的场景。出于出于个人好奇，我决定使用 遗传算法 .
来解决这个问题
一般来说，这样的算法需要我们把解看成是"genes"服从"evolution"在某个"fitness function"下的集合。碰巧的是，很容易识别这个问题中的基因和适应度函数：
- 基因：
x，y，r。
- 适应度函数：从技术上讲，圆的最大面积，但这相当于最大
r（或最小-r，因为算法需要一个函数来最小化).
- 特殊约束 - 如果
r 大于提供的最近点（即圆包含一个点）的欧式距离，则有机体 "dies".
下面是这样一个算法的基本实现（basic”因为完全没有优化，还有很大的优化空间在这个问题中没有双关语）。
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 代 运行 的 "time-lapse" 图：
（其中蓝点为当前世代，红叉为"insta-killed"种生物，绿色六芒星为"non-insta-killed"种生物，红圈为归宿）
我正在处理图像以检测并拟合包含分布式粒子的图像的任何自由区域中的最大可能圆：
（能够检测粒子的位置）。
一个方向是定义一个与任意三点组合相接的圆，检查圆是否为空，然后在所有空圆中找出最大的圆。然而，它会导致大量的组合，即 C(n,3)，其中 n 是图像中的粒子总数。
如果有人能提供任何提示或我可以探索的替代方法，我将不胜感激。
我不习惯图像处理，所以这只是一个想法：
实现类似于高斯滤波器（模糊）的东西，它将每个粒子（像素）转换为具有 r=image_size 的圆形渐变（它们全部重叠）。这样，您应该得到一张图片，其中最白的像素应该是最好的结果。不幸的是，gimp中的演示失败了，因为极度模糊使点消失了。
或者，您可以通过标记一个区域中的所有相邻像素来逐步扩展所有现有像素（例如：r=4），剩下的像素将是相同的结果（与任何像素距离最大的像素）
让我们做一些数学，我的朋友，因为数学总会走到尽头！
维基百科：
In mathematics, a Voronoi diagram is a partitioning of a plane into regions based on distance to points in a specific subset of the plane.
例如：
rng(1)
x=rand(1,100)*5;
y=rand(1,100)*5;
voronoi(x,y);
这张图的好处是，如果你注意到的话，那些蓝色区域的所有 edges/vertices 都与它们周围的点距离相等。因此，如果我们知道顶点的位置，并计算到最近点的距离，那么我们可以选择距离最大的顶点作为我们的圆心。
有趣的是，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
然后我们得到所有内部没有点的可能圆的圆心 (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] 内的中心。您可以轻松修改它以更改此约束。您还可以尝试在感兴趣的区域内找到完全适合的圆圈（而不是仅仅位于中心）。这将需要在获得最大值的末尾进行少量添加。
您可以使用图像处理工具箱中的 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');
我想提出另一种解决方案，该解决方案基于经过优化的网格搜索。它不像 Ander 的那样先进，也不像 rahnema1 的那样简短，但应该非常容易理解和理解。而且，它运行得相当快。
算法包含几个阶段：
- 我们生成一个均匀分布的网格。
- 我们找到网格中的点到所有提供的点的最小距离。
- 我们丢弃所有距离低于某个百分位数（例如第 95 个）的点。
- 我们选择包含最大距离的区域（如果我的初始网格足够好，这应该包含正确的中心）。
- 我们在所选区域周围创建一个新的网格并再次查找距离（这部分显然不是最优的，因为计算了所有点的距离，包括远点和不相关的点）。
- 我们在区域内迭代细化，同时关注前 5% 的值的方差 -> 如果它低于某个预设阈值，我们就会中断。
几个注意事项：
- 我假设圆圈不能超出散点的范围（即散点的边界正方形充当"invisible wall"）。
- 合适的百分位数取决于初始网格的精细程度。这也会影响
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%的距离点，只考虑包含最大距离的区域（那块面代表保留值）：
最后：
事实上这个问题可以用 "direct search" 来解决（在
一般来说，这样的算法需要我们把解看成是"genes"服从"evolution"在某个"fitness function"下的集合。碰巧的是，很容易识别这个问题中的基因和适应度函数：
- 基因：
x，y，r。 - 适应度函数：从技术上讲，圆的最大面积，但这相当于最大
r（或最小-r，因为算法需要一个函数来最小化). - 特殊约束 - 如果
r大于提供的最近点（即圆包含一个点）的欧式距离，则有机体 "dies".
下面是这样一个算法的基本实现（basic”因为完全没有优化，还有很大的优化空间在这个问题中没有双关语）。
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 代 运行 的 "time-lapse" 图：
（其中蓝点为当前世代，红叉为"insta-killed"种生物，绿色六芒星为"non-insta-killed"种生物，红圈为归宿）