肺实质分割matlab实现

前言

最近有个课程作业，肺实质分割，找了很多代码，大部分都不能用，最后找到一个以matlab为基础的代码，进行修改后，可以实现在matlab中读取单张肺部CT图片进行肺实质分割。完整代码附于最后，使用时，将待分割图片导入到matlab中，然后修改代码中读取的图片名字（改为自己要分割的图片名字）即可。

一、阈值分割

关于阈值分割，找到了多种方式，有全局阈值分割、基于边缘像素分割、基于有效像素分割。从图像上可以看到，肺实质在图像上是巨大的低灰度连通区域，而其他无关结构的灰度值较高，这为阈值法提供了基础，另外，待分割区域灰度值集中且与其他区域相差较大。但这一步我们只是要区别人体和背景，所以用最基础的全局分割就可以了。但为了对比效果，还是选择同时在matlab中保留这三种方法进行对比，代码中gb2gray是matlab内部一种处理图像的函数，通过消除图像色调和饱和度信息同时保留亮度实现将RGB图像或彩色图转换为灰度图像，即灰度化处理的功能，调用这个功能的格式是I = rgb2gray(RGB)，意思是将真彩色图像RGB转换为灰度强度图像。
全局阈值法是一个不断迭代求取最终阈值的方法，其处理过程如下：
（1）为全局阈值设置初始值T，一般而言，若图像中最大的灰度值为Tmax,最小的灰度值为Tmin,那么 ：T=(Tmax+Tmin)/2。
（2）以T为阈值将图像分为前景（灰度值大于等于T）和背景（灰度值小于T），计算前景的平均灰度值TF和背景的平均灰度值TB。
（3）更新阈值T：T=(TF+TB)/2。
（4）重复步骤（2）（3）知道T不在变化为止。
全局阈值、边缘像素分割、有效像素分割对比在matlab中代码如下：

clear; close all; clc
img = rgb2gray(imread('lung1.jpg'));%读取jpg格式 图片 并 转换为灰度图像 
imset=img;
%load imset % 数据集
showflag = 1; % 是否显示中间图像
imidx = 1; % 测试图像序号
 
im = im2double(imset);
 
%% 1.阈值分割
% 全局分割
imbi0 = imbinarize(im);
 
% 基于轮廓像素分割
Eidx = edge(im,'log'); % 拉普拉斯获取边缘
E1 = im(Eidx); % 获取边缘像素
E1 = E1(E1>0);
imbi1 = imbinarize(im,graythresh(E1));
 
% 基于有效像素分割
E2 = im(im>0.02); % 获取非0像素
imbi2 = imbinarize(im,graythresh(E2));

if showflag
    figure
    subplot(221),imshow(im),title('\fontsize{16}原图')
    subplot(222),imshow(imbi0),title('\fontsize{16}全局分割')
    subplot(223),imshow(imbi1),title('\fontsize{16}基于轮廓像素分割')
    subplot(224),imshow(imbi2),title('\fontsize{16}基于有效像素分割')
end

代码运行后的结果：

二、提取人体部分

通过观察不难发现，人体部分在图中所占面积最大（远远大于除人体之外的部分），所以这一步我们求最大连通分量即可。
求最大连通量的一般实现步骤如下：
(1)将输入图像阈值化，若输入图像为彩色图像，则先灰度化后再阈值化为二值图像。
(2)给二值图像添加一个是否已访问的属性，类型为Bool（避免死循环）
(3)找到第一个非零的像素点，将其入栈并将其是否已访问的属性置为真。
(4)以栈的大小是否为0作为结束条件，寻找栈顶元素相邻的八邻域非零像素点，并将它们入栈，结束后将栈顶元素删除。
(5)当栈为空时，表明一个连通区域已经遍历完成，需继续找到下一个非空且未访问过的像素点作为起点，重复4步骤，直到所有的非零像素点都被访问完成。
(6)当所有的连通区域求解完成之后，将像素点个数最大的连通区域标记出来即可。
matlab中有封装完成的bwlabel函数，可以直接调用，实现同样的效果，所以此处直接调用，不在自行封装函数。代码如下：

%% 2.提取人体部分
% 计算连通分量
[label,num] = bwlabel(imbi0);
% 计算最大连通分量
MAX = 0;
for k = 1:num
    maxtmp = sum(find(label==k));
    if maxtmp>MAX
        IDX = k;
        MAX = maxtmp;
    end
end
imbi = label==IDX;
 
if showflag
    figure
    subplot(121),imshow(imbi0),title('\fontsize{16}二值图')
    subplot(122),imshow(imbi),title('\fontsize{16}胸腔')
end

代码运行结果：

三、提取疑似肺质

经过上一步的结果，我们已经得到了分离出的人体部分，即躯干的掩膜，下一步将该掩膜进行孔洞填充获得新的掩膜，两个掩膜相减即可得到肺实质掩膜，其实现代码如下：

