matlab保留三维矩阵的最大连通域

对标matlab官方函数bwlabel

自写，封装成函数，拿来就用。

function [res] = cwlabel(image)
%CWLABEL 此处显示有关此函数的摘要
%   提取三维空间的最大连通域
    [w, l, h] = size(image);
    res = image;
    value = [];
    step = [[-1, -1, -1]; [-1, -1, 0]; [-1, -1, 1]; [-1, 0, -1]; [-1, 0, 0]; [-1, 0, 1]; [-1, 1, -1]; [-1, 1, 0]; [-1, 1, 1];...
[0, -1, -1]; [0, -1, 0]; [0, -1, 1]; [0, 0, -1]; [0, 0, 1]; [0, 1, -1]; [0, 1, 0]; [0, 1, 1];...
[1, -1, -1]; [1, -1, 0]; [1, -1, 1]; [1, 0, -1]; [1, 0, 0]; [1, 0, 1]; [1, 1, -1]; [1, 1, 0]; [1, 1, 1]];
    step = step';
    for i = 1 : w
        for j = 1 : l
            for k = 1 : h
                if image(i, j, k) ~= 0 && res(i, j, k) ~= 0
                    q = [];
                    q = [q; [i, j, k]];
                    res(i, j, k) = 0;
                    num = 0;
                    while ~isempty(q)
                        cur_x = q(1, 1);
                        cur_y = q(1, 2);
                        cur_z = q(1, 3);
                        q(1, :) = [];
                        num = num + 1;
                        for s = step
                            new_x = cur_x + s(1);
                            new_y = cur_y + s(2);
                            new_z = cur_z + s(3);
                            if new_x < 1 || new_x > w || new_y < 1 || new_y > l || new_z < 1 || new_z > h
                                continue;
                            end
                            if res(new_x, new_y, new_z) ~= 0 && image(new_x, new_y, new_z) ~= 0
                                res(new_x, new_y, new_z) = 0;
                                q = [q; [new_x, new_y, new_z]];
                            end
                        end
                    end
                    value = [value; [i, j, k, num]];
                end
            end
        end
    end
    
    max_value = 0;
    for i = 1 : length(value)
        if value(i, 4) > max_value
            max_value = value(i, 4);
            begin_node = [value(i, 1), value(i, 2), value(i, 3)];
        end
    end
    
    res = image;
    res(res > 0) = 0;
    vis = image;
    vis(begin_node(1), begin_node(2), begin_node(3)) = 0;
    res(begin_node(1), begin_node(2), begin_node(3)) = image(begin_node(1), begin_node(2), begin_node(3));
    q = [];
    q = [q; begin_node];
    while ~isempty(q)
        cur_x = q(1, 1);
        cur_y = q(1, 2);
        cur_z = q(1, 3);
        q(1, :) = [];
        for s = step
            new_x = cur_x + s(1);
            new_y = cur_y + s(2);
            new_z = cur_z + s(3);
            if new_x < 1 || new_x > w || new_y < 1 || new_y > l || new_z < 1 || new_z > h
                continue;
            end
            if vis(new_x, new_y, new_z) ~= 0 && image(new_x, new_y, new_z) ~= 0
                vis(new_x, new_y, new_z) = 0;
                res(new_x, new_y, new_z) = image(new_x, new_y, new_z);
                q = [q; [new_x, new_y, new_z]];
            end
        end
    end
    
end