在MATLAB中找到具有共同重叠区域的多个圆

Posted 2021-04-07

我有两个矩阵。每个矩阵的维数为3 * k，表示k个圆，每个列的形式为[x y r]，其中（x，y）是圆的中心，r是半径。所以一个有3个圆的矩阵是[x1 x2 x3; y1 y2 y3; r1 r2 r3]

我需要找到重叠的圆圈和两个矩阵之间的重叠区域。假设，考虑第一个矩阵中的圆。现在考虑第二个矩阵中的每个其他圆。现在，我需要找出第二个矩阵中的哪些圆与所考虑的圆重叠。我需要为第一个矩阵中的每个圆圈执行此操作。类似地，对于第一矩阵，第二矩阵中的每个圆。

所以我需要两个矩阵中的每个圆（其中有k1 + k2），与另一个矩阵重叠多少区域以及另一个矩阵中重叠的圆圈是什么。

显然，可能存在多个重叠的圆。对于两个矩阵，维度k可以是不同的。我有每个矩阵2个额外的矩阵，一个按x的坐标排序，另一个按y坐标排序，如果这有助于计算。问题是矩阵中有很多圆圈，我正在寻找一种有效的方法做到这一点。此外，我想扩展这个超过两个矩阵，然后一个有效的算法将大大改善执行时间。

这个链接给出了两个图像（我对应的圆圈的相应矩阵）的预览：https://www.dropbox.com/s/om5has5uw91dj6p/overlap.jpg

答案

你需要两件事：

1. 找出哪些圆圈重叠（距离<半径之和）
2. 计算这些对的重叠区域

EDITED我修改了代码以添加一些矢量化并捕获一个圆圈完全位于另一个圆圈内的情况 - 其中Wolfram formula正在分解。与1000 x 2000圈相比，显示速度。

以下是一些执行这些操作的代码（使用简单示例）。它似乎相当有效。花了1.4秒来比较1000与2000圈（166万重叠）。

% circles [x; y; r]
m1 = rand(3, 1000);
m2 = rand(3, 2000);

tic
% distances between circles:
dx = bsxfun(@minus, m1(1,:), m2(1,:)');
dy = bsxfun(@minus, m1(2,:), m2(2,:)');
sumr = bsxfun(@plus, m1(3,:), m2(3,:)');

% distance between centers:
dist = sqrt(dx.^2 + dy.^2);

% pairs that have overlap:
pairs = find(dist < sumr);

s1 = size(m1,2);
s2 = size(m2,2);

% calculate overlaps
c1 = repmat(1:s1, [s2 1]); % circle from m1 corresponding to value in pairs()
c2 = repmat(1:s2, [s1 1])'; % ditto for m2

o = overlap(m1(3, c1(pairs)), m2(3, c2(pairs)), dist(pairs)');
toc

将以下内容放在路径中自己的文件overlap.m中：这会计算重叠

function o = overlap(r1, r2, d)
r12 = r1.^2;
r22 = r2.^2;
d2 = d.^2;
a1 = acos((d2 + r12 - r22)./(2.*d.*r1));
a2 = acos((d2 + r22 - r12)./(2.*d.*r2));
o = r12.*a1 + r22.*a2-.5*sqrt((-d+r1+r2)*(d+r1-r2)*(d-r1+r2)*(d+r1+r2));

% fix situation where r1 - r2 > d:
% i.e. one circle fully inside the other
f = find(abs(imag(a1)) + abs(imag(a2)) > 0);
o(f) = pi * min(r1(f),r2(f)).^2;

另一答案

我想过使用rasterization来近似圆圈之间的重叠。想法是创建N点以近似圆形路径，然后将它们传递给poly2mask函数。这将返回在定义圆的位置填充的二进制掩码图像（算法应以子像素精度工作）。对所有圆对进行此操作，我们可以通过仅将逻辑AND应用于两个掩码来计算交集，并计算结果中剩余的“开”像素的数量（或者更好地使用bwarea进行稍微更好的估计）。这又是一个近似值，因为我们使用的是离散的像素网格......

从好的方面来说，我们不必担心数学方程式。实际上，这种方法适用于任何多边形形状，在这种形状下很难为交叉点提出封闭形式的解决方案。

不幸的是，这种蛮力方法被证明比我预期的要慢（没有像@Floris的代码那么快），我仍在发布我的尝试。我借用了检查两个圆圈是否对其他代码感兴趣的想法，如果您拥有的大多数圆圈相距很远，这应该会显着缩短时间。

只是为了好玩，我添加了一些代码来可视化过程！

ANIMATION = true;

% size of image we are working in
%sz = [180 360];
sz = [100 100];

% generated two sets of random circles (specified as columns [x;y;r])
% (I am just trying to create circles wholly visible within the box)
k1 = 100; k2 = 80;
M1 = bsxfun(@times, rand(3,k1), [sz(:)*0.6+sz(:)*0.2;0.1*min(sz)+0.1*min(sz)]);
M2 = bsxfun(@times, rand(3,k2), [sz(:)*0.6+sz(:)*0.2;0.1*min(sz)+0.1*min(sz)]);

% animation
if ANIMATION
    clf
    hImg = imshow(zeros(sz), 'InitialMag','fit');
    hLine(1) = line(NaN, NaN, 'Color','r', 'LineWidth',2);
    hLine(2) = line(NaN, NaN, 'Color','b', 'LineWidth',2);
    axis on
end

% used to approximate circles
num = 50;
t = linspace(0, 2*pi, num);
ct = cos(t); st = sin(t);

% test to find which circles intersect
dist = pdist2(M1(1:2,:)', M2(1:2,:)');
sumr = bsxfun(@plus, M1(3,:)', M2(3,:));
skipIdx = (sumr < dist);

% compute overlap between all pairs
overlap = zeros(k1,k2);
for i=1:k1
    for j=1:k2
        % early skip if circles dont interset
        if skipIdx(i,j), continue, end

        % compute each circle points
        x1 = M1(3,i)*ct + M1(1,i); y1 = M1(3,i)*st + M1(2,i);
        x2 = M2(3,j)*ct + M2(1,j); y2 = M2(3,j)*st + M2(2,j);

        % rasterize circles
        BW1 = poly2mask(x1, y1, sz(1), sz(2));
        BW2 = poly2mask(x2, y2, sz(1), sz(2));

        % compute area of intersection of the two masks
        %overlap(i,j) = sum(BW1(:)&BW2(:));
        overlap(i,j) = bwarea(BW1&BW2);

        if ANIMATION
            % update animation
            set(hImg, 'CData',BW1&BW2)
            set(hLine(1), 'XData',x1, 'YData',y1)
            set(hLine(2), 'XData',x2, 'YData',y2)
            title(sprintf('%g',overlap(i,j)))
            drawnow
        end
    end
end

---

为了了解近似的精确程度，我捕获了一个迭代，其中一个圆圈完全位于另一个圆圈内。这样我们就可以将实际交叉区域与我们的近似值进行比较。

对于下面显示的图像，我们有：

K>> M1(:,i)
ans =
   58.2106    % x
   24.0996    % y
    8.9387    % r

K>> pi*M1(3,i)^2
ans =
  251.0122

K>> overlap(i,j)
ans =
  252.5000

我认为这足够接近:)使用更高分辨率的像素网格将改善近似值，但肯定会减慢它的速度。