这是图像处理中的常见问题。有很多变化,例如洪水填充图像中的某个区域,或者查找哪些像素属于同一区域。一种常见的方法是使用深度优先搜索。这个想法是,从左到右、从上到下遍历图像,对于遇到的任何等于 1 的像素,将它们添加到堆栈中。对于堆栈中的每个像素,您从堆栈中弹出,然后查看该像素周围的相邻像素。添加到堆栈中的任何为 1 的像素。您需要保留一个附加变量,其中任何像素您已经访问过,您不将它们添加到堆栈中。当堆栈为空时,我们发现那些像素是整个区域,因此您可以使用唯一的 ID 来标记它们。然后重复此过程,直到用完图像中的区域。
因此,假设您的矩阵存储在A
,这是基本算法:
初始化一个大小相同的数组A
那是logical
。这将记录我们检查或访问过哪些像素。还初始化输出数组B
全部为零,为您提供您正在寻找的所有连接组件。任何最终为零的位置都不属于任何连接的组件。还要初始化一个 ID 计数器,用于跟踪每个连接组件的标签。
-
对于矩阵中的每个位置:
A。如果位置是0
,将此位置标记为已访问并继续。
b.如果我们已经访问过该位置,则继续。
C。如果我们还没有访问过该位置...请转到步骤#3。
-
将此未访问的位置添加到堆栈中。
A。虽然这个堆栈不为空...
b.从堆栈中弹出此位置
C。如果我们访问过该位置,则继续。
d.否则,将此位置标记为已访问,并使用连接的组件 ID 标记此位置。
e.给定这个位置,查看 8 个相邻像素。
F。删除此列表中那些已访问过、不等于 1 或超出矩阵范围的像素
G。无论剩余位置如何,都将它们添加到堆栈中。
一旦堆栈为空,就增加计数器,然后返回到步骤 #2。
继续下去,直到我们访问完数组中的所有位置。
话不多说,这是代码。
%// Step #1
visited = false(size(A));
[rows,cols] = size(A);
B = zeros(rows,cols);
ID_counter = 1;
%// Step 2
%// For each location in your matrix...
for row = 1 : rows
for col = 1 : cols
%// Step 2a
%// If this location is not 1, mark as visited and continue
if A(row,col) == 0
visited(row,col) = true;
%// Step 2b
%// If we have visited, then continue
elseif visited(row,col)
continue;
%// Step 2c
%// Else...
else
%// Step 3
%// Initialize your stack with this location
stack = [row col];
%// Step 3a
%// While your stack isn't empty...
while ~isempty(stack)
%// Step 3b
%// Pop off the stack
loc = stack(1,:);
stack(1,:) = [];
%// Step 3c
%// If we have visited this location, continue
if visited(loc(1),loc(2))
continue;
end
%// Step 3d
%// Mark location as true and mark this location to be
%// its unique ID
visited(loc(1),loc(2)) = true;
B(loc(1),loc(2)) = ID_counter;
%// Step 3e
%// Look at the 8 neighbouring locations
[locs_y, locs_x] = meshgrid(loc(2)-1:loc(2)+1, loc(1)-1:loc(1)+1);
locs_y = locs_y(:);
locs_x = locs_x(:);
%%%% USE BELOW IF YOU WANT 4-CONNECTEDNESS
% See bottom of answer for explanation
%// Look at the 4 neighbouring locations
% locs_y = [loc(2)-1; loc(2)+1; loc(2); loc(2)];
% locs_x = [loc(1); loc(1); loc(1)-1; loc(1)+1];
%// Get rid of those locations out of bounds
out_of_bounds = locs_x < 1 | locs_x > rows | locs_y < 1 | locs_y > cols;
locs_y(out_of_bounds) = [];
locs_x(out_of_bounds) = [];
%// Step 3f
%// Get rid of those locations already visited
is_visited = visited(sub2ind([rows cols], locs_x, locs_y));
locs_y(is_visited) = [];
locs_x(is_visited) = [];
%// Get rid of those locations that are zero.
is_1 = A(sub2ind([rows cols], locs_x, locs_y));
locs_y(~is_1) = [];
locs_x(~is_1) = [];
%// Step 3g
%// Add remaining locations to the stack
stack = [stack; [locs_x locs_y]];
end
%// Step 4
%// Increment counter once complete region has been examined
ID_counter = ID_counter + 1;
end
end %// Step 5
end
通过你的示例矩阵,这就是我得到的B
:
B =
1 1 0 0 0 0 0
1 1 0 0 2 2 0
1 1 0 0 0 2 0
1 1 0 0 0 0 0
0 0 0 0 0 3 0
0 0 0 0 0 0 0
在 4 个连通邻域中搜索
修改代码以在 4 个连通区域中搜索,即只有北、东、西和南,您看到的部分%// Look at the 8 neighbouring locations
, 那是:
%// Look at the 8 neighbouring locations
[locs_y, locs_x] = meshgrid(loc(2)-1:loc(2)+1, loc(1)-1:loc(1)+1);
locs_y = locs_y(:);
locs_x = locs_x(:);
要以 4 连接方式进行搜索,您只需修改此代码以仅给出这些基本方向:
%// Look at the 4 neighbouring locations
locs_y = [loc(2)-1; loc(2)+1; loc(2); loc(2)];
locs_x = [loc(1); loc(1); loc(1)-1; loc(1)+1];
其余代码保持不变。
为了匹配 MATLAB 的bwlabel
功能
如果你想匹配 MATLAB 的输出bwlabel功能,bwlabel
按列主序或 FORTRAN 顺序搜索连通分量。上面的代码按行专业或 C 顺序搜索。因此,您只需首先沿着列而不是行进行搜索,就像上面的代码所做的那样,您可以通过交换两者的顺序来完成此操作for
循环。
具体来说,不要这样做:
for row = 1 : rows
for col = 1 : cols
....
....
你会这样做:
for col = 1 : cols
for row = 1 : rows
....
....
现在应该复制以下输出bwlabel
.