基于Matlab的3-D胸部扫描CT切片的肺部分割——使用活动轮廓(snakes)进行三维分割及建模

最新推荐文章于 2024-07-18 10:54:00 发布

亦我飞也

最新推荐文章于 2024-07-18 10:54:00 发布

阅读量4k

点赞数 11

分类专栏： Matlab 医学图像分割图像建模文章标签： Matlab 肺部分割建模活动轮廓（snakes）

本文链接：https://blog.csdn.net/qq_32809093/article/details/112429365

版权

医学图像分割同时被 3 个专栏收录

34 篇文章 30 订阅

订阅专栏

Matlab

12 篇文章 10 订阅

订阅专栏

图像建模

5 篇文章 3 订阅

订阅专栏

目录：

一、准备数据

二、肺部分割

三、创建种子MASK掩膜并使用活动轮廓（snakes）分割肺部

四、计算分割肺的体积

例程完整源码：

参考链接

本例程配套完整源码和数据下载

此示例显示了如何使用活动轮廓(snakes)执行三维分割。您可以使用the Volume Viewer app查看结果

一、准备数据

将人体胸部CT扫描数据加载到工作空间中。要运行此示例，您必须使用附加浏览器从MathWorks下载样本数据。

使用加载项资源管理器安装示例数据

图像处理工具箱使样本3-D MRI (三维磁共振)胸部扫描数据集作为一个可选功能。要获取此数据集，请使用附加组件浏览下载它

1、打开Matalb——>主页——>选择附加功能 ——> 获取附加功能

2、在打开的附加功能资源管理器 中搜索 Image Processing Toolbox Image Data

3、在搜索结果中单击数据包——>安装。遵循安装人员提供的说明

4、下载完成后得到的安装文件夹MathWorks

5、在Windows平台进行安装：

在Matlab的安装目录下选中目录——>键盘cmd——>回车

输入SupportSoftwareInstaller.exe -archives C:\Users\Administrator\Downloads\MathWorks\SupportPackages\R2020b\

（注Matlab2017----Matlab2020版本可能是install_supportsoftware.exe）还有后面刚下载安装文件夹目录自行更改

勾选——>下一步——>我接受——>下一步——>安装完毕

6、关闭Matlab后重新打开，命令行窗口输入

load chestVolume
whos

加载数据到Matlab

7、将CT扫描数据从int16转换为single，以将值标准化到范围[0，1]

V = im2single(V);

8、使用matlab中 APP——> 图像处理与计算机视觉——>Volume Viewer查看胸部扫描。也可以在命令行窗口使用volumeViewer(V)命令打开应用程序。使用ct-骨骼(ct-bone)，获得胸部扫描的最佳视图

volumeViewer(V)

二、肺部分割