%% 3.提取疑似肺质
imbiFull = imfill(imbi,'hole'); % 填充
objtmp = imbiFull-imbi;% 求差
if showflag
    figure
    subplot(131),imshow(imbi),title('\fontsize{16}胸腔')
    subplot(132),imshow(imbiFull),title('\fontsize{16}填充')
    subplot(133),imshow(objtmp),title('\fontsize{16}疑似肺质')
end

代码运行结果：

四、去除非肺质

根据上一步的运行效果图，可以看到其中的非肺质部分已经肉眼不可见，但是为了确保最终生成的结果图达到正确效果，这里还是要进行去除气管等小型非肺质部分操作。操作方法仍可利用连通分量的面积，去除一些小的伪目标。我们可以提出连通区域面积小于P的连通区域即可，我这里设置为P=2000。肺部MASK为去除非肺实质后图像，其代码如下：

%% 4.去除非肺质
P = 2000;%设置连通域面积
MASK = bwareaopen(objtmp,P,4);  % 删除面积小于P的连通分量
if showflag
    figure
    subplot(131),imshow(objtmp),title('\fontsize{16}疑似肺质')
    subplot(132),imshow(MASK),title('\fontsize{16}肺质MASK')
    subplot(133),imshow(im),title('\fontsize{16}原图')
end

代码运行结果：

五、最终输出肺质图

经过上述操作，最终得到了接近理想状态下的肺实质图像，在上一步的操作之后还可以采取一些形态学操作或二值中值滤波等进行各种优化。将前文所述操作生成的有效图像合并输出。其代码如下：

%% 5.最终结果输出
figure
subplot(221),imshow(im),title('\fontsize{16}原图')
subplot(222),imshow(imbi),title('\fontsize{16}胸腔')
subplot(223),imshow(objtmp),title('\fontsize{16}疑似肺质')
subplot(224),imshow(MASK),title('\fontsize{16}肺质MASK')

代码运行结果：

完整代码

% （图像分割）MATLAB胸部肺实质提取
%   by cz  2021.4.20
clear; close all; clc
img = rgb2gray(imread('lung1.jpg'));%读取jpg格式 图片 并 转换为灰度图像 
imset=img;
%load imset % 数据集
showflag = 1; % 是否显示中间图像
imidx = 1; % 测试图像序号
 
im = im2double(imset);
 
%% 1.阈值分割
% 全局分割
imbi0 = imbinarize(im);
 
% 基于轮廓像素分割
Eidx = edge(im,'log'); % 拉普拉斯获取边缘
E1 = im(Eidx); % 获取边缘像素
E1 = E1(E1>0);
imbi1 = imbinarize(im,graythresh(E1));
 
% 基于有效像素分割
E2 = im(im>0.02); % 获取非0像素
imbi2 = imbinarize(im,graythresh(E2));
 
if showflag
    figure
    subplot(221),imshow(im),title('\fontsize{16}原图')
    subplot(222),imshow(imbi0),title('\fontsize{16}全局分割')
    subplot(223),imshow(imbi1),title('\fontsize{16}基于轮廓像素分割')
    subplot(224),imshow(imbi2),title('\fontsize{16}基于有效像素分割')
end
 
 
%% 2.提取人体部分
% 计算连通分量
[label,num] = bwlabel(imbi0);
% 计算最大连通分量
MAX = 0;
for k = 1:num
    maxtmp = sum(find(label==k));
    if maxtmp>MAX
        IDX = k;
        MAX = maxtmp;
    end
end
imbi = label==IDX;
 
if showflag
    figure
    subplot(121),imshow(imbi0),title('\fontsize{16}二值图')
    subplot(122),imshow(imbi),title('\fontsize{16}胸腔')
end
 
 
%% 3.提取疑似肺质
imbiFull = imfill(imbi,'hole'); % 填充
objtmp = imbiFull-imbi;
 
if showflag
    figure
    subplot(131),imshow(imbi),title('\fontsize{16}胸腔')
    subplot(132),imshow(imbiFull),title('\fontsize{16}填充')
    subplot(133),imshow(objtmp),title('\fontsize{16}疑似肺质')
end
 
 
%% 4.去除非肺质
P = 2000;%设置连通域面积
MASK = bwareaopen(objtmp,P,4);  % 删除面积小于P的连通分量
if showflag
    figure
    subplot(131),imshow(objtmp),title('\fontsize{16}疑似肺质')
    subplot(132),imshow(MASK),title('\fontsize{16}肺质MASK')
    subplot(133),imshow(im),title('\fontsize{16}原图')
end
 
 
%% 5.最终结果输出
 
figure
subplot(221),imshow(im),title('\fontsize{16}原图')
subplot(222),imshow(imbi),title('\fontsize{16}胸腔')
subplot(223),imshow(objtmp),title('\fontsize{16}疑似肺质')
subplot(224),imshow(MASK),title('\fontsize{16}肺质MASK')