实战：使用yolov3完成肺结节检测（Luna16数据集）及肺实质分割
实战：使用yolov3完成肺结节检测（Luna16数据集）
yolov3是一种被广泛应用于目标检测领域的深度学习模型，在此被应用于肺结节检测任务中。基于其三维特性，我们通过切片获取二维图像，并将其转换为yolov3能够处理的形式进行分析。
1. yolov3代码及原理
- 代码
- 原理
2. Luna16数据集 - 数据集介绍
- 转换成voc数据集格式
3. 检测效果预览
更新： 2019.9.10
对数据集进行肺实质分割。
更新： 2019.9.14
使用分割后的数据进行训练
更新： 2019.9.23
修改anchor size后重新训练
更新： 2019.12.5
发现一个注释的错误：
a = image_array.transpose(1,2,0)[:,:,20] #transpose是将(z,x,y)的三维矩阵转为(x,y,z)的矩阵
代码解读应该是transpose是将(z,y,x)的三维矩阵转为(y,x,z)的矩阵
更新 ：2019.12.11
发现应该坐标换算错误：
#如果图像进行旋转了，那么坐标也要旋转
if isflip :
x_max = 512 - x_min#修改前是 x_min = 512 - x_min 以下类似
y_max = 512 - y_min
x_min = 512 - x_max
y_min = 512 - y_max
x = 512 - x
y = 512 - y
代码解读补充 2020.11.19
博主临近毕业，可能有消息无法及时回复，见谅，有问题请留言。
很多ER向我求取数据集与代码。
其中代码的核心部分已在博客中详细讲解，并已上传至平台（其中存在一些小问题）。
如需提升自身能力，则建议先深入理解后再动手实践。
如果是想要快速上手，则可直接测试相关数据集。
目前我已经将VOC数据集分享至百度云存储服务，
并提供了相应的访问权限及提取密码zc0p。
赠人玫瑰，手留余香，如果我的博文/数据集/代码对你有所帮助，麻烦点个赞 ~~
补充说明：在解释数据集大小不为512×512的原因时，请注意以下几点。首先，在标准化处理过程中采用了分辨率分别为...的具体方法进行归一化处理；其次，在此过程中我特意选择了...作为基准单位来进行计算...]
1. yolov3代码及原理
1.1代码
采用C语言开发：
https://github.com/AlexeyAB/darknet
采用PyTorch框架构建：
https://github.com/eriklindernoren/PyTorch-YOLOv3
1.2原理
可查看：
PyTorch | YOLOv3基础与实现细节第一篇
[ "可查看"]
PyTorch | YOLOv3基础与实现细节第二篇
[ "可查看"]
PyTorch | YOLOv3基础与实现细节第三篇
[ "可查看"]
PyTorch | YOLOv3基础与实现细节第四篇
通过本篇文章，请深入理解YOLOv3算法基础及其对应的PyTorch实现细节。如需进一步探讨，请访问NotFound1911在博客上的相关文章
2.Luna16数据集
2.1 Luna16数据集介绍
LUNA16数据集概述_pursuit_zhangyu
补充说明
2.1.1.图像的翻转
.mhd文件中有CT图像的属性，其中有一条是：
TransformMatrix = 1 0 0 0 1 0 0 0 1
代码解读- TransformMatric:图像矩阵是否翻转的标志。(在现实中CT扫描的时候有的人是正卧，有的人是仰卧的，所以会导致图像会出现翻转的情况。)
- TransformMatrix = 1 0 0 0 1 0 0 0 1是正卧的，这里举一个反卧的例子：TransformMatrix = -1 0 0 0 -1 0 0 0 1 。当图像反卧的时候我们把图像的x，y坐标进行倒序调整，使其正卧。如下面代码所示（加载CT图像）：
def load_itk_image(filename):
with open(filename) as f:
contents = f.readlines()
line = [k for k in contents if k.startswith('TransformMatrix')][0]
offset = [k for k in contents if k.startswith('Offset')][0]
EleSpacing = [k for k in contents if k.startswith('ElementSpacing')][0]
# 把值进行分割
offArr = np.array(offset.split(' = ')[1].split(' ')).astype('float')
eleArr = np.array(EleSpacing.split(' = ')[1].split(' ')).astype('float')
CT = CTImage(offArr[0],offArr[1],offArr[2],eleArr[0],eleArr[1],eleArr[2])
transform = np.array(line.split(' = ')[1].split(' ')).astype('float')
transform = np.round(transform)#round() 方法返回浮点数x的四舍五入值
if np.any(transform != np.array([1, 0, 0, 0, 1, 0, 0, 0, 1])):#判断是否相等
isflip = True
else:
isflip = False
itkimage = sitk.ReadImage(filename)
numpyimage = sitk.GetArrayFromImage(itkimage)
if(isflip == True):
numpyimage = numpyimage[:,::-1,::-1] #::-1 倒序
return (numpyimage, CT, isflip)
代码解读此注释是对上述操作的时间记录：首次修改时间为20世纪4年代末至5年代初；后续修改时间则为6世纪初
此注解为我对上述操作的时间记录：首次修改时间为公元4年代末至公元5年代初；后续修改时间则为公元6世纪初
#如果图像进行旋转了，那么坐标也要旋转
if isflip :
x_max1 = 512 - x_min
y_max1 = 512 - y_min
x_min1 = 512 - x_max
y_min1 = 512 - y_max
x_max = x_max1
y_max = y_max1
x_min = x_min1
y_min = y_min1
x = 512 - x
y = 512 - y
代码解读2.1.2图像的预处理
如前所述，CT图像中的体素值范围为[-1000, +400]，因此我们进行了预处理步骤，并对该数值范围进行了裁剪处理以避免超出标准范围的情况。\n需要注意的是，在直接进行裁剪操作以保存单张图片时会出现问题——无法正确显示——因此建议在执行此操作之前将图像转换为RGB格式以便正确保存
def truncate_hu(image_array):
image_array[image_array > 400] = 400
image_array[image_array <-1000] = -1000
def normalazation(image_array):
max = image_array.max()
min = image_array.min()
#归一化
image_array = (image_array - min)/(max - min)*255
#image_array = image_array.astype(int)#整型
image_array = np.round(image_array)
return image_array
代码解读 #这个图的x，y坐标是反的，在生成标注的时候需要注意。
case_path = 'data/subset3/1.3.6.1.4.1.14519.5.2.1.6279.6001.123697637451437522065941162930.mhd'
numpyimage, _, _ = load_itk_image(case_path)
truncate_hu(numpyimage)
image_array = normalazation(numpyimage)
#3维图像切割
a = image_array.transpose(1,2,0)[:,:,20] #transpose是将(z,y,x)的三维矩阵转为(y,x,z)的矩阵
#转换成RGB类型并保存
im = Image.fromarray(a)
im = im.convert("RGB")
im.save("error/test1.png")
代码解读2.2 Luna16数据处理成VOC格式
2.2.1 遍历Luna16的标注文件
所有Luna16的标注都包含在annotations.csv文件内。该注释的内容本博已有详细说明。请注意以下几点：首先确保数据格式一致性；其次验证转换逻辑无误即可。
1.同一个CT图可能具有多个肺结节。
所以，在处理文件时，必须确保每个CT图像中的所有肺部异常体征都被彻底检查以避免遗漏任何可能影响后续处理的情况。
#查询文件 返回路径 -1表示无
def search(path=".", name="",fileDir = []):
for item in os.listdir(path):
item_path = os.path.join(path, item)
if os.path.isdir(item_path):
search(item_path, name)
elif os.path.isfile(item_path):
if name in item:
fileDir.append(item_path)
#print("fileDir:",fileDir)
return fileDir
代码解读2.坐标换算
上述博客中也提到的了，要完成世界坐标到图像坐标的换算。
#世界坐标转换到图像中的坐标
def worldToVoxelCoord(worldCoord, offset, EleSpacing):
stretchedVoxelCoord = np.absolute(worldCoord - offset)
voxelCoord = stretchedVoxelCoord / EleSpacing
return voxelCoord
代码解读肺结节的三个维度x、y、z之间的换算相对容易理解；此外还包括了w和h这两个维度的计算方法。w维度沿着x轴延伸表示，在y轴方向上则为h维度。需要注意的是,w和h均为长度单位,因此无需考虑坐标原点处的偏移问题。
w = np.round(r/CT.x_ElementSpacing).astype(int)
h = np.round(r/CT.y_ElementSpacing).astype(int)
代码解读这里的CT是我建的一个类，目的是方便CT各项数据的传递。
#定义CT图像类保存 中心点
class CTImage(object):
"""docstring for Hotel"""
def __init__(self, x_offset, y_offset, z_offset, x_ElementSpacing, y_ElementSpacing, z_ElementSpacing):
self.x_offset = x_offset
self.y_offset = y_offset
self.z_offset = z_offset
self.x_ElementSpacing = x_ElementSpacing
self.y_ElementSpacing = y_ElementSpacing
self.z_ElementSpacing = z_ElementSpacing
代码解读2.2.2 生成.xml文件
VOC数据集格式介绍：
Pascal Voc数据集深入解析_持久决心的博客-博客_pascal voc
VOC包含多种属性，在进行肺结节检测的过程中。其中某些属性无需考虑。我们关注的重点仅限于物体名称和坐标信息。具体实现方式如下所示：一个对应的XML标签示例。
<annotation>
<folder>VOC</folder>
<filename>0003.png</filename>
<size>
<width>512</width>
<height>512</height>
<depth>3</depth>
</size>
<object>
<name>nodule</name>
<bndbox>
<xmin>381</xmin>
<ymin>266</ymin>
<xmax>389</xmax>
<ymax>389</ymax>
</bndbox>
</object>
</annotation>
代码解读生成.xml的函数为（Luna16的CT图切割之后都是3512512的）：
def beatau(e,level=0):
if len(e)>0:
e.text='\n'+'\t'*(level+1)
for child in e:
beatau(child,level+1)
child.tail=child.tail[:-1]
e.tail='\n' + '\t'*level
def ToXml(name, x, y, w, h):
root = Element('annotation')#根节点
erow1 = Element('folder')#节点1
erow1.text= "VOC"
erow2 = Element('filename')#节点2
erow2.text= str(name)
erow3 = Element('size')#节点3
erow31 = Element('width')
erow31.text = "512"
erow32 = Element('height')
erow32.text = "512"
erow33 = Element('depth')
erow33.text = "3"
erow3.append(erow31)
erow3.append(erow32)
erow3.append(erow33)
erow4 = Element('object')
erow41 = Element('name')
erow41.text = 'nodule'
erow42 = Element('bndbox')
erow4.append(erow41)
erow4.append(erow42)
erow421 = Element('xmin')
erow421.text = str(x - np.round(w/2).astype(int))
erow422 = Element('ymin')
erow422.text = str(y - np.round(h/2).astype(int))
erow423 = Element('xmax')
erow423.text = str(x + np.round(w/2).astype(int))
erow424 = Element('ymax')
erow424.text = str(y + np.round(h/2).astype(int))
erow42.append(erow421)
erow42.append(erow422)
erow42.append(erow423)
erow42.append(erow424)
root.append(erow1)
root.append(erow2)
root.append(erow3)
root.append(erow4)
beatau(root)
return ElementTree(root)
代码解读2.2.3 生成训练集、测试集、验证集
Luna16总计有1186个肺结节。主要流程开始时，将其中70%划分为训练集数据集，并从剩下的中选择约25%作为测试集数据集和验证数据集。
import random
import numpy as np
getTrain = []
getTest = []
getVal = []
arrTrain = []
arrTest = []
arrVal = []
arr = np.arange(1,1186)
#训练集
num = 0
for i in np.random.permutation(arr):
if num >829:
break
arrTrain.append(i)
num = num + 1
getTrain=np.array(arrTrain)
getTrain.sort()
f=open('train.txt','w')
for i in getTrain:
f.write('{:04d}\n'.format(i))
f.close()
#测试集
num = 0
for i in np.random.permutation(arr):
if num >178:
break
arrTest.append(i)
num = num + 1
getTest=np.array(arrTest)
getTest.sort()
f=open('test.txt','w')
for i in getTest:
f.write('{:04d}\n'.format(i))
f.close()
#验证集
num = 0
for i in np.random.permutation(arr):
if num >178:
break
arrVal.append(i)
num = num + 1
getVal=np.array(arrVal)
getVal.sort()
f=open('val.txt','w')
for i in getVal:
f.write('{:04d}\n'.format(i))
f.close()
代码解读生成之后.txt就有对应的图片ID。
下一步将采用yolov3自带的voc_label.py来进行解析，并生成相应的label。在处理过程中，默认将下列代码中的参数difficult设置为默认值0。
import xml.etree.ElementTree as ET
import pickle
import os
from os import listdir, getcwd
from os.path import join
sets=[('', 'train'), ('', 'val'), ('', 'test')]
classes = ["nodule"]
def convert(size, box):
dw = 1./(size[0])
dh = 1./(size[1])
x = (box[0] + box[1])/2.0 - 1
y = (box[2] + box[3])/2.0 - 1
w = box[1] - box[0]
h = box[3] - box[2]
x = x*dw
w = w*dw
y = y*dh
h = h*dh
return (x,y,w,h)
def convert_annotation(year, image_id):
in_file = open('E:/LUNA/VOCdevkit/VOC%s/Annotations/%s.xml'%(year, image_id))
out_file = open('E:/LUNA/VOCdevkit/VOC%s/labels/%s.txt'%(year, image_id), 'w')
tree=ET.parse(in_file)
root = tree.getroot()
size = root.find('size')
w = int(size.find('width').text)
h = int(size.find('height').text)
for obj in root.iter('object'):
difficult = 0
cls = obj.find('name').text
if cls not in classes or int(difficult)==1:
continue
cls_id = classes.index(cls)
xmlbox = obj.find('bndbox')
b = (float(xmlbox.find('xmin').text), float(xmlbox.find('xmax').text), float(xmlbox.find('ymin').text), float(xmlbox.find('ymax').text))
bb = convert((w,h), b)
out_file.write(str(cls_id) + " " + " ".join([str(a) for a in bb]) + '\n')
wd = getcwd()
for year, image_set in sets:
if not os.path.exists('VOCdevkit/VOC%s/labels/'%(year)):
os.makedirs('VOCdevkit/VOC%s/labels/'%(year))
image_ids = open('VOCdevkit/VOC%s/ImageSets/Main/%s.txt'%(year, image_set)).read().strip().split()
list_file = open('%s_%s.txt'%(year, image_set), 'w')
for image_id in image_ids:
list_file.write('%s/VOCdevkit/VOC%s/JPEGImages/%s.png\n'%(wd, year, image_id))
convert_annotation(year, image_id)
list_file.close()
os.system("cat 2007_train.txt 2007_val.txt 2012_train.txt 2012_val.txt > train.txt")
os.system("cat 2007_train.txt 2007_val.txt 2007_test.txt 2012_train.txt 2012_val.txt > train.all.txt")
代码解读这一步完成之后可以看到根目录下有以下文件：
打开后，里面包含的是路径：
对应训练文件的Labels也出现了：
2.2.4 修改anchor size
anchor size needs to be adjusted according to the specific dataset, but I haven't made any changes here. If someone wants to modify it, they can refer to this link.
anchor size needs to be adjusted according to the specific dataset, but I haven't made any changes here. If someone wants to modify it, they can refer to this link.
YOLO-v3模型参数设置中的锚框配置_m_buddy的博客-博客中的锚框配置
该文详细探讨了基于K-means算法实现目标检测中的锚框生成过程
3. 检测效果预览
3.1 训练
在执行检测任务前，必须对生成的VOC格式中的Luna16数据集进行训练。请确保确认第二部分所生成的所有相关文件（如.xml、.jpgs和.txt）是否已经完成。随后，则是需要修改yolov3网络结构的工作。
每个yolo层的class参数设置为1（仅用于肺结节检测），而filters参数则设置为3 \times (1 + 5)。同时对.data文件和.names文件进行相应的调整。
待完成修改后的.cfg、.data和.nam文件配置即可开始训练；在当前阶段（即即将更新）进行中的是基于C语言实现的yolov3模型的学习过程
3.2 训练结果显示
~~还在训练中，训练完了更新。 ~~
训练了大概3天，迭代了48000次，loss值约0.03（GTX1080）。
原图：
预测：
标注图（根据标注生成）：
在应用yolov3算法进行肺结节检测时，发现其在实际应用中仍存在一定局限性。具体表现为：1）迭代训练不足导致模型尚未达到充分收敛的状态；2）有必要对训练数据集进行聚类分析以优化参数设置，并调整anchor框大小与比例参数设置以进一步提升检测精度。
3.3 存在的问题
1. 数据不平衡问题。
在本次训练期间，仅涉及了标准病例（具有肺结节）作为训练案例，并未包含属于其他类别但具有相似特征的情况
例如这张图片是根据标注生成的。观察到肺结节确实很小。然而实际上在图像左下方还有许多类似肺结节的节点但这些节点并不是真正的肺结层归因于我仅使用了正例即真实存在的肺结节数据进行标注
当class标记为1时，则表明存在肺结节。然而，在仅关注类似肺结节的区域分析时，则需要分阶段处理：第一阶段是识别潜在的可疑区域；第二阶段是从这些可疑区域中筛选出真实的病变区域。对于正负样本比例失衡的问题而言，在模型训练过程中适当采用focal loss损失函数是一种可行的选择。
2.数据结构信息不完整
这一现象主要源于CT图像从三维转换为二维的过程，并导致结构信息的丢失。Luna16数据集共计存储了约112GB的数据，在经二维处理后体积缩减至约247MB左右。所有图片均由标注生成（labelImage）。通过观察第一张图片发现其仅在右下角显示肺结节特征；而观察第二张图片则显示其左上角存在肺结节特征。然而这两者的表现却极为相似的原因在于：由于CT图像具有三维特性，在切割过程中采用Z轴方向进行操作；这也说明标注的肺结节在Z轴方向上的位置存在差异性；然而令人意外的是，在后续分析中发现这两组切片在Z轴坐标上的数值却极为接近；因此造成上述现象的原因尚不明确需进一步探究。我认为解决这一问题的有效方法应是对每个病灶体赋予一个概率值参数；当计算出的概率值超过设定阈值时即可实现清晰显示。
以上，简单的把yolov3用于实际数据集（Luna16）处理。
把Luna转换成VOC格式的脚本已上传：
<>
里面包含了两个脚本，使用到的所有函数已在第二部分讲解。
这两个脚本不是最新的。
将L Luna 1 6 数据 集 转 换 为 VOC 格 式 并 使 用 于 进 行 肺 实 质 的 深 度 学 习 分 割 和 产 生 ， Mat 文件 在 C S D N 平 台 上 的 深 度 学 习 代 码 库 中 下 载
更新： 2019.9.10
对数据集进行肺实质分割。
1.标准化数据，查看数值分布：
2.计算位于肺周围区域的平均像素值，并将被遮挡或损坏的图像进行重新归一化处理。通过K-means算法分别识别出目标物体（不透明组织）与背景区域（透明组织），其中背景区域重点关注的是肺部区域。这种处理方法仅针对图像中心区域实施，并尽力避免对未参与区域产生影响
- 采用腐蚀与膨胀的操作进行处理。如需了解详细内容可参考https://www.jianshu.com/p/05ef50ac89ac。腐蚀操作能够有效去除某些区域的颗粒状结构；随后采用较大的膨胀作用于肺部血管及其周围组织，并以消除血管周围的黑色噪声为主，在此过程中尤其是由于不透明射线导致的肺部阴影区域也会得到适当的清理与处理。
- 获取肺部掩膜。选择CT图像中心区域的目的在于去除周围非感兴趣区域（如床板等），随后对图像进行连通域分析并确定较大的连通域作为肺部掩膜基础。在此基础上实施一次较大范围的膨胀操作用于填充和提取肺部掩膜区域。**此时边界可能出现不规则现象：可以通过腐蚀操作优化边界细节；然而需要注意的是，在此步骤之前已经进行了充分的大范围膨胀操作；因此可能会导致内部空洞扩大风险（因为后续的腐蚀处理可能进一步放大原有空洞）。
参考：https://github.com/booz-allen-hamilton/DSB3Tutorial
更新：2019.9.14
使用分割后的数据进行训练
经过大约四天的训练后确实比之前有了显著提升
下图左边是预测，右边是根据标注生成的ground truth。
更新： 2019.9.23
修改anchor size后重新训练
对数据集进行kmean 是聚类，修改anchor size，这里需要主要几点：
在YOLACT框架中采用了9个k-means聚类算法来处理图像分割任务。然而，在针对肺癌结节检测的数据集中观察到这一设置存在不合理之处：该算法几乎不可能达到数量为9的状态。当采用该数量时会导致算法无法收敛或计算过程陷入死循环。
数据集中的bounding boxes所占的比例较小，在经过归一化处理后（即将其缩放到0~1区间），这些数值普遍较低（只有约0.01x）。这可能会导致计算负担加重并可能影响收敛性。在这里，我采用的方法是先将这些bounding boxes的尺寸放大十倍，在完成聚类操作之后再将其缩小回原比例。
基于标注生成的标注图像，在实际应用中无法获取到发现肺结节的bounding boxes大致呈现比例为1:1。然而，在yolov3算法中由于采用了3 个特征尺度来进行位置预测，在此处我们选择聚类中心的数量定为**3。
通过反复进行多次聚类运算，可以获得较为合适的anchor大小（但有时可能出现不收敛的情况，并导致计算出现问题）：
这个之84.07%已经很高了，比voc 67%左右高出不少。
当然也测试了一下其他聚类中心（1个或者2个）：
但都没有使用3个聚类中心的好。
训练：
检测效果：