使用活动轮廓(active contour)技术对CT扫描数据进行肺部分割。主动轮廓是一种区域生长算法，需要初始种子点(initial seed points)。该示例使用图像分割APP通过分割两个正交的二维切片来创建种子掩膜(seed mask)，一个在XY平面，另一个在XZ平面。然后，该示例将这两个分段插入到三维掩膜(3-D mask)中。然后将该遮罩传递给活动轮廓函数，以创建胸腔中肺部的三维分割。(此示例使用活动轮廓方法，但您也可以使用其他分割技术来实现相同的目标，例如flood-fill。

1、在XY和XZ维度提取中心切片

XY = V(:,:,160);
XZ = squeeze(V(256,:,:));   %squeeze函数三维转化为二维切片，即去除一个维度

2、使用图像显示功能查看二维切片

查看XY切片

figure
imshow(XY,[],'Border','tight');

查看XZ切片

imshow(XZ,[],'Border','tight');

3、可以在matlab中 APP——> 图像处理与计算机视觉——>Image Segmenter 中执行分割。或使用图像分割器命令，指定一个二维切片作为参数进行分割。如：imageSegmenter(XY)

imageSegmenter(XY)   %对XY切片进行分割

4、第3小节实现也可以通过命令行实现

BW = XY > 5.098000e-01;   %提取XY切片中阈值大于0.5098的像素，并保存为二值图像BW

5、在这个初始肺部分割之后，使用 APP——> 图像处理与计算机视觉——>Image Segmenter——>细化掩膜菜单中的选项来对掩膜BW进行进一步优化

加载掩膜：将粗分割的肺部掩膜加载到Image Segmenter APP中

反转掩膜(Invert Mask)：反转掩膜图像，以便肺部位于前景中

清理边框(Clear Borders)：移除除肺部之外的其他分割杂质

填充孔(Fill Holes)：填充肺部分割内的孔洞

形态学(Morphology)：使用形态学选项来平滑肺部分割的边缘

选择“形态学”——>选择腐蚀“erode”——>应用——>关闭形态学——>显示二元——>导出

便可将优化后的掩膜图像保存到工作区

形态学处理后二元图（即二值化图）：

第5小节也可以通过源码来实现：

BW = imcomplement(BW);    %反转掩膜（使肺部位于前景中）
BW = imclearborder(BW);   %清理边框（移除除肺部之外的其他分割杂质）
BW = imfill(BW, 'holes'); %填充孔（填充肺部分割内的孔洞）
%形态学处理
radius = 3;     %设置结构元素半径为3        
decomposition = 0;
se = strel('disk',radius,decomposition);  %创建盘形结构元素
BW = imerode(BW, se);   %形态学的腐蚀操作，得到的结果图像覆盖为BW
maskedImageXY = XY;     %XY切片赋值给maskedImageXY（即复制一个图像maskedImageXY）
maskedImageXY(~BW) = 0; %将maskedImageXY中BW图像中的非1区域都清0（即只保留maskedImageXY图像中的肺部区域）
imshow(maskedImageXY)   %显示图像maskedImageXY

6、XZ切片的掩膜优化处理

在XZ切片上执行上述相同的操作。对于XZ切片，全局阈值可创建了一个适当的分段(调用imbinarize)。与XY切片一样，在侵蚀操作中，指定半径为13以移除无关的小对象

BW = imbinarize(XZ);      %对XZ切片进行全局阈值提取
BW = imcomplement(BW);    %反转掩膜（使肺部位于前景中）
BW = imclearborder(BW);   %清理边框（移除除肺部之外的其他分割杂质）
BW = imfill(BW,'holes');  %填充孔（填充肺部分割内的孔洞）
%形态学处理
radius = 13;              %设置结构元素半径为13
decomposition = 0;
se = strel('disk',radius,decomposition);  %创建盘形结构元素
BW = imerode(BW, se);     %形态学的腐蚀操作，得到的结果图像覆盖为BW
maskedImageXZ = XZ;       %XZ切片赋值给maskedImageXZ（即复制一个图像maskedImageXZ）
maskedImageXZ(~BW) = 0;   %将maskedImageXZ中BW图像中的非1区域都清0（即只保留maskedImageXZ图像中的肺部区域）
imshow(maskedImageXZ)     %显示图像maskedImageXZ

三、创建种子MASK掩膜并使用活动轮廓（snakes）分割肺部

创建三维种子掩膜mask，并使用活动轮廓功能来分割肺部。创建与输入体积大小相同的逻辑三维体积（volume），并在适当的空间位置插入掩膜_XY和掩膜_XZ

%创建3-D种子mask
mask = false(size(V));    %创建与V矩阵维数大小相同的，元素都为0矩阵mask作为3-D种子mask
mask(:,:,160) = maskedImageXY; %maskedImageXY中的值赋值给mask的第三维度Z的第160切片
%reshape将maskedImageXZ重构为1*512*318矩阵，再与mask的第一维度切片进行或运算，目的同XY切片
mask(256,:,:) = mask(256,:,:)|reshape(maskedImageXZ,[1,512,318]);

使用此三维种子掩膜mask，通过活动轮廓方法（snakes）在三维体积中分割肺部。此操作可能需要几分钟时间。要获得高质量的分割，可使用直方图均衡化（histeq）在可用范围内扩展体素值

%使用活动轮廓方法（snakes）在三维体积中分割肺部，得到三维二值分割结果BW，三维原数据分割结果segmentedImage
V = histeq(V);             %对原始数据V进行直方图均衡化，增强对比度
BW = activecontour(V,mask,100,'Chan-Vese'); %调用活动轮廓方法，设置迭代100次，得到分割二值三维图像BW
segmentedImage = V.*single(BW);             %原数据V与分割二值数据BW进行点乘运算，将分割数据部分以原数据像素显示

之后可以通过运行命令volumeViewer(segmentedImage)可在Volume Viewer APP中查看三维分段肺效果。通过操纵渲染编辑器中的alphamap设置，可以获得肺部的良好视图

四、计算分割肺的体积

1、使用regionprops3函数来计算肺的体积个数

volLungsPixels = regionprops3(logical(BW),'volume'); %计算分割肺的体积

2、根据从原始文件元数据中收集的x、y和z维度中体素的间距，计算分割肺的体积大小（单位： $mm^{3}$ ）

spacingx = 0.76;     %单个体素（volume）X轴大小（单位mm）
spacingy = 0.76;     %单个体素（volume）Y轴大小（单位mm）
spacingz = 1.26*1e-6;%单个体素（volume）Z轴大小（单位mm）
unitvol = spacingx*spacingy*spacingz;   %单个体素（volume）的体积大小（单位:mm³）

volLungs1 = volLungsPixels.Volume(1)*unitvol;  %左肺总体积大小（单位:mm³）
volLungs2 = volLungsPixels.Volume(2)*unitvol;  %右肺总体积大小（单位:mm³）
volLungsLiters = volLungs1 + volLungs2         %分割肺部总体积大小（单位:mm³）

单个体素（volume）体积大小unitvol为7.2778*10^-7 $mm^{3}$

分割肺部所有体素（volume）总体积大小volLungsLiters为5.7726 $mm^{3}$

注：由于官方例程并未给定单个体素的大小单位，所有我觉得单个体素Z轴大小应该是有问题的

例程完整源码：

close all,clear all; %清空数据
%%%准备数据%%%%%%%%%%%%%%%%%%%
load chestVolume  %加载样本数据
V = im2single(V); %将CT扫描数据V从int16转换为single，以将其值标准化到范围[0，1]
volumeViewer(V);  %使用Volume Viewer APP查看胸部扫描图像

%%%肺部分割%%%%%%%%%%%%%%%%%%%
%在XY和XZ维度提取中心切片
XY = V(:,:,160);            
XZ = squeeze(V(256,:,:));   %squeeze函数三维转化为二维切片，即去除一个维度
figure(1),
imshow(XY,[],'Border','tight'); %使用图像显示功能查看二维切片XY
figure(2),
imshow(XZ,[],'Border','tight'); %使用图像显示功能查看二维切片XZ

%%%创建掩膜Mask_XY
%imageSegmenter(XY);   %对XY切片进行分割
%阈值提取
BW = XY > 5.098000e-01;   %提取XY切片中阈值大于0.5098的像素，并保存为二值图像BW
%创建XY切片的掩膜图像maskedImageXY
BW = imcomplement(BW);    %反转掩膜（使肺部位于前景中）
BW = imclearborder(BW);   %清理边框（移除除肺部之外的其他分割杂质）
BW = imfill(BW, 'holes'); %填充孔（填充肺部分割内的孔洞）
%形态学处理
radius = 3;     %设置结构元素半径为3        
decomposition = 0;
se = strel('disk',radius,decomposition);  %创建盘形结构元素
BW = imerode(BW, se);   %形态学的腐蚀操作，得到的结果图像覆盖为BW
maskedImageXY = XY;     %XY切片赋值给maskedImageXY（即复制一个图像maskedImageXY）
maskedImageXY(~BW) = 0; %将maskedImageXY中BW图像中的非1区域都清0（即只保留maskedImageXY图像中的肺部区域）
imshow(maskedImageXY)   %显示图像maskedImageXY

%创建XZ切片的掩膜图像maskedImageXZ
BW = imbinarize(XZ);      %对XZ切片进行全局阈值提取
BW = imcomplement(BW);    %反转掩膜（使肺部位于前景中）
BW = imclearborder(BW);   %清理边框（移除除肺部之外的其他分割杂质）
BW = imfill(BW,'holes');  %填充孔（填充肺部分割内的孔洞）
%形态学处理
radius = 13;              %设置结构元素半径为13
decomposition = 0;
se = strel('disk',radius,decomposition);  %创建盘形结构元素
BW = imerode(BW, se);     %形态学的腐蚀操作，得到的结果图像覆盖为BW
maskedImageXZ = XZ;       %XZ切片赋值给maskedImageXZ（即复制一个图像maskedImageXZ）
maskedImageXZ(~BW) = 0;   %将maskedImageXZ中BW图像中的非1区域都清0（即只保留maskedImageXZ图像中的肺部区域）
imshow(maskedImageXZ)     %显示图像maskedImageXZ

%创建3-D种子mask
mask = false(size(V));    %创建与V矩阵维数大小相同的，元素都为0矩阵mask作为3-D种子mask
mask(:,:,160) = maskedImageXY; %maskedImageXY中的值赋值给mask的第三维度Z的第160切片
%reshape将maskedImageXZ重构为1*512*318矩阵，再与mask的第一维度切片进行或运算，目的同XY切片
mask(256,:,:) = mask(256,:,:)|reshape(maskedImageXZ,[1,512,318]);

%使用活动轮廓方法（snakes）在三维体积中分割肺部，得到三维二值分割结果BW，三维原数据分割结果segmentedImage
V = histeq(V);             %对原始数据V进行直方图均衡化，增强对比度
BW = activecontour(V,mask,100,'Chan-Vese'); %调用活动轮廓方法，设置迭代100次，得到分割二值三维图像BW
segmentedImage = V.*single(BW);             %原数据V与分割二值数据BW进行点乘运算，将分割数据部分以原数据像素显示

volumeViewer(segmentedImage)        %在Volume Viewer APP中查看分割效果

volLungsPixels = regionprops3(logical(BW),'volume'); %计算整个分割肺的体积volume个数，主要由两部分肺（左右肺）组成

spacingx = 0.76;     %单个体素（volume）X轴大小（单位mm）
spacingy = 0.76;     %单个体素（volume）Y轴大小（单位mm）
spacingz = 1.26*1e-6;%单个体素（volume）Z轴大小（单位mm）
unitvol = spacingx*spacingy*spacingz;   %单个体素（volume）的体积大小（单位:mm³）

volLungs1 = volLungsPixels.Volume(1)*unitvol;  %左肺总体积大小（单位:mm³）
volLungs2 = volLungsPixels.Volume(2)*unitvol;  %右肺总体积大小（单位:mm³）
volLungsLiters = volLungs1 + volLungs2         %分割肺部总体积大小（单位:mm³）