七、联邦学习和医疗保健
随着越来越多的计算机和硬件技术变得如此容易获得,从临床机构到保险公司,从患者到制药行业,不同的医疗保健利益相关方提供了大量的分析数据。这些海量数据是一座金矿,有助于发现有助于设计人工智能集成医疗系统的见解,旨在以合理的成本提供更好的结果和质量。
无论医疗保健数据产生的数量有多少,它们仍然是碎片化的,法律、伦理和隐私方面的考虑阻碍了用于稳健研究的大规模数据分析。例如,正如您在第 3 和 4 章中所看到的,贝斯伊斯雷尔女执事医疗中心收集的 EHR 数据虽然仍然是一个大型数据集,但缺少白人和非白人人口分布、年龄分布差异等方面的信息。而包含更多这种无代表群体的数据可能存在于其他地方。因此,需要一种更巧妙的思维方式。
联邦学习通过将模型引入数据,而不是相反,帮助我们解决隐私和法律限制等问题。在这一章中,你将深入研究联邦机器学习。什么是 TensorFlow Federated?有哪些不同的隐私机制?本章的目的不是向您介绍一个新颖的案例研究,而是更多地了解 TensorFlow 联合生态系统(联合、隐私和加密)及其功能。
介绍
联邦学习(FL)是一个分布式机器学习概念,它允许对分散的数据进行模型训练,同时解决每个利益相关者的数据传输、隐私和安全问题。一个外语系统有四个主要组成部分,它们同步工作以进行联邦学习:
-
中央服务器/节点:协调本地模型的培训和部署,并作为创建全球模型的平台。本地模型是在本地节点/边缘设备上训练的模型,而全局模型是使用来自本地节点的权重更新权重的模型。
-
本地服务器/设备/节点:这是真实世界数据所在的地方。它们通常是收集客户数据的已安装机器的边缘设备。
-
本地模型:这是任何类型的机器学习模型,它对本地服务器中存在的数据进行训练。这些模型学习特定于本地设备的数据。
-
全局模型:由不同局部模型的信息组合而成的最终模型。
联邦学习是如何工作的?
联邦学习培训有四个关键步骤。
步骤 1:从中心节点转移初始模型(见图 7-1
图 7-1
第一步
-
从中央服务器获得的初始模型根据模型所有者可用的数据进行训练(即,该模型根据中央服务器上可用的数据进行训练)。
-
这个全局模型然后通过网络被传送到所有的本地节点。
第二步:模特训练
-
任何类型的机器学习模型,从朴素贝叶斯和 SVM 这样的基本模型到 DeepNets,都可以被训练。
-
一小部分客户被选择用于本地模型训练,因为选择大量的客户在性能和成本上具有递减的回报。
-
本地节点的计算资源用于训练,这节省了中央服务器的计算时间和资源。
-
有时,本地节点上的数据不足,这可能使该节点对全局模型的贡献无效,因此需要安全聚合等技术,这种技术允许使用公钥-私钥在节点之间共享数据。此外,这种技术有助于防止个人数据泄露问题。
步骤 3:本地模型转移到中心节点(见图 7-2
图 7-2
第三步
-
训练之后,所有模型都可以传回中央服务器。对于边缘设备,这可能导致巨大的网络开销(跨设备培训),而在跨孤岛的联邦培训(作为本地节点的组/机构)中,这种影响不太明显。
-
有时,模型可能会受到敌对攻击,这些攻击有助于识别用于训练模型的用户敏感数据。因此,为了防止这种攻击,可以使用实现差分隐私或安全聚合等技术的隐私保护层。注意,差分隐私原则上也可以应用于本地,而不是全局。在后面的章节中会有更多的介绍。
步骤 4:中心节点聚集来自所有本地模型的结果。
-
联合平均不仅仅是输出概率或多数表决的简单平均。
-
需要学习的任何参数,例如深度学习模型对权重更新起作用。因此,全局权重向量是通过对损失度量进行加权并用观察到的样本数进行归一化来决定的。通过这种方式,您可以获得更多的权重表示,这在统计上(样本数)会带来更好的性能。
-
根据结果如何从本地节点传输,可以有许多其他平均技术。
联邦学习的类型
根据数据在 FL 训练过程中如何分布在多个本地节点上,您可以将 FL 分为三个主要类别。
水平联邦学习
在水平联邦学习中,不同本地节点的数据集具有相同的特征空间集,但是样本的重叠量是最小的。
这是跨设备设置的自然划分,不同的节点/用户试图改进一项共同的任务,比如在移动应用程序上使用 GBoard 打字时的键盘建议,或者使用可穿戴设备数据进行疾病的风险预测。见图 7-3 。
图 7-3
横向联邦学习。来源:Kurupathi 等人的“面向隐私保护人工智能的联邦学习调查”
垂直联邦学习
这里,不同本地节点的数据集具有相同的样本/人集合,但是特征空间的重叠量可以根据组织数据而不同。
当多个组织进行协调时,他们可以期待实施垂直 FL。使用特征对齐方法来对齐不同个体的特征,然后训练单个模型。对齐是隐私保护的,这意味着不容易识别受保护的信息。这可以通过加密来实现。您将在安全聚合讨论中了解更多相关信息。
一些例子是保险公司和银行公司在共享客户的公共数据上的联合协作。标签可以是默认利率或任何欺诈交易。在医疗保健领域,不同的医院可以共享他们擅长的不同检查的信息,以绘制出患者的综合病史。见图 7-4 。
图 7-4
垂直联邦学习。来源:Kurupathi 等人的“面向隐私保护人工智能的联邦学习调查”
联邦迁移学习
这是在特征空间和样本不同的情况下实现的。比如说有一群医院想做乳腺癌研究。每家医院都有一组不同的患者(样本),他们可能会捕获不同的指标(特征空间),在两个维度上都有一些最小的捕获。
通常,使用有限的公共样本集在两个特征空间之间学习公共表示,然后将其应用于获得仅单侧特征样本的预测。见图 7-5 。
图 7-5
联邦迁移学习。来源:Kurupathi 等人的“面向隐私保护人工智能的联邦学习调查”
隐私机制
只有当除了介绍中讨论的流程之外,还实现了多种隐私机制时,FL 才在现实世界中被广泛接受。虽然数据集驻留在本地节点上,但是可以对模型参数进行重新设计,以获得有关数据的信息。此外,可以一起应用多种隐私机制技术,以确保个人身份/数据的更强大的安全性。
如图 7-6 所示,恶意攻击者可以尝试将梯度与本地梯度更新进行匹配,并重建数据。
图 7-6
梯度匹配攻击
有许多隐私机制技术,但在这一章中,我们将讨论两种在当前的 FL 系统中最常用的技术。
要更深入地了解隐私机制和衡量其有效性的方法,请参考瓦格纳等人在 2018 年发布的论文“技术隐私指标:系统调查”。
安全聚合
安全聚合是一种保护隐私的机器学习技术,当以安全的方式从各个用户设备更新时,它依赖于多方计算来计算模型参数的和。
2017 年,谷歌最初在题为“隐私保护机器学习的实用安全聚合”的论文中提出了一种安全聚合技术。关于数学细节,你可以看看论文,但现在让我们直观地理解它。
-
公钥和私钥是使用模式生成的。
-
公共密钥与每个本地节点共享。
-
这些密钥用于加密模型参数的更改。
-
所有的本地节点使用像加法或乘法这样的数学运算来累积模型权重。
-
累积的更改被发送到中央服务器,中央服务器使用私钥解密数据。
在上述过程中需要注意两件事
-
我们可以对加密数据本身进行 ALU(算术和逻辑)运算,因为加密本质上是同态的。我们可以在不解密的情况下对数据进行 ALU 运算。
-
中央服务器看到累积的结果,这些结果可以使用私钥解密。
此外,某些用户可能会由于网络问题而突然退出。只有当总和来自至少 n 个本地节点时,中央服务器上发生的任何变化才会发生。
在内部,TFF (TensorFlow Federated)使用 TensorFlow Encrypted 来执行这个练习,但为了简单起见,让我们使用 pallier 包来看看这是如何工作的。
您将使用 python-pallier,它使用 Paillier 加密系统(一种同态加密方案,参见 https://blog.openmined.org/the-paillier-cryptosystem/
了解同态加密如何工作)。
import phe
import numpy as np
# Generate Public and Private Key
public_key, private_key = phe.generate_paillier_keypair(n_length=1024)
weight1 = np.random.rand(10)
weight2 = np.random.rand(10)
# Note : This is a simple addition but it can be more complex as well
sum_of_local_weights = np.add(weight1, weight2)
print("Addition Of w1 and w2: " + str(sum_of_local_weights))
encrypted_w1 = [public_key.encrypt(i) for i in weight1]
encrypted_w2 = [public_key.encrypt(j) for j in weight2]
encrypted_sum_of_w1_and_w2 = [i+j for i,j in zip(encrypted_w1, encrypted_w2)]
decryped_sum_of_w1_and_w2 = [private_key.decrypt(k) for k in encrypted_sum_of_w1_and_w2]
print("Addition Of Encrypted Number: " + str(decryped_sum_of_w1_and_w2))
输出
Addition Of w1 and w2: [0.01965569 1.38181926 0.95724207 1.40539024 0.56162914 1.26444545
0.84660776 0.55585975 1.60470971 0.74662359]
Addition Of Encrypted Number: [0.01965569240712506, 1.381819260975988, 0.957242068080129, 1.4053902417875905, 0.5616291366817605, 1.2644454455590868, 0.8466077626079891, 0.5558597475342251, 1.604709707486859, 0.7466235859816883]
您可以看到,通过共享聚合结果本身,同态加密和多方(以确保健壮性和更好的正常运行时间)可以多么容易地保护个人数据。
正如您现在可能已经想到的,这种技术在计算方面包含大量开销,这些开销可以随着本地节点的数量和参数向量的大小而扩展。
差异隐私
差分隐私是一种隐私机制,它试图量化由于在本地节点(本地 DP)或在聚集级别(全局 DP)向数据添加噪声而导致的隐私量,使得最终分析保持不变。我们通过一个例子来了解一下。
假设你在班上做了一项调查,看看有多少学生对绿色色盲。你计划在解释一些概念时包括大量可能使用绿色的视觉效果。
目的:大多数人对绿色不是色盲吗?
假设您管理调查,结果如图 7-7 所示。
Note
为简单起见,我们使用非常小的样本量。
图 7-7
差异隐私
让我们假设你的二次研究告诉你,来自某个种族的人倾向于对绿色表现出色盲。因此,如果给你提供数据 D*,你可以以某种方式识别教室中的那些特定个人。但是,如果你添加一些噪音,使数字变得不直观,无法精确到教室的某个部分,那该怎么办呢?*
*这正是差别隐私所保证的。它保护了参与分析的个人,但不影响结果,就像上面发现班上大多数人对绿色不敏感的情况一样。
差分隐私引入了一个称为ε的度量,它量化了数据分布的接近程度:
)
如果 ε = 0,那么你就有了一个精确的分布,你就达到了峰值隐私。f(D n )表示数据函数,f(D’n)表示加入噪声后的数据函数。
在实践中,拉普拉斯分布和正态分布被用于生成对查询的回答,因为这些函数更有可能预测更接近平均值的数字(< =1 与平均值的标准偏差;对于标准常态是 68%,而对于拉普拉斯是 74% (b=1),但没有给出正确的答案。这里的平均值就是你的真实值。
如你所知,如果从任何分布中抽取足够多的随机样本,你可以用中心极限定理来估计它的平均值。同样,如果向包含来自本地节点的数据的数据库发出多个查询,则可以形成平均值的估计。
举个例子,
-
调查中非色盲学生的数量是多少?
-
属于“这个”族裔的学生人数是多少?
-
总人数是多少?
因此,每次您抛出一个数字,您就给了对手一个更好的机会来猜测来自本地节点或其任何特征的正确数字/数据。
在实现差分隐私时,您必须确保两个概率分布尽可能接近。在小样本情况下,噪声可以完全改变数据;对于大量样本,噪声的影响有限(因为引入了更多的变化,所以单个噪声无法很好地掩盖所有样本)。因此,设计噪声函数有时是一项极其困难的任务。
由于噪声有时会淹没小样本,您可以引入另一个名为 δ 的参数,这是一个帮助您丢弃罕见类别的阈值。因此,独特的差分隐私机制实际上是两个因素的函数:
-
阈值( δ
-
噪音量( ε
TensorFlow Privacy 是 TensorFlow 生态系统中的一个库,用于训练具有训练数据隐私的机器学习模型。该库提供了三个不同的特性:
-
训练算法,特别是梯度下降
-
它通过剪切梯度来限制单个数据点在结果梯度计算中的影响。
-
它通过向剪切的梯度添加随机噪声,使得梯度值与训练批次中的任何特定点无关。
-
-
隐私机制的选择和配置(超参数调整)以应用于收集的每个集合(模型梯度、批量标准化权重更新、度量)
-
绩效指标
-
隐私预算
-
埃普西隆
-
Note
差分隐私是一种独立的隐私维护技术,可用于 FL 架构,在这种情况下,更新来自多方。
TensorFlow 联邦
TensorFlow Federated (TFF)是一个通过模拟实验在本地应用联邦学习的开源框架。TFF 使开发人员能够在他们的模型和数据上模拟包含的联邦学习算法,以及试验新的算法。
TFF 的界面分为两层:
-
联邦学习(FL) API :这一层提供了一组高级接口,允许开发人员将联邦培训和评估的实现应用到他们现有的 TensorFlow 模型中。
-
联邦核心(FC) API :系统的核心是一组底层接口,通过结合 TensorFlow 和分布式通信操作符来表达新颖的联邦算法。
输入数据
您将使用疟疾数据集,该数据集包含总共 27,558 个细胞图像,其中包含来自分割细胞的薄血涂片图像的相同寄生和未感染细胞实例。数据集可以从 https://ceb.nlm.nih.gov/proj/malaria/cell_images.zip
.
中获得
TensorFlow 联合图书馆生态系统中存在模拟数据集,但疟疾数据集接近医疗保健领域。在下一章中,您将看到医学图像分析如何处理 2D 和 3D 图像数据,因此这是一个良好的开端。
疟疾数据集包含两个类,如图 7-8 所示:
图 7-8
寄生和未感染细胞的例子
-
被寄生的(即被感染的细胞)
-
非寄生细胞(又称未感染细胞)
首先从本地目录加载数据,并查看感染和未感染图像样本的分布。
import os
import glob
BASE_DIR = os.path.join('./Data')
parasitized_dir = os.path.join(BASE_DIR,'Parasitized')
uninfected_dir = os.path.join(BASE_DIR,'Uninfected')
parasitized_files = glob.glob(parasitized_dir+'/*.png')
uninfected_files = glob.glob(uninfected_dir+'/*.png')
len(parasitized_files), len(uninfected_files)
输出
(13779, 13779)
看起来你对这两个类有一个平衡的表示。
联邦学习需要一个联合数据集(来自多个用户的数据集合,也称为本地节点)。任何联邦数据都应该是非 iid 的,这意味着不同的客户机应该至少有一些合理的相似分布(特定于本地节点的特征会影响每个系统上的数据分布)。
在您的情况下,您不会有不同的数据集分布,但通过可视化来探索它们会很好。
自定义数据加载流水线
如果您使用的是 tff 库中已经存在的模拟数据集,您可以简单地调用函数load_data().
_train, _test = tff.simulation.datasets.<dbname>.load_data()
load_data()
返回的数据集是tff.simulation.ClientData
的实例,它枚举本地节点的集合来构造一个代表特定节点数据的tf.data.Dataset
,并查询各个数据元素的结构。
因为您没有使用预先模拟的数据集,所以您需要自己构建一个。因为您的目录结构是按照以下方式组织的
Data/
...Parasitized/
......image_1.png
......image_2.png
...Uninfected/
......image_1.png
......image_2.png
您可以利用 tf.keras 预处理函数image_dataset_from_directory
。
调用image_dataset_from_directory(data_directory, labels="inferred")
将返回一个tf.data.Dataset
,它从子目录Parasitized
和Uninfected
中产生一批图像,以及标签 0 和 1 (0 对应于Parasitized
,1 对应于Uninfected
)。
tf.keras.preprocessing.image_dataset_from_directory(
BASE_DIR, labels='inferred', label_mode='int',
class_names=None, color_mode='rgb', batch_size=32, image_size=(256,
256), shuffle=True, seed=None, validation_split=None, subset=None,
interpolation='bilinear', follow_links=False
)
在上面的函数中,还有一个调整图像大小的选项,但是要调整图像大小,您需要知道正确的调整后的形状。由于这些是细胞图像,它们可能有不同的形状。让我们快速检查一下,然后使用预处理函数加载数据。
由于您有大约 30k 个图像,顺序加载每个图像会花费一些时间,因此您应该尝试在不同的 CPU 内核上并行化操作,并使用 OpenCV 库返回每个图像的形状。
首先加载库,并使用内置的 os 库来计算 CPU 数量。
from joblib import Parallel, delayed
import os
nprocs = os.cpu_count()
为了不中断其他应用程序的计算资源,您使用的 CPU 比总数少一个。这只是一个可以遵循的好习惯。
您将使用 OpenCV 3 来读取图像。可以通过运行以下命令来下载它:
pip install opencv-python==3.4.6.27
def load_image_shape(img):
return cv2.imread(img).shape
results = Parallel(n_jobs=nprocs-1)(delayed(load_image_shape)(img_file) for img_file in parasitized_files + uninfected_files)
print('Min Dimensions:', np.min(results, axis=0))
print('Avg Dimensions:', np.mean(results, axis=0))
print('Median Dimensions:', np.median(results, axis=0))
print('Max Dimensions:', np.max(results, axis=0))
Min Dimensions: [40 46 3]
Avg Dimensions: [132.98345308 132.48715437 3\. ]
Median Dimensions: [130\. 130\. 3.]
Max Dimensions: [385 394 3]
注意,这个过程只是加快了加载速度,但是您仍然要在内存中加载完整的数据,对于大型数据集,通常不建议这样做。在这种情况下,可以使用生成器,它会在需要时加载数据。
因此,图像形状的中值尺寸为 130,因此您可以安全地将所有图像缩放为标准形状(128,128,3)。
此外,为了整形,Keras 预处理库将使用双线性插值,这是默认选项,因此您只需使用它( bi 在这里表示图像的二维(x,y))。
import numpy as np
import pandas as pd
import tensorflow as tf
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
import tensorflow_federated as tff
IMG_HEIGHT = 128
IMG_WIDTH = 128
BATCH_SIZE = 32
train_ds = tf.keras.preprocessing.image_dataset_from_directory(BASE_DIR,
seed=123,
labels='inferred',
label_mode='int',
image_size=(IMG_HEIGHT, IMG_WIDTH),
color_mode='rgb',
subset="training",
shuffle=True,
validation_split = 0.2,
batch_size= BATCH_SIZE)
Found 27558 files belonging to 2 classes.
Using 22047 files for training
.
val_ds = tf.keras.preprocessing.image_dataset_from_directory(BASE_DIR,
seed=123,
labels='inferred',
label_mode='int',
image_size=(IMG_HEIGHT, IMG_WIDTH),
color_mode='rgb',
subset="validation",
shuffle=True,
validation_split = 0.2,
batch_size= BATCH_SIZE)
Found 27558 files belonging to 2 classes.
Using 5511 files for validation.
Note
到目前为止,TF 2.2.0 是我们旅程中的版本,它不支持image_dataset_from_directory
函数,所以建议使用最新的 TensorFlow 联邦库,它默认安装 TF 2.3.0。在 TF 2.3 和更高版本中,支持image_dataset_from_directory
。
您还可以看到标签映射到的类名。
class_names = train_ds.class_names
print(class_names)
输出:
['Parasitized', 'Uninfected']
这意味着整数 0 是类Parasitized
,1 是类Uninfected
。
疟疾数据集是一个大型数据集。根据您的机器设置,您可以将全部数据加载到内存中,也可以不加载。为了避免运行时出现任何问题,您将把联邦数据创建放在一个try-catch
块中。
NUM_CLIENTS = 10 # Local Nodes
CLIENT_LR = 1e-2
SERVER_LR = 1e-2 # Central Node
NUM_BATCH_CLIENT = int(len(train_ds)/NUM_CLIENTS)
import collections
client_train_dataset = collections.OrderedDict()
skip = 0
try :
for i in range(1, NUM_CLIENTS+1):
client_name = "Client_" + str(i)
take = NUM_BATCH_CLIENT
client_data = train_ds.skip(skip).take(take)
x_train, y_train = zip(*client_data)
print(f"Adding data from Batch No {skip} to {take*i} for client : {client_name}")
# We are going to unbatch and load the data to prevent data dropping in creating client data later on
data = collections.OrderedDict((('label', [y for x in y_train for y in x]),
('pixels', [y for x in x_train for y in x])))
client_train_dataset[client_name] = client_data
skip = take*i
except Exception as e:
print("Memory Error - Client Data creation stopped")
print(f"Total number of clients created are {len(client_train_dataset)}")
NUM_CLIENTS = len(client_train_dataset)
输出
Adding data from Batch No 0 to 68 for client : Client_1
Adding data from Batch No 68 to 136 for client : Client_2
Adding data from Batch No 136 to 204 for client : Client_3
Adding data from Batch No 204 to 272 for client : Client_4
Adding data from Batch No 272 to 340 for client : Client_5
Memory Error - Client Data creation stopped
Total number of clients created are 4
在上面的代码中,您试图创建一个有序字典,以便在从张量切片创建客户数据时保持客户的顺序。
正如我所说的,您可以从预期数量的客户端开始,但是根据本地可用的计算资源,您也可以预期更少的客户端数量。在这里,您最终只剩下四个客户端,因此减少了用于训练的数据。现在你不应该担心这个问题,因为对于一个更大的机器来说,这样的问题很容易解决。
此外,根据 TFF 团队的说法,“我们近期的未来路线图包括一个高性能的运行时,用于非常大的数据集和大量客户端的实验。”
接下来,通过传递客户端数据的键值对,在模拟环境中创建客户端数据(参见图 7-9 )。
图 7-9
来自联合数据的图像
train_dataset = tff.simulation.FromTensorSlicesClientData(client_train_dataset)
sample_dataset = train_dataset.create_tf_dataset_for_client(train_dataset.client_ids[0])
sample_element = next(iter(sample_dataset))
本地节点的训练样本总数为
len(sample_dataset)
输出
2176
plt.imshow(sample_element['pixels'].numpy().astype('uint8'))
plt.grid(False)
plt.show()
此时,一旦有了可用的联邦数据,因为这是一个模拟环境,所以可以做几个测试来检查客户机数据的非 iid 行为的强度。我将把这个练习留给您去探索和试验,但是请记住,在现实世界中,这种类型的分析是不可能的,因为数据不是集中可用的。
预处理输入数据
对于预处理,您必须确保以下几点:
-
数据质量
-
通过将像素强度乘以 1/255 来重新调整所有通道的比例,您已经标准化了像素值。
-
合适的比例:装载时已经保证
-
增强以创建更多数据并避免过拟合 OOB/验证数据集。因为案例研究是为了发现联邦原则,所以您现在将跳过这一步,回到第八章。
-
-
培训改进:
-
使用梯度下降创建训练批次
-
随机产生随机性,使损失与样本选择无关
-
预取某些样本以减少训练滞后的可能性,因为样本用于训练,所以您必须运行预处理
-
SHUFFLE_BUFFER = len(sample_dataset) # How much data to shuffle
EPOCHS = 5 # Number of epochs to run for training @ individual node
PREFETCH_BUFFER = 100 # Preloading some number of samples to aid faster training.
# Normalizing the pixel values
normalization_layer = tf.keras.layers.experimental.preprocessing.Rescaling(1.0/255)
def preprocess(dataset):
def batch(sample):
_x = normalization_layer(sample['pixels'])
return collections.OrderedDict(
x = _x
y = tf.reshape(sample['label'], [-1, 1]))
return dataset.repeat(EPOCHS).shuffle(SHUFFLE_BUFFER).batch(
BATCH_SIZE).map(batch).prefetch(PREFETCH_BUFFER)
创建联合数据
既然已经准备好了预处理函数,就可以通过创建客户机数据集的迭代器来创建最终的联邦数据。
此外,在现实世界的设置中,您通常从大量客户端中选择一个客户端样本,因为其中只有一小部分可用(跨设备设置)。
selected_clients = np.random.choice(train_dataset.client_ids,NUM_CLIENTS, replace = False)
federated_train_data = (preprocess(train_dataset.create_tf_dataset_for_client(i)) for i in selected_clients)
您还可以使用前面创建的示例批处理来创建预处理联邦数据集的示例,因为它可以在以后用于输入规范。
sample_federated_dataset = preprocess(sample_dataset)
联合通信
在 TFF 框架内,任何在本地训练的模型都需要包装在tff.learning.Model
接口中。这允许两件事:
-
帮助计算各个节点的联合指标和性能
-
一组变量在每个本地节点上的筒仓中受到影响。
首先创建一个训练模型函数,它构建您正在使用的 NNet 体系结构。
-
Conv2D 层对输入图像进行卷积运算,以捕捉每个像素的局部效应。
-
池层是为了降低维度,集中信息。
-
为了防止过度适应,在训练过程中你会去掉随机的神经元。
-
最后,在展平用于预测的下降图层的 2-D 输出后,添加一个密集图层。
def train_model():
model = Sequential([
tf.keras.layers.InputLayer(input_shape=(IMG_HEIGHT,IMG_WIDTH, 3)),
# Ingesting a 2-d Image with 3 channels
tf.keras.layers.Conv2D(16, 3, padding='same', activation='relu'),
# Max pooling to reduce dimensions
tf.keras.layers.MaxPooling2D(),
tf.keras.layers.Conv2D(32, 3, padding='same', activation='relu'),
tf.keras.layers.MaxPooling2D(),
tf.keras.layers.Conv2D(64, 3, padding='same', activation='relu'),
tf.keras.layers.MaxPooling2D(),
# Dropout to prevent over-fitting
tf.keras.layers.Dropout(0.2),
# Flattening to feed data for sigmoid activation
tf.keras.layers.Flatten(),
tf.keras.layers.Dense(128, activation='relu'),
tf.keras.layers.Dense(len(class_names)-1, activation = 'sigmoid')
])
return model
def federated_train_model():
local_train_model = train_model()
return tff.learning.from_keras_model(
local_train_model,
input_spec=sample_federated_dataset.element_spec,
loss=tf.keras.losses.BinaryCrossentropy(),
metrics=[tf.keras.metrics.AUC()])
接下来,您将为中央服务器创建流程,以便使用来自所有本地节点的参数更新来更新中央模型。
parameter_iteration_process = tff.learning.build_federated_averaging_process(
federated_train_model,
client_optimizer_fn = lambda: tf.keras.optimizers.SGD(learning_rate= CLIENT_LR),
server_optimizer_fn = lambda: tf.keras.optimizers.SGD(learning_rate= SERVER_LR))
TFF 构建了一对联邦计算,并将它们打包到一个tff.templates.IterativeProcess
中,这些计算作为一对称为initialize
和next
的属性可用。
-
initialize
表示服务器上联合平均过程的状态。它包括-
模型:分配给所有设备的初始参数
-
优化器状态:为联邦指标计算和平均而维护。它跟踪梯度更新。
-
三角洲聚集物
-
-
next_fn
将利用client_update
和server_update
,代表一个联合平均周期。
state = parameter_iteration_process.initialize()
state, metrics = parameter_iteration_process.next(state, federated_train_data)
print('round 1, metrics={}'.format(metrics))
输出
round 1, metrics=OrderedDict([('broadcast', ()), ('aggregation', OrderedDict([('value_sum_process', ()), ('weight_sum_process', ())])), ('train', OrderedDict([('auc', 0.5897039), ('loss', 0.6823319)]))])
同样,你可以有多个回合。
NUM_ROUNDS = 6 # Total 5 rounds of training
for round_num in range(2, NUM_ROUNDS):
state, metrics = parameter_iteration_process.next(state, federated_train_data)
print('round {:2d}, metrics={}'.format(round_num, metrics))
输出
round 2, metrics=OrderedDict([('broadcast', ()), ('aggregation', OrderedDict([('value_sum_process', ()), ('weight_sum_process', ())])), ('train', OrderedDict([('auc', 0.60388386), ('loss', 0.67804503)]))])
round 3, metrics=OrderedDict([('broadcast', ()), ('aggregation', OrderedDict([('value_sum_process', ()), ('weight_sum_process', ())])), ('train', OrderedDict([('auc', 0.61434853), ('loss', 0.6752475)]))])
round 4, metrics=OrderedDict([('broadcast', ()), ('aggregation', OrderedDict([('value_sum_process', ()), ('weight_sum_process', ())])), ('train', OrderedDict([('auc', 0.62443274), ('loss', 0.67076266)]))])
round 5, metrics=OrderedDict([('broadcast', ()), ('aggregation', OrderedDict([('value_sum_process', ()), ('weight_sum_process', ())])), ('train', OrderedDict([('auc', 0.6333971), ('loss', 0.6674127)]))])
你们中的一些人可能会发现训练过程(收敛)有点慢。实际上,这是由于较低的服务器学习率。我把它保持在 0.1。如果将它保持为 1,这意味着每次迭代都全力贡献给中心模型的参数。换句话说,更新是完全学习的。
Note
如果您在 Jupyter 笔记本中运行相同的代码,您必须允许异步操作。在 Python 中,可以通过调用
import nest_asyncio
nest_asyncio.apply()
估价
TensorFlow 库提供了build_federated_evaluation,
,它允许通过联邦通信(跨本地节点)来聚合指标。
def evaluate(train_fn, state, train_data, test_data):
# Print training metrics
evaluation = tff.learning.build_federated_evaluation(train_fn)
train_metrics = evaluation(state.model, train_data)
print("Training Metrics: AUC : {}, Binary Cross Entropy Loss: {}".format(
train_metrics['auc'],
train_metrics['loss']))
# Print testing metrics
test_metrics = evaluation(state.model, test_data)
print("Validation Metrics: AUC: {}, Binary Cross Entropy Loss: {}".format(
test_metrics['auc'],
test_metrics['loss']))
您必须通过与训练数据相同格式的验证集。为此,您需要创建一个client_test_dataset
,它是一个字典,包含每个本地节点或服务器节点的验证数据。
然后使用上面定义的preprocess()
函数对所有的验证进行评估处理。
val_dataset = tff.simulation.FromTensorSlicesClientData(client_test_dataset)
federated_val_data = [preprocess(val_dataset.create_tf_dataset_for_client(i)) for i in selected_clients]
evaluate(federated_train_model, state, federated_train_data, federated_val_data)
输出
Training Metrics: AUC : 0.6697379946708679, Binary Cross Entropy Loss: 0.6773737072944641
Validation Metrics: AUC: 0.6535744071006775, Binary Cross Entropy Loss: 0.6790395379066467
在这一节中,我讨论了 TF learning API。TFF 还提供了核心 API,您可以在其中修改 TFF 提供的几个不同组件,如联合平均技术和联合通信(跨设备网络负载和本地处理)。
结论
联邦学习是一个不断发展的领域,随着保护私有和昂贵数据的需求变得普遍,它一定会发展壮大。在这一章中,您详细介绍了差分隐私和多方通信的隐私机制,但是新的研究还在不断出现。秦斌等人的“联邦学习系统调查:数据隐私和保护的愿景、宣传和现实”是一篇出色的论文,揭示了联邦学习的不同层面。
话虽如此,联邦学习并不是进行受保护学习的唯一方式。人们也在研究没有中央服务器来协调工作的对等系统;相反,它是自治的。这种系统在现实世界中的可靠性还有待证实。
最后,Owkin、Google 和 Apple 等几家公司正在积极投资联合技术,特别是围绕患者的药物发现、打字推荐和改进聊天机器人。在我看来,ML 产品进入市场解决国家间本地问题的速度意味着它是一项重要的技术。*
八、医学成像
医学图像分析在过去的三十年里有了巨大的发展。最初,该领域的分析被视为应用模式识别和精算计算机视觉方法,但随着高级图像处理和基于深度学习的方法的广泛使用,该领域不仅在算法进步方面,而且在处理各种数据方面都发展迅速,因为在此期间出现了不同的模式。
在这个案例研究中,你将接触到医学成像的许多不同方面。您将特别关注不同类型的医疗数据,以及这些医疗图像数据是如何捕获、数字化存储和分发的。你将不会触及这些图像是如何基于组织-能量相互作用和相关统计数据而形成的物理学。
您将使用二维和三维图像深入研究图像分割和分类的两个终端应用。最后,您将探索当前存在的各种挑战,如图像质量、可解释性和对抗性攻击。
什么是医学影像?
医学成像涉及对不同成像模式(例如 X 射线、CT、MRI 等)上的生物医学图像的科学分析。监测健康(通过筛查)、疾病和伤害的诊断和治疗。
这些生物医学图像是在宏观、中观和微观等不同尺度上对人体、器官或组织的测量。这些标尺的穿透深度和图像分辨率不同,如图 8-1 和 8-2 所示。
图 8-2
基于尺度的光学成像技术比较。资料来源:Subhamoy Mandal 等人,“将生物成像扩展到第五维度”
图 8-1
光学分辨技术综述。资料来源:光学学会
生物医学图像来源于测量人体不同物理特性的不同成像模式。
图像模态
图像模态是通过利用与技术/设备中使用的能量类型的相互作用,以 n 维图像的形式捕捉器官/组织特征的各种方式。举个例子,
-
X 射线成像中的辐射吸收
-
超声波中的声压
-
磁共振成像中的射频信号幅度
MRI、超声、X 射线和 CT 是一些主要的成像方式,但还有更多,如图 2 所示。存在如此多的模态是因为一个简单的原因,即单一的技术不足以捕捉人体解剖和生理。
为了提供这些技术的简要概述,表 8-1 对主要模态进行了比较和对比。
表 8-1
比较不同的图像形态
|南号码
|
形式
|
应用
|
主要特点
|
缺点
|
辐射
|
| — | — | — | — | — | — |
| one | x 射线 | 非均匀组成的材料,如骨骼。这些图像有助于评估是否存在疾病、损伤或异物。 | 获得图像通过使用 x 光片。无创无痛。 | 有时结构重叠会造成翻译上的问题。 | 电离 |
| Two | 计算机化 X 线体层照相术 | 非均匀组成的材料,如骨骼。这些图像有助于评估是否存在疾病、损伤或异物。 | 扫描是使用 x 光和后来的计算机被用来构建一个系列横截面的图像。这就消除了叠加。 | 高剂量的电离辐射,因此将来会引起致癌的疾病。 | 电离 |
| three | 核磁共振成像 | 一般用于分析撕裂的韧带和肿瘤。也有助于检查大脑和脊髓。 | 使用磁信号和无线电波。 | 强信号会导致幽闭恐惧症倾向。 | 非电离的 |
| four | 超声 | 主要是胎儿成像。也用于腹部器官、心脏、乳房、肌肉、肌腱、动脉和静脉的成像。 | 使用高频率声音信号对诸如器官的内部结构成像,软组织和未出生的婴儿。 | 容易产生噪声,并且该过程由放射科医师驱动,因此容易出现人为错误。 | 非电离的 |
那么,我们为什么会对理解这些模式感兴趣呢?
首先,要理解依赖于手头的用例,我们必须仔细选择要使用的模态。
图 8-3
核磁共振扫描比 CT 更清楚地显示受伤的脑组织
-
提高发现问题的敏感性(异物/血管问题等)。)
-
与像 X 射线这样的 2-D 成像模态相比,3-D 成像模态允许更好的定位。
-
更好地描绘组织类型。例如,如图 8-3 所示,如果目标是发现中风中受伤的脑组织,您可以看到,与 CT 相比,MRI 图像清楚地显示了受损区域,而 CT 的大部分区域是黑暗的。
其次,这些模式在获取价值的方式上可能有所不同。因为数字图像有像素强度,所以有不同的度量来测量数字医学图像中的信息值。
CT 扫描和 X 射线 Hounsfield 单位(HU)用于测量电离辐射的强度。更高的 HU 意味着辐射更难通过,因此有更高的衰减。参见图 8-4 。
图 8-4
Hounsfield 标度范围从-1000 到+ 1000。资料来源:Osborne 等人,www . southsudanmedica l journal . com/
注意许多灰色的阴影。人眼甚至计算机在某些情况下都不可能在如此小的梯度上工作,因此使用一种称为开窗的技术来观察感兴趣的区域,如软组织、肺和骨骼。确定窗位 L 和宽度 W。然后,仅在范围L-W/2 到 L + w /2 内保持梯度,其余部分完全变成黑色(小于 L - W/2)和白色(大于 L + W/2)。所有这些重要的决定都是在我们对这些图像建模之前做出的。
最后,不同的模式可以使用不同的造影剂来突出某些组织区域。由于不同组织对药物的吸收率不同,某些组织就显得特别突出。CT 成像使用碘基,而磁共振成像使用钆基,通常口服(如片剂)或静脉注射(直接泵入血流)。如图 8-5 所示,由于使用了对比剂,经过一段时间(20-30 秒延迟)后,您可以看到致癌结节在肝脏中突出显示。
图 8-5
造影剂突出某些组织区域
数据存储
在我们深入探讨如何处理多维图像格式并在其上建立机器学习模型之前,让我们快速了解一下医学图像数据将采用的标准文件格式以及不同的组成部分。
典型的医学图像由四个基本部分组成:
-
像素深度
-
像素数据
-
[计]元数据
-
光度解释
我们用一个很简单的例子来理解这一点。假设你有一张黑白图像,其中有各种不同的灰度强度,并且你知道各种其他信息,如这张图像的人物、内容和时间。
根据此描述,各种元素对应于您刚刚了解到的基本组件。
-
黑白图像告诉我们所使用的通道以及图像的光度解释。图像可以是单色的,也可以是彩色的。
-
不同的强度等级暗示了两件事。首先,它告诉我们像素深度,即用于编码信息的位数。例如,一个 8 位像素可以包括 28 个亮度级(无符号为 0 到 255)。其次,它告诉我们像素强度值以及像素值的范围。
-
其他信息是元数据,如研究日期、图像模态、患者性别、形状等。
特殊数据需要特殊格式。从射线照相设备收集的图像主要有六种不同的格式。
-
医学数字成像和通信
-
神经成像信息学技术倡议
-
PAR/REC (Philips MRI 扫描仪格式)
-
分析(梅奥医学成像)
-
NRRD(接近原始栅格数据)
在这六种不同的格式中,DICOM 和 NIFTI 是使用最广泛的。DICOM 和 NIFTI 之间的主要区别在于,DICOM 中的原始图像数据存储为二维切片文件的集合,这使得该结构对于三维数据分析来说有点麻烦,而在 NIFTI 中,我们有完整的三维图像。
除了处理、存储、打印和传输信息之外,所有这些不同的格式还可以帮助您获得您的机器学习模型所需的所有功能,作为一名数据科学家,这是您应该集中精力的地方:哪种格式可以提供哪种信息。
在前面的小节中,为了让您熟悉这些格式,您将使用 DICOM 的一个例子和 NIFTI 的另一个例子。
处理二维和三维图像
在我们讨论过的各种方法中,只有 X 射线成像能产生二维图像。CT 和 MRI 创建三维图像,因为它们从不同的角度捕捉信息。还有另外两种形式可以帮助我们捕捉二维数据:
-
眼底成像:用于扫描眼部微小血管的健康状况。通常用于识别糖尿病视网膜病变(DR)。
-
病理成像:通过对细胞进行染色,使不同的细胞结构呈现不同的颜色,然后进行数字化处理,从而获得的细胞水平成像(记得上一章)。
然而,在本章中,我们将讨论通过 X 射线设备进行二维图像分析。
类似地,3d 图像不限于 CT 和 MRI。其他形式如超声波和 PET/SPECT 扫描也能产生 3d 图像,用于了解人体的不同部位。我们可以将 3d 图像视为 2d 图像的堆叠,使得这些图像从不同角度拍摄,然后拼接在一起以创建全面的 3d 视图。
在这一章中,我们将讨论通过 MRI 设备进行三维图像分析。
你可能有时会听到术语四维图像。好吧,别惊讶。它只是在不同的时间或不同的子模型中捕获的几个三维图像,就像 MRI ??、?? 等的情况。我不会报道它,但如果你有兴趣,我建议你看看李等人的论文,题为“医学影像和放射治疗的进展”
处理二维图像
你将接受北美放射学会的 RSNA 肺炎检测挑战。它组织在 Kaggle 上,拥有 DICOM 格式的数据。尽管该挑战旨在定位胸片上的肺活量,但 DICOM 元数据文件还包含以下图像标签:
-
标准
-
没有肺部阴影/不正常
-
肺部阴影
因此,您也将使用相同的标签进行影像分类。
RSNA 每年都组织医学影像竞赛。查看不同数据集和比赛的空间: www.rsna.org/education/ai-resources-and-training/ai-image-challenge
。
Python 中的 DICOM
您的目录应该如下所示:
-
从 Kaggle 笔记本的数据部分下载单独的 ZIP 文件,并按照上面所示的格式创建目录。
-
stage_2_detailed_class_info.csv
包含每个患者 id 的标签,而train
和test
文件夹中的每个 DICOM 文件被命名为patient-id
。 -
stage_2_train_labels.csv
包含来自train
和test
文件夹的每个患者 id 的目标标签肺炎或非肺炎。
Data/
...2d_lung_opacity_challenge/
......Train/
.........000db696-cf54-4385-b10b-6b16fbb3f985.dcm
.........000fe35a-2649-43d4-b027-e67796d412e0.dcm
......Test/
.........00b4e593-fcf8-488c-ae55-751034e26f16.dcm
.........00f376d8-24a0-45b4-a2fa-fef47e2f9f9e.dcm
......stage_2_detailed_class_info.csv
DICOM 文件包含标题元数据和原始图像像素阵列的组合。在 Python 中,你可以使用一个名为pydicom
的库来处理 DICOM 文件。
还记得我们讨论过除了典型的图像分析流水线之外,不同的图像模态如何引入新的预处理步骤吗?幸运的是,在您的案例中,您已经有了预处理数据。RSNA 共享的数据在两个方面进行预处理:
-
将高动态范围转换为 8 位编码,其值范围为 0 到 255 灰度。
-
图像通常以较高的分辨率捕获,但出于实用目的,图像被调整为 1024 x 1024 矩阵。
对于那些还在考虑如何进行窗口和大小调整的人来说,如果这样的预处理还没有完成的话,这里有一些代码:
def windowed_image(img, center, width):
img_min = center - width // 2
img_max = center + width // 2
windowed_image = img.copy()
windowed_image[windowed_image < img_min] = img_min
windowed_image[windowed_image > img_max] = img_max
return windowed_image
您可以合并目标和类数据,以便更好地理解分布。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import pydicom
import glob
import os
BASE_DIR = "./Data"
DATA_DIR = os.path.join(BASE_DIR,"2d_lung_opacity_challenge/")
classes = pd.read_csv(glob.glob(os.path.join(DATA_DIR,"*.csv"))[0])
target = pd.read_csv(glob.glob(os.path.join(DATA_DIR,"*.csv"))[1])
train_labels = pd.merge(classes, target[["patientId","Target"]], on = "patientId", how="left")
因为有多个patientId
,所以单个图像可以有多个边界框。扔掉重复的。
assert train_labels.drop_duplicates().shape == train_labels.drop_duplicates('patientId').shape:
train_labels = train_labels.drop_duplicates().reset_index(drop = True)
print(train_labels.groupby(['class', 'Target']).size().reset_index(name='Patient Count').to_markdown())
| | class | Target | Patient Count |
|---:|:-----------------------------|---------:|----------------:|
| 0 | Lung Opacity | 1 | 6012 |
| 1 | No Lung Opacity / Not Normal | 0 | 11821 |
| 2 | Normal | 0 | 8851 |
哪里有肺部阴影,哪里就有肺炎。然而,医学上肺部阴影不能完全和唯一地确定肺炎,因为诊断需要其他临床信息,如实验室数据、症状等。但为了简单起见,所有肺部阴影都称为肺炎。然而,在现实世界中,你不能做同样的假设;你必须咨询适当的医学研究人员和放射科医生来做出这样的假设。
接下来,非肺炎可以分为无肺部阴影/不正常和正常。正常的图像是健康胸部的图像。你不能说没有肺部阴影/不正常。让我们来看看其中的几个。
def draw(input_ids):
# A maximum of 3 images in a row
ncols, nrows = min(3,len(input_ids)), len(input_ids)//min(3,len(input_ids)) +1 if len(input_ids)%min(3,len(input_ids)) !=0 else len(input_ids)//min(3,len(input_ids))
# figure size, inches
figsize = [10, 8]
# create figure (fig), and array of axes (ax)
fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=figsize)
# plot image for single sub-plot
for i, axi in enumerate(ax.flat):
try:
dicom_path = input_ids[i]
data = pydicom.read_file(dicom_path)
# one can also use plt.cm.bone
axi.imshow(data.pixel_array, cmap="gray")
# get indices of row/column
rowid = i // ncols
colid = i % ncols
except IndexError as e:
continue
# For some of you who want to add bounding box info to plots as
# well can access by row-id and col-id on the array of axes
# ax[row-id][col-id].plot()
plt.tight_layout(True)
plt.show()
np.random.seed(123)
examples_non_normal = np.random.choice(train_labels[train_labels["class"].\
isin(["No Lung Opacity / Not Normal"])].patientId,
size = 3,
replace = False)
examples_non_normal = [os.path.join(DATA_DIR,"Train",x+".dcm") for x in examples_non_normal]
draw(examples_non_normal)
图 8-6 中的一些观察结果。
-
肺部阴影图像和无肺部阴影/不正常图像具有一些相似的特征。
-
电线和管子的存在,这表明除了目标值为 0 的肺炎之外,可能还有其他一些观察到的疾病。
-
在大多数情况下,间隙/不透明性(充满液体/病原体的间隙等)的性质。)对于两种类型是不同的,尽管由于类似于 COPD 或哮喘的外来物质在肺中的扩散,它可以重叠。
-
-
由于胸腔积液,液体或异物的积聚可以
-
渗出来让肺看起来更小。参见无肺部阴影/不正常行中的样本 3。
-
上述情况很容易与肺部阴影混淆,因此在这种情况下可能需要咨询几位放射科医生才能得出结论。
-
图 8-6
三种不同标签的样品
对类别标签进行分析而不仅仅是盲目跟踪目标标签的目的是让您意识到,医学图像分析需要一定的领域知识来理解和实现一个强大的图像分析系统。尤其是如果你计划将模型投入使用,FDA 将会调查与你的模型相关的风险,在这种情况下,这种微妙的理解将会派上用场。
基于 DICOM 元数据的 EDA
您可以定义一个函数来从 DICOM 文件中选择重要的元数据。
-
患者年龄
-
病人性别:只有两类,男性和女性
-
像素间距:较高的像素间距意味着图像质量较低
-
查看位置:AP(射线从胸部到背部,卧位;通常用于病人或老年人)和 PA(射线从背部到胸部,站立位置)
def get_metadata(patient_id):
"""
Returns metadata from each dicom file
"""
data = pydicom.read_file(os.path.join(DATA_DIR,"Train",patient_id+".dcm"),
stop_before_pixels=False)
_id = data.PatientID
_age = data.PatientAge
_sex = data.PatientSex
# col_spacing (horizontal)
_pixelspacing_x = data.PixelSpacing[1]
# row_spacing (vertical)
_pixelspacing_y = data.PixelSpacing[0]
_viewpos = data.ViewPosition
_mean = np.mean(data.pixel_array)
_min = np.min(data.pixel_array)
_max = np.max(data.pixel_array)
return pd.DataFrame([[_id, _age, _sex, _pixelspacing_x, _pixelspacing_y, _viewpos ,_min, _max, _mean]],
columns = ["patientId","age","sex","pixel_spacing_x","pixel_spacing_y","view_pos",
"min_pixint","max_pixint","mean_pixint"])
使用并行处理,您可以捕获所有元数据,以查看其相关性和对目标变量的影响。
from joblib import Parallel, delayed, parallel_backend
from tqdm import tqdm
train_dicom = Parallel(n_jobs=os.cpu_count()-1, backend="threading")(delayed(get_metadata)(pt_id) for pt_id in tqdm(train_labels.patientId))
然后连接从每个 DICOM 返回的单个数据帧。
train_dicom_df = pd.concat(train_dicom, axis = 0)
最后,将目标/标签数据集与元数据数据帧合并,并创建用于分析的数据。参见图 8-7 。
图 8-7
患者在不同观察位置的分布
# Train Labels with Metadata
train_labels_w_md = pd.merge(train_labels, train_dicom_df, on = "patientId", how="left")
查看位置
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
sns.countplot(x='view_pos', hue='class', data=train_labels_w_md, ax=axes[0])
sns.countplot(x='view_pos', hue='Target', data=train_labels_w_md, ax=axes[1])
基于以下原因,视图位置看起来是一个重要的变量:
-
虽然 PA 或 AP 位置的患者数量相似,但 AP 位置的肺炎患者更多(Target = 1)。
-
此外,对于 AP 视图,肺部不透明度标记很明显,而对于 PA 视图,正常标记更明显。
年龄
根据目标和类别标签绘制年龄分布图,以检查分布情况。
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
p = sns.distplot(train_labels_w_md[train_labels_w_md['class']=='No Lung Opacity / Not Normal']['age'],
hist=True,
kde=False,
color='red',
label='No Lung Opacity / Not Normal', ax=axes[0])
p = sns.distplot(train_labels_w_md[train_labels_w_md['class']=='Normal']['age'],
hist=True,
kde=False,
color='cornflowerblue',
label='Normal', ax=axes[0])
p = sns.distplot(train_labels_w_md[train_labels_w_md['class']=='Lung Opacity']['age'],
hist=True,
kde=False,
color='lime',
label='Lung Opacity', ax=axes[0])
_ = p.legend()
p = sns.distplot(train_labels_w_md[train_labels_w_md['Target']==0]['age'],
hist=True,
kde=False,
color='gray',
label='0', ax=axes[1])
p = sns.distplot(train_labels_w_md[train_labels_w_md['Target']==1]['age'],
hist=True,
kde=False,
color='lime',
label='1', ax=axes[1])
_ = p.legend()
正如您在图 8-8 中看到的,年龄类别对于任何目标(0 或 1)都没有显示出任何明显的特征,对于类别标签也是如此。然而,可以形成间距为 20 的某些组。
图 8-8
患者的年龄分布
性
你在性专栏做了两项分析。
-
检查患者在不同目标和不同性别标签中的分布。
-
Although the distribution is not structurally different, women generally have a higher pneumonia percentage (see Figure 8-9).
图 8-9
患者在性别栏中的分布
-
-
年龄和性别相关性
- 两性显示相同的分布模式,正常患者的平均年龄通常低于无肺部阴影/不正常患者(见图 8-10 )。
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
sns.countplot(x='sex', hue='class', data=train_labels_w_md, ax=axes[0])
sns.countplot(x='sex', hue='Target', data=train_labels_w_md, ax=axes[1])
train_labels_w_md["age"] = train_labels_w_md.age.apply(lambda x:int(x))
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
sns.boxplot(x='sex', y = 'age', hue='class', data=train_labels_w_md, ax=axes[0])
sns.boxplot(x='sex', y = 'age', hue='Target', data=train_labels_w_md, ax=axes[1])
图 8-10
用箱线图来观察不同性别的年龄差异
像素间距
像素间距代表每个像素的大小。每个像素代表图像上的某个区域。不同的像素间距会导致空间信息的不均匀分布。让我们看看这种差异有多明显。
首先,将像素间距四舍五入两点。
train_labels_w_md["pixel_spacing_x_norm"] = train_labels_w_md.pixel_spacing_x.apply(lambda x: round(float(x),2))
train_labels_w_md["pixel_spacing_y_norm"] = train_labels_w_md.pixel_spacing_y.apply(lambda x: round(float(x),2))
接下来,绘制患者计数。
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
plot = sns.countplot(x='pixel_spacing_x_norm', hue='class', data=train_labels_w_md, ax=axes[0][0])
plot.set_xticklabels([x for x in plot.get_xticklabels()],rotation=90)
plot.legend(loc='upper right')
plot = sns.countplot(x='pixel_spacing_x_norm', hue='Target', data=train_labels_w_md, ax=axes[0][1])
plot.set_xticklabels([x for x in plot.get_xticklabels()],rotation=90)
plot.legend(loc='upper right')
sns.countplot(x='pixel_spacing_y_norm', hue='class', data=train_labels_w_md, ax=axes[1][0])
plot.set_xticklabels([x for x in plot.get_xticklabels()],rotation=90)
plot.legend(loc='upper right')
sns.countplot(x='pixel_spacing_y_norm', hue='Target', data=train_labels_w_md, ax=axes[1][1])
plot.set_xticklabels([x for x in plot.get_xticklabels()],rotation=90)
plot.legend(loc='upper right')
图 8-11 中有两个主要观察结果:
图 8-11
按像素间距划分的患者分布
-
像素间距有很大的变化,范围从 0.13 到 0.2。
-
对于 0.15 至 0.17 之间的间距,更容易观察到肺部阴影。
尽管您不会直接使用像素间距,但您将使空间信息一致,为此,您会将图像重新采样为 1 毫米 X 1 毫米。
平均强度
最后,您将看到标签和目标的平均强度是如何变化的。如果您发现目标或标签有不同的模式,您肯定会尝试包含元数据
在图 8-12 中,肺炎和非肺炎患者遵循相似的分布,因此他们没有给出任何类别的任何特殊信息。
图 8-12
平均强度下的患者分布
Note
同样,您可以使用 DICOM 元数据进行更多的分析。我希望您对处理 DICOM 文件和元数据有了相当详细的了解。
处理三维图像
处理三维图像与处理二维图像没有太大的不同,但是为了让事情变得更有挑战性,让我们学习一下 NIFTI 中的三维图像处理。与上一节一样,这将帮助您为前面的图像分割案例研究做好准备。
您将使用位于 www.med.upenn.edu/cbica/brats2020/data.html
的佩雷尔曼医学院的 BRATS 2020 数据集。BRATS 2020 数据包含各种模式的 NIFTI 文件(通过改变 MRI 机器中的脉冲序列制作),即
-
本地(??)
-
增强后 ?? 加权(T1CE)
-
?? 加权(??)
-
?? 流体衰减反演恢复(??-FLAIR)
如果你有兴趣了解脉冲序列之间的差异,可以访问 https://radiopaedia.org/articles/mri-pulse-sequences-1?lang=us
了解更多。
既然您的数据将是三维数据,那么让我们来理解是什么使它成为三维数据。在医疗系统中,我们的身体可以分为三个平面:
-
轴向/横向:从上到下
-
矢状面:身体从左到右
-
冠状:从后到前(从后到前)
这就是核磁共振成像的三维图像。在上面的 DICOM 图像中,图像仅在冠状面拍摄,因此您看到的是二维图像。要了解更多关于这个主题的信息,请参考 https://teachmeanatomy.info/the-basics/anatomical-terminology/planes/
。
尽管 MRI 数据来自不同的临床协议和多台扫描仪,但它们已经以三种方式进行了预处理:
图 8-13
医学图像的配准
- 与相同的解剖模板共同配准:由于我们正在采集多种模式的 MRI,如果患者在采集这些图像之间移动(即使是轻微的移动),当所有序列组合在一起用于分割任务时,可能会导致未对准,因此需要进行一个称为配准的过程来避免这种错误。既然这个已经给你做好了,就不需要你操心了。见图 8-13 。
图 8-14
去除头骨的图像示例。(a)核磁共振成像;(b)颅骨剥离 MRI 图像。来源:Arun 等人在 2017 年发表的“SVM-LWT 实现了基于模糊聚类的脑肿瘤检测图像分析”
-
以相同的分辨率进行插值:这仅仅意味着所有四个序列的空间信息在整个 3d 体积上是一致的。
-
颅骨剥离:在 MRI 图像中,当解决诸如脑肿瘤分割等任务时,去除颅骨边界是一个好的做法,原因很简单,颅骨边界不提供任何有助于解决分割问题的信息,因此我们只是将其剥离。见图 8-14 。
Note
如果明天你开始使用一个不同的数据集,确保你检查了这三件事。
nifti 格式
NIFTI 文件格式不是由扫描仪生成的,因此元数据信息不如 DICOM 文件格式丰富,但它仍然有一些元数据。此外,它将图像系列表示为单个文件。
与 DICOM 相比,NIFTI 文件中的坐标系略有不同。了解这种差异是有好处的,因为通常您可能希望以 NIFTI 格式保存文件,因为整个图像系列都存在于单个文件中,这与 DICOM 不同,因此更易于共享和维护。参见图 8-15 。
图 8-15
NIFTI 图像中的坐标系
最后,NIFTI 中的测量单位可以不同,不像 DICOM 那样固定为毫米。NIFTI 单独存储测量单位的信息(例如,您在 DICOM 中看到的像素间距信息)。
因为与 DICOM 相比,NIFTI 头没有那么复杂,所以实际上我们很少使用 NIFTI 头信息。下面是一些需要注意的重要信息:
-
像素间距
-
三个平面的尺寸
-
XYZ-T 单位
磁共振图像处理简介
让我们快速设置您的输入流水线。
BASE_DIR = "./Data/3d_brain_tumor_segmentation/MICCAI_BraTS2020_TrainingData/"
label_paths = glob.glob(os.path.join(BASE_DIR,"**","*seg.nii"))
flair_paths = glob.glob(os.path.join(BASE_DIR,"**","*flair.nii"))
t1_paths = glob.glob(os.path.join(BASE_DIR,"**","*??.nii"))
t1ce_paths = glob.glob(os.path.join(BASE_DIR,"**","*t1ce.nii"))
t2_paths = glob.glob(os.path.join(BASE_DIR,"**","*??.nii"))
# Let's create a dictionary of dictionary to order the data
full_data = {i:{'label':label,
'flair':flair,
'??':??,
't1ce':t1ce,
'??':??} for i, (label,flair,??,t1ce,??) \
in enumerate(zip(label_paths,
flair_paths,
t1_paths,
t1ce_paths,
t2_paths))}
你已经知道有四个不同的序列,每个序列都可以在轴向、矢状和冠状面上以三种不同的方式观察。
patient_id = 5
k=1
plt.figure(figsize=(20,20))
for i,seq in enumerate(["flair","??","t1ce","??"]):
img = io.imread(full_data[patient_id][seq], plugin='simpleitk')
for j in range(3):
if (j==0):
plt.subplot(4,3,k)
plt.imshow(img[100,:,:])
# x-y plane
plt.title("Axial/Traverse View")
plt.ylabel(seq.upper())
k=k+1
elif (j==1):
plt.subplot(4,3,k)
plt.imshow(img[:,100,:])
plt.title("Coronal View")
k+=1
else:
plt.subplot(4,3,k)
plt.imshow(img[:,:,100])
plt.title("Sagittal View")
k+=1
在图 8-16 中,您可以清楚地看到不同的模态如何在不同的视图中突出大脑的不同部分,并提供补充信息。
图 8-16
磁共振成像模式和视图的交叉
对于同一个患者,我们也来看看目标标签。
# For the same patient let's also have a look at the target label
img = io.imread(full_data[patient_id]['label'], plugin='simpleitk')
plt.figure(figsize = (20,20))
k = 1
for i in [50, 75,100, 125]:
plt.subplot(1,4,k)
plt.imshow(img[i,:,:])
plt.title("Labels:- " + ", ".join([str(i) for i in np.unique(img[i,:,:])]))
k+=1
如图 8-17 所示,分段任务有四种不同的标签。肿瘤部分用绿色、黄色和蓝色(1、2 和 4)标记,而背景用紫色(0)标记
图 8-17
分段标签(轴向视图)
从图 8-16 和 8-17 中的上述曲线图可以清楚地看出
-
不是所有的切片都重要。
-
对于不同的序列,像素强度不是均匀分布的。
-
分割标签的像素强度高度不平衡(因为图像的大部分是紫色,其次是黄色标签 2)。
让我们探讨一下上面的第 1 点和第 2 点,看看你从这些观察中得到了什么。
非均匀像素分布
让我们快速看看不同的序列是如何使像素强度发生变化的。这将帮助您决定每个序列的规范化策略。
import seaborn as sns
sns.set_style('whitegrid')
fig, axes = plt.subplots(nrows=4, ncols=4)
fig.tight_layout(pad=1, w_pad=1, h_pad=0.5)
fig.set_size_inches(20,20)
k=1
for patient_id in [5,10,20,50]:
for i,seq in enumerate(["flair","??","t1ce","??"]):
img = io.imread(full_data[patient_id][seq], plugin='simpleitk')
if (i==0):
plt.subplot(4,4,k)
plt.hist(x= img.reshape(-1,1))
plt.title(seq.upper())
plt.ylabel("Patient "+str(patient_id))
k=k+1
else:
plt.subplot(4,4,k)
plt.hist(x= img.reshape(-1,1))
plt.title(seq.upper())
k+=1
从图 8-18 中可以清楚地看到
图 8-18
患者不同磁共振成像序列的强度变化
-
大多数像素的亮度为 0,并且是右偏的。
-
对于不同的序列,对于异常值处理可以观察到不同的截止值。对于 ?? 和 T1CE,它大约是 500,而对于 FLAIR 和 ??,它从 300 到 600 不等。
-
您必须对此进行规范化,并处理偏斜。
相关性检验
为了分析您是否需要考虑所有的切片,您可以跨深度维度进行相关性测试。您将遍历图像的深度并计算相关性。为了使您的工作更容易,您只需将它转换成熊猫数据帧,然后计算相关性。其思想是,如果相邻或接近相邻切片的像素强度没有变化,它们将产生 NA 形式的相关性,因为它们的协方差为 0。
让我们快速绘制一些图表,看看哪些切片是不相关的。
k = 1
from itertools import chain
fig, axes = plt.subplots(nrows=1, ncols=4)
fig.tight_layout(pad=0.5, w_pad=2, h_pad=0.5)
fig.set_size_inches(13,5)
for i,seq in enumerate(["flair","??","t1ce","??"]):
_indices = []
for patient_id in range(5,85,10):
img = io.imread(full_data[patient_id][seq], plugin='simpleitk')
depth_dimension = img.shape[0]
_slice = np.array([list(img[i,:,:].reshape(-1,1)) for i in range(depth_dimension)])
_slice = np.squeeze(_slice,axis = 2).T
slice_df = pd.DataFrame(_slice)
# correlation matrix
_df = slice_df.corr()
# indices or slice numbers whose correlation is nan
_indices.append([y for x in np.argwhere(_df.isnull().all(axis=1).values) for y in x])
plt.subplot(1,4,k)
plt.hist(x= list(chain.from_iterable(_indices)))
plt.title(seq.upper())
k+=1
从图 8-19 中,你可以清楚地观察到
图 8-19
跨切片的不变强度
- 总体趋势中从 0-5 和 140-154 的切片数显示所有序列的强度没有变化,因此这些切片可以忽略。
裁剪和填充
还有其他类型的预处理可以完成,例如将切片裁剪到一个较低的维度,然后将它们填充到一个标准维度。这通常是为了减小卷的大小。您将遵循另一种方法来减少图像体积上不必要的卷积。但是,你必须确保在这样做之后没有错位,如图 8-20 所示。
图 8-20
裁剪原始图像以减少尺寸
虽然您可以明显地看到大小已经减小而没有任何信息损失,但我认为最有效的方法是在创建体训练数据的 3d 块时处理它,方法是在块体上设置阈值,比如至少 10%或 5%的非零像素强度。在后面的章节中会有更多的介绍。
二维图像的图像分类
名为“处理 2-D 图像”的部分详细介绍了许多数据属性以及您可以用它们做什么。还记得关于像素间距以及如何对图像进行重新采样以均匀分布空间信息的讨论吗?这个预处理步骤帮助您有效地使用 CNN 进行处理,使得内核从一个图像单元中学习相同的信息。核是用于从图像中提取特征的过滤器(二维矩阵)。
图像预处理
直方图均衡
有时由于对比度差,X 射线图像需要增强以突出小的纹理和细节。这样做基本上是为了扩大图像像素值的范围。参见图 8-21
图 8-21
直方图均衡。来源:维基
现在,如果您的整个图像被限制在一个像素范围内,您可以简单地将当前的像素分布映射到一个更宽更均匀的分布,但是如果已经存在高强度和低强度的区域(也就是更大范围的像素值),您必须局部应用一种叫做自适应直方图均衡化的东西。
具体来说,您将使用 CLAHE 方法。它可以增强图像的局部对比度,增强图像各部分的边缘和曲线的可视性。
-
对比度限制:如果该区域的任何直方图超过对比度限制,它们将被剪切。
-
自适应直方图均衡化:将图像分成小块,称为瓦片。这些图块是直方图均衡化的。
您将使用 OpenCv 来直方图规范化您的图像。参见图 8-22 中的结果。
图 8-22
直方图均衡的结果
def histogram_equalization(img, clip_limit, grid_size):
"""
Histogram Equalization
"""
clahe = cv2.createCLAHE(clipLimit = clip_limit,
tileGridSize = grid_size)
img_clahe = clahe.apply(img)
return img_clahe
像素的各向同性均衡
为了确保均匀的像素间距,您必须对图像进行插值和重新采样。
def resample(img, x_pixel, y_pixel):
new_size = [1, 1]
size = np.array([x_pixel, y_pixel])
img_shape = np.array(img.shape)
new_shape = img_shape * size
new_shape = np.round(new_shape)
resize_factor = new_shape / img_shape
resampled_img = scipy.ndimage.interpolation.zoom(img, resize_factor)
return resampled_img
虽然你的像素现在间隔适当,但这导致了形状较低的图像,这意味着不同的像素间隔将导致不同的图像大小,所以现在你必须通过裁剪/填充/插值来重塑它们。我通常更喜欢插值而不是固定大小,因为原始形状和目标形状之间的差异并不大。同样,您将使用 Opencv。
模型创建
由于您在 DICOM 元数据中发现了一些重要信息,您将创建一个双输入单输出神经网络,其中一个分支接收一批二维图像,另一个分支接收缩放的特征列。参见图 8-23
图 8-23
医学图像分类模型
你们中的一些人可能不知道我接下来要用到的不同层次和术语。我推荐阅读杜默林等人分享的这本优秀指南,名为《深度学习卷积算法指南》。
我们先从 TensorFlow 库导入相关库开始。
from tensorflow.keras.models import Model, Sequential
from tensorflow.keras.layers import (
Input,
Conv2D,
Dropout,
MaxPooling2D,
concatenate,
BatchNormalization,
Flatten,
Dense
)
from tensorflow.keras.optimizers import Adam
METRICS = [
tf.keras.metrics.AUC(name='auc'),
]
接下来,创建一个函数作为卷积块(如图 8-23 中的虚线所示)。您提供了几个控件来创建此卷积块,即
-
一个具有任意数量滤镜的卷积层和一个特定
kernel_size.
的卷积 -
一个
BatchNormalization
层,用于批量标准化数据,从而减少数据之间的协变量偏移。 -
pooling
利用核中像素亮度的最大值来压缩信息,从而缩小图像的特征空间。 -
dropout
通过在训练时随机丢弃神经元来防止过拟合。
def convolution_block(input_layer, num_filters, kernel_size,
strides, padding = 'valid',
activation = 'selu',
batch_normalization = False,
pool_kernel = None, dropout_rate = None):
layer = Conv2D(num_filters, kernel_size, strides = strides,
padding=padding, activation=activation)(input_layer)
if batch_normalization:
layer = BatchNormalization()(layer)
if pool_kernel:
layer = MaxPooling2D(pool_kernel)(layer)
if dropout_rate:
layer = Dropout(dropout_rate)(layer)
return layer
现在您可以构建您的主函数来创建您想要的网络。
-
首先创建两个输入层,它告诉模型预期的输入大小。
-
然后,根据您想要训练的参数数量,您可以不断添加到卷积块中,并相应地选择内核和池。我通常更喜欢从大内核和无池开始。然后我会更深入地介绍这两者。
def build_model():
input_img = Input(TARGET_SHAPE+(1,))
input_feats = Input((6,))
cb1 = convolution_block(input_img, num_filters = 128, kernel_size = 8,
strides = 1, padding = 'valid',
batch_normalization = True,
activation = 'selu',
pool_kernel = None, dropout_rate = None)
cb2 = convolution_block(cb1, num_filters = 32, kernel_size = 8,
strides = 1, padding = 'valid',
activation = 'selu',
pool_kernel = 2, dropout_rate = None)
cb3 = convolution_block(cb2, num_filters = 8, kernel_size = 8,
strides = 1, padding = 'valid',
activation = 'selu',
pool_kernel = 2, dropout_rate = 0.2)
cb4 = convolution_block(cb3, num_filters = 4, kernel_size = 8,
strides = 1, padding = 'valid',
activation = 'selu',
pool_kernel = 2, dropout_rate = 0.2)
conv_flat = Flatten()(cb4)
cl1 = Dense(128, activation='selu')(conv_flat)
cl2 = Dense(64, activation='selu')(cl1)
cl3 = Dense(32, activation='selu')(cl2)
# Feature block
fl1 = Dense(4, activation='selu')(input_feats)
concat_layer = concatenate([cl3, fl1], axis = 1)
# prediction block
pl1 = Dense(16, activation = 'selu')(concat_layer)
pl2 = Dense(8, activation = 'selu')(pl1)
output = Dense(1, activation = 'sigmoid')(pl2)
return Model([input_img, input_feats], output)
准备输入数据
Keras 中的现成生成器不支持这样的多输入,因此您必须创建自己的定制生成器。
首先对视图位置和年龄箱进行一次性编码,这是您发现的与目标变量相关的两个最重要的特征变量。
bin_labels = ['0_20', '20_40', '40_60', '60_80', '80_plus']
train_labels_w_md['age_bucketed'] = pd.cut(train_labels_w_md['age'].astype(int),
bins = [0, 20, 40, 60, 80, max(train_labels_w_md['age'].astype(int))],
labels = bin_labels)
视图位置已经是一个分类变量,因此可以直接进行一键编码。
from sklearn.preprocessing import LabelBinarizer
age_binarizer = LabelBinarizer()
age_binarizer.fit(train_labels_w_md['age_bucketed'])
transformed_age = age_binarizer.transform(train_labels_w_md['age_bucketed'])
transformed_age_ohe = pd.DataFrame(transformed_age)
transformed_age_ohe.columns = ["age_bin_trans_"+str(i) for i in range(len(age_binarizer.classes_))]
view_pos_binarizer = LabelBinarizer()
view_pos_binarizer.fit(train_labels_w_md['view_pos'])
transformed_view_pos = view_pos_binarizer.transform(train_labels_w_md['view_pos'])
transformed_view_pos_ohe = pd.DataFrame(transformed_view_pos)
transformed_view_pos_ohe.columns = ["view_pos_trans"]
data = pd.concat([train_labels_w_md, transformed_age_ohe, transformed_view_pos_ohe], axis=1)
接下来,为图像数组定义一个预处理函数。除了直方图均衡化和各向同性均衡化,您还可以将图像转换为标准形状,并通过将每个像素除以 255(图像像素的最大值)来归一化像素值。
您还添加了另一个充当通道的维度。这样做是为了满足 Conv2D 层的要求。
def get_train_images(dicom_path, target_shape):
img = pydicom.read_file(dicom_path)
img_equalized = histogram_equalization(img.pixel_array, 4, (8,8))
img_isotropic = resample(img.pixel_array, img.PixelSpacing[1], img.PixelSpacing[0])
img_standardized = cv2.resize(img_isotropic, target_shape, interpolation = cv2.INTER_CUBIC)
# Pixel Standardization
img_standardized = np.array(img_standardized)/255
res = np.expand_dims(img_standardized, axis = 2)
return res
创建训练集和验证集。
from sklearn.model_selection import train_test_split
train, val = train_test_split(data,test_size=0.25, random_state=42)
TARGET_SHAPE = (224,224)
BATCH_SIZE = 32
最后,您创建您的生成器,类似于您在第四章中创建的生成器。你把你的多输入让给了网络。
def get_data_generator(df, target_shape, shuffle = True, batch_size=32):
"""
Generator function which yields the input data and output for different clusters
"""
img, feat_set, y = [], [], []
if shuffle:
df = df.sample(frac=1).reset_index(drop=True)
while True:
for i,row in df.iterrows():
feat_set.append(np.array(row[[x for x in df.columns if "_trans" in x]].tolist()))
img.append(get_train_images(os.path.join(DATA_DIR,"Train",row['patientId'] + ".dcm"), TARGET_SHAPE))
y.append(np.array([row['Target']]))
if len(feat_set) >= batch_size:
yield (np.array(img), np.array(feat_set)), y
img, feat_set, y = [], [], []
培养
在训练中,您分别为训练集和验证集调用generator
函数。请注意,并不总是建议以这种方式创建生成器,因为数据流水线并没有针对预取和创建批处理数据时的许多数据操作进行优化。如果有更多的计算和 RAM 来预处理数据并以期望的格式存储数据,就可以避免这种情况。
train_generator = get_data_generator(train, TARGET_SHAPE, True, BATCH_SIZE)
val_generator = get_data_generator(val, TARGET_SHAPE, True, BATCH_SIZE)
model = build_model()
model.compile(optimizer= 'adam',
loss = tf.keras.losses.BinaryCrossentropy(from_logits=True),
metrics=METRICS)
如图 8-24 所示,基于架构,训练参数的数量为 526,217,与我们对 ImageNet 等超大型图像模型的预期相差甚远。因此,根据您的计算资源,可以随意构建不同的架构,并试验性能和收敛速度。
图 8-24
模型摘要
history = model.fit(train_generator,
steps_per_epoch= len(train)//BATCH_SIZE,
epochs=10,
validation_data=val_generator,
validation_steps= len(val)//BATCH_SIZE)
三维图像的图像分割
当我讲述每种方法的主要挑战和主要发展/解决方案时,我已经深入讨论了各种图像分析方法。在本节中,您将重点关注三维图像的图像分割问题。
让我们快速回顾一下图像分割的关键点:
-
**什么事?**图像分割基于训练数据将给定图像分割成各种片段,也称为感兴趣区域。
-
生物医学细分的主要挑战:
-
捕获图像中的噪声会导致不均匀的强度。
-
目标器官或病变的大小和形状可能因患者而异。
-
病变/目标器官占据整个图像的非常小的区域的类别不平衡可以导致 ML 模型学习更多关于背景或局部最小值的知识。
-
图像预处理
根据您在上一节中对 BRATS 数据的分析以及为 MRI 图像推荐的其他常规预处理,您将执行以下预处理:
-
偏场校正
-
移除不需要的切片
创建用于训练的面片时,您将标准化中心像素强度并忽略空体积。
偏场校正
当捕获 MRI 图像时,偏置场可以通过减少诸如边缘的高频内容来模糊图像。它还影响像素的强度,使得相同的组织显示灰度变化。
对于肉眼来说,这种差异并不意味着什么,但对于 ML 算法来说,它可以产生巨大的差异。让我们纠正这种偏见。为此,您将使用 SimpleITK 库,它提供了 N4 字段校正。N4 偏置场校正算法是一种用于校正 MRI 图像数据中存在的低频强度不均匀性的流行方法,称为偏置场或增益场。更多详情请见 https://simpleitk.readthedocs.io/en/master/link_N4BiasFieldCorrection_docs.html
由于边缘/轮廓受偏置场的影响,您必须使用阈值算法来分离背景和前景像素。为此,您将使用 Otsu 的方法。注意 SimpleITK 库中还有其他自动阈值算法,如最大熵、三角形等。我敦促你尝试这些不同的变化。
Otsu 的方法计算量很大,因此可能需要很长时间才能完成,所以您将把这种偏差校正的结果保存在一个单独的文件夹中。
NEW_BASE_DIR = os.path.join(os.path.split(BASE_DIR)[0],
"PROCESSED_IMAGE")
你首先读取图像,然后使用 Otsu 的方法创建一个遮罩。这个遮罩只包含 1 和 0 来分隔前景和背景像素。之后,使用阈值遮罩对输入图像进行实地校正。
def correct_bias_field(input_path, output_path):
inputImage = sitk.ReadImage(input_path)
maskImage = sitk.OtsuThreshold( inputImage,
0, # Background Value
1, # Foreground Value
250 # Number of Histograms
)
# Casting to allow real pixel value
inputImage = sitk.Cast( inputImage, sitk.sitkFloat32 )
corrector = sitk.N4BiasFieldCorrectionImageFilter()
output = corrector.Execute( inputImage, maskImage)
# Since our original image followed the 16-bit pixel format
outputCasted = sitk.Cast(output,sitk.sitkVectorUInt16)
sitk.WriteImage(outputCasted,output_path)
您可以为每个患者和每个图像序列调用上述函数,并相应地保存结果。
processed_full_data = {}
for patient_id,v in full_data.items():
processed_full_data[patient_id] = {}
for seq,input_path in v.items():
print(f"Started Bias Correction for Patient {patient_id} and Sequence {seq.upper()}")
folder_name = os.path.split(os.path.split(input_path)[0])[-1]
file_name = os.path.split(input_path)[-1]
output_path = os.path.join(NEW_BASE_DIR,folder_name, file_name)
# Automatically create the directory that doesn't exist
if not os.path.exists(os.path.join(NEW_BASE_DIR,folder_name)):
os.makedirs(os.path.join(NEW_BASE_DIR,folder_name))
# Updating the new paths for 4 sequences
if seq == "label":
processed_full_data[patient_id].update({seq:input_path})
else:
processed_full_data[patient_id].update({seq:output_path})
correct_bias_field(input_path, output_path)
break
移除不需要的切片
这是创建训练数据之前预处理流水线的最后一步。将结果保存为 HDF5 文件格式,这样可以将单个数据拼接到一个文件中。
您将保存所有序列堆叠在一起的最终图像文件,以创建大小为4,135,240,240
的 4-D 卷和大小为135,240,240.
的标签卷。您将使用 h5py Python 包来保存 HDF5 文件。
import h5py
NEW_BASE_DIR = os.path.join(os.path.split(BASE_DIR)[0],
"PROCESSED_IMAGE","SLICE_CORRECTED")
# Automatically create the directory that doesn't exist
if not os.path.exists(NEW_BASE_DIR):
os.makedirs(NEW_BASE_DIR)
您运行连续的for
循环来遍历路径,并将 h5py 文件保存在上面创建的新目录中。
for patient_id,v in processed_full_data.items():
image_vol_w_seq = {}
image_mask = []
for seq,input_path in v.items():
image_volume = io.imread(input_path, plugin='simpleitk')
slices_to_keep = np.array([_slice for i,_slice in enumerate(image_volume) if i not in (list(range(5))+list(range(140,155)))])
if seq == "label":
# To enable one-hot encoding of these categories
# we make a continous range of classes from 0 to 3
slices_to_keep[slices_to_keep == 4] = 3
image_mask = np.copy(slices_to_keep)
else:
image_vol_w_seq[seq] = slices_to_keep
final_image = np.stack((image_vol_w_seq['flair'],
image_vol_w_seq['??'],
image_vol_w_seq['t1ce'],
image_vol_w_seq['??'])).astype('float')
# Check individual size of mask and train images
assert image_mask.shape == slices_to_keep.shape
assert final_image.shape == (4, ) + slices_to_keep.shape
# Initialize the HDF5 File
_path = os.path.join(NEW_BASE_DIR, f'{str(patient_id+1).zfill(3)}.h5')
_hf = h5py.File(_path, 'w')
# Use create_dataset to give dataset name and provide numpy array
_hf.create_dataset('X', data = final_image)
_hf.create_dataset('Y', data = image_mask)
# Close to write to the disk
_hf.close()
模型创建
3D U-Net 体系结构的灵感来自 U-Net 体系结构,该体系结构在成为医学图像分割的 SOTA 一段时间后积累了巨大的人气。该架构是由弗赖堡大学与谷歌的 Deepmind 团队合作在题为“3D U-Net:从稀疏注释中学习密集体积分割”的论文中介绍的。图 8-25 是同一篇论文中展示 U-Net 架构的图片。
图 8-25
3D U-Net 架构。蓝框代表要素地图。通道的数量在每个特征图上标出
U-Net 架构由一条**收缩路径(左侧)和一条扩张路径(右侧)**组成。
收缩路径遵循卷积网络的典型架构。
卷积层之后是非线性激活和汇集操作,以防止过拟合。有时添加 BatchNormalization 或其变体,以确保一批数据中的协变量变化不会突然影响梯度学习过程。协变量变化是观察到的批次间数据分布的变化。
在每个下采样步骤,特征通道被加倍,而扩展路径包括上采样和连接,随后是常规卷积运算。在此路径中,您尝试通过扩展特征尺寸来恢复压缩特征。从图 8-25 中绿色箭头指示的收缩路径中,以符合所需特征图形状的方式进行上采样。
U-Net 的主要优势在于,在进行上采样的同时,您还可以连接来自编码器/收缩网络的特征图。
让我们开始编码吧。首先,导入创建模型所需的相关层。
from tensorflow.keras.models import Model
from tensorflow.keras.layers import (
Input,
Activation,
Conv3D,
Conv3DTranspose,
MaxPooling3D,
UpSampling3D,
SpatialDropout3D,
concatenate,
BatchNormalization
)
from tensorflow.keras.optimizers import Adam
您创建了卷积块,它基本上创建了收缩和扩展路径。就像前几章一样,您继续使用 SELU 作为您的首选激活方式。
您还将使用批处理规范化层,它处理一批数据的协变量变化。您可以通过使用标志变量来控制它的使用。
需要注意的一点是data_format = 'channels_first'
的使用。这样做是为了告诉层,输入图像具有作为第一维的通道。不要迷茫;你的输入图像实际上是一个五维张量。
5+D tensor with shape: batch_shape + (channels, conv_dim1, conv_dim2, conv_dim3) if data_format='channels_first' or 5+D tensor with shape: batch_shape + (conv_dim1, conv_dim2, conv_dim3, channels) if data_format='channels_last'
def convolution_block(input_layer, n_filters, batch_normalization=False,
kernel=(3, 3, 3), activation='selu',
padding='same', strides=(1, 1, 1)):
"""
Creates Convolutional Block
"""
layer = Conv3D(n_filters, kernel, activation = 'selu', data_format = 'channels_first', padding = padding, strides = strides)(input_layer)
if batch_normalization:
layer = BatchNormalization(axis=1)(layer)
return layer
有时,为了防止过拟合,您可能希望在最大池化之后添加一个 dropout 层,但为了简单起见,我们不介绍它。对于那些想要的人,可以通过将 max-pooling 输出传递给 Dropout 层来添加SpatialDropout3D layer
。
以类似的方式,对于扩展路径,定义一个上卷积运算。要得到同样大小的图像,有各种方法。在下面的函数中,你可以看到其中的两个,分别是反卷积和上采样。
-
去卷积:使用滤波器、内核、填充和步长,就像卷积层一样,以获得所需大小的图像。
-
上采样:通过传递用于压缩图像的池大小,将图像调整到所需的大小。
def up_convolution(n_filters, pool_size, kernel_size = (2, 2, 2),
strides = (2, 2, 2),
deconvolution = False):
if deconvolution:
return Conv3DTranspose(filters=n_filters, data_format = 'channels_first',
kernel_size=kernel_size, strides=strides)
else:
return UpSampling3D(size=pool_size, data_format = 'channels_first')
创建如图 8-25 所示的架构,并返回模型进行训练。
def unet_model_3d(loss_function, input_shape=(4, 24, 160, 160),
pool_size = 2, n_labels = 3,
initial_learning_rate = 0.001,
deconvolution=False, depth = 4, n_base_filters = 32, metrics=[],
batch_normalization = True):
"""
U-Net 3D Model
"""
# Input Layer for the Image patch
inputs = Input(input_shape)
current_layer = inputs
levels = list()
# add levels with max pooling
for layer_depth in range(depth):
layer1 = convolution_block(input_layer = current_layer,
n_filters = \
n_base_filters * (2 ** layer_depth),
batch_normalization = \
batch_normalization)
layer2 = convolution_block(input_layer=layer1,
n_filters = \
n_base_filters * (2 ** layer_depth)* 2,
batch_normalization = \
batch_normalization)
# Do Max-Pooling until reaching the bridge
if layer_depth < depth - 1:
current_layer = MaxPooling3D(pool_size = pool_size, data_format = 'channels_first')(layer2)
levels.append([layer1, layer2, current_layer])
else:
current_layer = layer2
levels.append([layer1, layer2])
# add levels with up-convolution or up-sampling
for layer_depth in range(depth - 2, -1, -1):
up_convolution_layer = up_convolution(pool_size = pool_size,
deconvolution = deconvolution,
n_filters = \
current_layer.shape[1])(current_layer)
# Concatenate Higher and Lower Dimensions
concat = concatenate([up_convolution_layer, levels[layer_depth][1]], axis=1)
current_layer = convolution_block(
n_filters = levels[layer_depth][1].shape[1],
input_layer = concat, batch_normalization = batch_normalization)
current_layer = convolution_block(
n_filters=levels[layer_depth][1].shape[1],
input_layer=current_layer,
batch_normalization=batch_normalization)
final_convolution = Conv3D(n_labels, (1, 1, 1),
data_format = 'channels_first',
activation = 'sigmoid')(current_layer)
model = Model(inputs = inputs, outputs = final_convolution)
if not isinstance(metrics, list): metrics = [metrics]
model.compile(optimizer=Adam(lr = initial_learning_rate),
loss = loss_function,
metrics=metrics)
return model
准备输入数据
为了准备输入数据,您需要了解在给定上述模型的情况下,分段在理想情况下是如何工作的。
理想情况下,您希望获得全部信息,无论是通道/序列、深度维度等。,用于创建训练标签,但正如您可以从一个 HDF5 文件的文件大小中看到的,解压缩时图像大小将非常大。此外,这将导致卷积时更大的核大小,因此需要学习更多的参数。
所以你只能选择这两种方法之一:
-
通过一次输入堆叠图像中的每个切片和标签图像中的相应切片来训练分割模型,就像进行 2-D 卷积一样,但随后您会有意放弃沿深度维度呈现的空间信息,这些信息在确定肿瘤类型方面起着关键作用。一些肿瘤在轴位视图中看起来很小,但在矢状位视图中可能非常明显,因此这种技术不能提供良好的结果。
-
您还可以从体积立方体创建小的三维体积块,并捕获所有维度的空间信息。现在,这不是一个完美的技术,因为您可能会错过一些托管信息,但它允许您一次捕获更多的信息,因此它是首选的。
如图 8-26 所示,您将创建您的培训数据。根据每个图像卷需要的最大补丁数,重复此过程 n 次。
图 8-26
创建训练数据补丁的详细流程
首先创建标准化函数,该函数将像素亮度集中到平均值 0 和标准偏差 1。对于每个序列,你在深度上循环,在(240,240)
图像上居中缩放,然后将它们堆叠在同一个地方。
import tensorflow as tf
def standardize(image):
"""
Centers the image with mean of zero and sd = 1
"""
# initialize to array of zeros, with same shape as the image
standardized_image = np.zeros(image.shape)
# iterate over sequences
for c in range(image.shape[0]):
# iterate over the depth dimension
for z in range(image.shape[1]):
image_slice = image[c,z,:,:]
# subtract the mean from image_slice
centered = image_slice - np.mean(image_slice)
# divide by the standard deviation (only if it is different from zero)
if np.std(centered):
centered_scaled = centered / np.std(centered)
standardized_image[c, z, :, :] = centered_scaled
else:
standardized_image[c, z, :, :] = image_slice
return standardized_image
接下来,为训练数据创建所需的修补程序。为此,您需要为新图像创建一个新文件夹。
NEW_BASE_DIR = os.path.join(os.path.split(BASE_DIR)[0],
"FINAL_TRAIN_IMAGE")
# Automatically create the directory that doesn't exist
if not os.path.exists(NEW_BASE_DIR):
os.makedirs(NEW_BASE_DIR)
函数中发生了很多事情,但这些是正在发生的主要步骤。
-
首先,获取所需的补片大小、标签图像中的类别数量,以及从患者图像中获取所需补片数量的尝试次数。您需要尝试多次,因为您希望阈值至少有 4%的肿瘤区域。
-
通过为所有轴(x、y 和 z)选择一个随机起点,可以选择一个随机面片。
-
由于您选择了多个修补程序,因此修补程序的起点可能会相同。为此,您可以为已经选择的起点维护一个单独的列表,以避免任何重叠。
-
您可以比较轴元组或单个轴点。我选择了前者,但你可以尝试任何一个。
-
-
在标签图像上,您引入了一个新的维度,对肿瘤标签进行一次性编码。正如您在下面的代码示例中看到的,引入了一个新的维度。
您将使用tf.keras.utils.to_categorical
来完成这项工作。
Converts a class vector (integers) to binary class matrix.
Read more at https://www.tensorflow.org/api_docs/python/tf/keras/utils/to_categorical
tf.keras.utils.to_categorical([0,1,2,3], num_classes=4)
输出:-
array([[1., 0., 0., 0.],
[0., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.]], dtype=float32)
tf.keras.utils.to_categorical([[0,1,0,3],[0,0,2,3]], num_classes=4)
输出:
-
您删除了背景类,因为您对预测它不感兴趣。但是请记住,这仍然不会改变由于不平衡导致的稀疏性。
-
如果通过了背景比率,则标准化输入图像并将其与标签一起保存。
array([[[1., 0., 0., 0.],
[0., 1., 0., 0.],
[1., 0., 0., 0.],
[0., 0., 0., 1.]],
[[1., 0., 0., 0.],
[1., 0., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.]]], dtype=float32)
Note
对于那些想知道如何决定输出尺寸的人,请重温裁剪图像练习。你可以看到,在对头骨进行裁剪后,你剩下的尺寸是 164,123,因此你可以选择一个想要的尺寸(180, 160).
def get_multiple_patchs(image, label, patient_id,
save_dir,
out_dim = (180,160,24),
num_classes = 4,
max_tries = 1000,
num_patches = 5,
background_threshold=0.96):
"""
Extract random sub-volume from original images.
"""
num_channels, orig_z, orig_x, orig_y = image.shape
out_x, out_y, out_z = out_dim
all_patches = []
tries = 0
# try until you fail :P
prev_start = []
while (tries < max_tries) and (len(all_patches) < num_patches):
# Start from the corner randomly sample a voxel (volume box)
start_x = np.random.randint(0, orig_x - out_x + 1)
start_y = np.random.randint(0, orig_y - out_y + 1)
start_z = np.random.randint(0, orig_z - out_z + 1)
# Make sure you are choosing a unique starting point each time
while (start_x,start_y,start_z) in prev_start:
start_x = np.random.randint(0, orig_x - out_x + 1)
start_y = np.random.randint(0, orig_y - out_y + 1)
start_z = np.random.randint(0, orig_z - out_z + 1)
# extract relevant area of label
y = label[start_z: start_z + out_z,
start_x: start_x + out_x,
start_y: start_y + out_y]
# One-hot encode the tumor categories to add a 4-th dimension
y = tf.keras.utils.to_categorical(y,num_classes)
# compute the background ratio
bgrd_ratio = np.sum(y[:,:,:,0])/(out_x*out_y*out_z)
# increment tries counter
tries += 1
# check if background ratio is less than the maximum background
# threshold
if bgrd_ratio < background_threshold:
# make copy of the sub-volume and take all the channels/seq
X = np.copy(image[:,
start_z: start_z + out_z,
start_x: start_x + out_x,
start_y: start_y + out_y])
X_std = standardize(X)
# we will also make sure that we bring the num class dimension
# as the first axis
y = np.moveaxis(y, 3, 0)
# Exclude the background class as we don't want to predict it
y = y[1:, :, :, :]
all_patches.append([X_std, y])
# Initialize the HDF5 File
_path = os.path.join(save_dir, f'{str(patient_id).zfill(3) + "_" + str(len(all_patches))}.h5')
_hf = h5py.File(_path, 'w')
# Use create_dataset to give dataset name and provide numpy array
_hf.create_dataset('X', data = X_std)
_hf.create_dataset('Y', data = y)
# Close to write to the disk
_hf.close()
return all_patches
最后,将所有补丁保存到培训目录中。
processed_path = glob.glob(os.path.join(os.path.join(os.path.split(BASE_DIR)[0],
"PROCESSED_IMAGE","SLICE_CORRECTED"),"*.h5"))
for _path in processed_path:
with h5py.File(_path, 'r') as f:
_image = f.get("X")
_label = f.get("Y")
_patient_id = int(os.path.split(_path)[-1].replace(".h5",""))
x = get_multiple_patchs(_image, _label, _patient_id, NEW_BASE_DIR,
out_dim = (180, 160,24),
num_classes = 4,
max_tries = 1000,
num_patches = 5,
background_threshold=0.96)
您现在已经准备好训练您的模型了。
培养
由于您正在处理大量的大尺寸图像,通常不建议一次加载所有的图像。您将使用发电机加载它们。
并非 TensorFlow-Keras 中的所有现成生成器都需要图像文件。有一些出色的生成器,但遗憾的是,它们不适用于 HDF5 数据,因此您必须自己编写一个。
Note
感兴趣的,可以查看本次回购的flow_from_* functions
:https://keras.io/api/preprocessing/image/
.
我正在利用这个优秀教程中的代码。跟随它来拓宽你的理解: https://stanford.edu/~shervine/blog/keras-how-to-generate-data-on-the-fly
.
发电机的代码在本书的 GitHub repo 中共享。一定要去看看。
为了训练你的模型,最后一步是确定一个损失函数。最有可能的是,在处理多个类时,您会倾向于选择交叉熵损失,但交叉熵损失在处理高度不平衡的数据集时并不那么有效。
我们来了解一下。图 8-27 显示了前景像素用 1 表示,背景像素用 0 表示的图像块。
图 8-27
实际和预测的图像像素
您可以看到,随着非重叠区域的增加,二进制交叉熵几乎呈线性增加,而 dice 系数即使没有重叠也不会达到 0,也没有线性减少,因此它可以更好地处理不平衡。
邹等人的题为“基于空间重叠指数的图像分割质量的统计验证”的论文详细讨论了 dice 系数和统计验证。我强烈建议你看一看。
您可以根据您的用例进一步调整损失函数。因为您的 DL 模型实际上输出的不仅仅是概率,所以您可以修改实际的 Dice 函数来处理概率而不是实际的二进制值。
该书的 GitHub repo 中分享了 Soft Dice 函数的代码,可以查看一下。另外,请注意,在阅读一些论文时,您不会发现任何关于 epsilon 的内容,但在实施过程中您会发现。不要害怕;在现实世界的实现中,使用拉普拉斯平滑来避免除法误差是一种常见的做法。
最终确定损失函数后,您现在需要一个性能指标来评估训练好的模型。为此,您将使用 dice 系数,这是评估细分模型性能的标准指标。
from tensorflow.keras import backend as K
K.set_image_data_format("channels_first")
def dice_coefficient(y_true, y_pred, axis=(1, 2, 3),
laplace_smoothing_factor=0.00001):
dice_numerator =2 *K.sum(y_pred*y_true,axis) + laplace_smoothing_factor
dice_denominator = K.sum(y_pred,axis) + K.sum(y_true,axis) + laplace_smoothing_factor
# For multiple classes take the mean across each axis
dice_coefficient = K.mean(dice_numerator/dice_denominator)
return dice_coefficient
def soft_dice_loss(y_true, y_pred, axis=(1, 2, 3),
laplace_smoothing_factor=0.00001):
"""
Compute mean soft dice loss over all Multiple classes.
"""
dice_numerator =2 *K.sum(y_pred*y_true,axis) + laplace_smoothing_factor
dice_denominator = K.sum(y_pred**2,axis) + K.sum(y_true**2,axis) + laplace_smoothing_factor
dice_loss = 1 - K.mean(dice_numerator / dice_denominator)
return dice_loss
详见 https://en.wikipedia.org/wiki/S%C3%B8rensen%E2%80%93Dice_coefficient
。
model = unet_model_3d(depth = 3,
pool_size= 2,
input_shape=(4,180, 160,24),
n_base_filters = 32,
loss_function=soft_dice_loss, metrics=[dice_coefficient])
import h5py
NEW_BASE_DIR = os.path.join(os.path.split(BASE_DIR)[0],
"FINAL_TRAIN_IMAGE")
all_patches = glob.glob(os.path.join(NEW_BASE_DIR,"*.h5"))
from sklearn.model_selection import train_test_split
train_data, val_data = train_test_split(all_patches, test_size=0.33, random_state=42)
BATCH_SIZE = 5
# Get generators for training and validation sets
train_generator = BatchDataGenerator(train_data, batch_size = BATCH_SIZE, dim = (180, 160, 24))
valid_generator = BatchDataGenerator(val_data, batch_size = BATCH_SIZE, dim = (180, 160, 24))
为了训练一个生成器,你要传递一个叫做steps per epoch.
的东西,这个步骤数就是一批中要训练的样本数。如果有 100 个训练样本,您想要 5 个批次,那么steps per epoch = ceil(num_samples/batch_size).
steps_per_epoch = len(train_data)//BATCH_SIZE
n_epochs=10
validation_steps = len(val_data)//BATCH_SIZE
history = model.fit(train_generator,
steps_per_epoch=steps_per_epoch,
epochs=n_epochs,
use_multiprocessing=False,
validation_data=valid_generator,
validation_steps=validation_steps)
Epoch 1/10
120/120 [==============================] - 860s 7s/step - loss: 0.3559 - dice_coefficient: 0.4725 - val_loss: 0.4446 - val_dice_coefficient: 0.4008
Epoch 2/10
120/120 [==============================] - 802s 7s/step - loss: 0.3350 - dice_coefficient: 0.5026 - val_loss: 0.3232 - val_dice_coefficient: 0.5164
Epoch 3/10
120/120 [==============================] - 4914s 41s/step - loss: 0.3214 - dice_coefficient: 0.5237 - val_loss: 0.5040 - val_dice_coefficient: 0.3589
Epoch 4/10
120/120 [==============================] - 505s 4s/step - loss: 0.3194 - dice_coefficient: 0.5305 - val_loss: 0.3561 - val_dice_coefficient: 0.4605
Epoch 5/10
120/120 [==============================] - 517s 4s/step - loss: 0.3023 - dice_coefficient: 0.5525 - val_loss: 0.4168 - val_dice_coefficient: 0.4328
Epoch 6/10
120/120 [==============================] - 518s 4s/step - loss: 0.3066 - dice_coefficient: 0.5552 - val_loss: 0.3478 - val_dice_coefficient: 0.5206
Epoch 7/10
120/120 [==============================] - 521s 4s/step - loss: 0.3054 - dice_coefficient: 0.5596 - val_loss: 0.3900 - val_dice_coefficient: 0.4615
Epoch 8/10
120/120 [==============================] - 238s 2s/step - loss: 0.2914 - dice_coefficient: 0.5746 - val_loss: 0.3515 - val_dice_coefficient: 0.5378
Epoch 9/10
120/120 [==============================] - 294s 2s/step - loss: 0.3022 - dice_coefficient: 0.5697 - val_loss: 0.3471 - val_dice_coefficient: 0.5286
Epoch 10/10
120/120 [==============================] - 374s 3s/step - loss: 0.2864 - dice_coefficient: 0.5864 - val_loss: 0.3279 - val_dice_coefficient: 0.5300
性能赋值
参见图 8-28 获取您的模型的性能和损耗图。注意,该模型表现良好,dice 系数随着每个时期而提高。有多种方法可以通过使用更大的过滤器尺寸和漏失来进一步改进模型,以防止过拟合和更大的深度。
图 8-28
三维 U-Net 模型的性能和损耗图
医学图像的迁移学习
迁移学习是一种从大型预训练网络中释放知识的方法,该网络在巨大的带注释的数据集上进行训练,以解决领域或其训练目的之外的任务。现在几乎所有的计算机视觉问题都使用迁移学习,实现 SOTA。
但与自然图像不同,自然图像始终由三个通道组成,不同集合之间的差异较小,医学图像可以有根本的不同。它们不仅可以具有变化的通道长度,而且这些图像的像素强度是基于医疗设备及其应用的物理特性来决定的。
因此,让我们了解迁移学习是否适用于医学图像,如果不适用,我们可以做些什么来实现它。
这一部分主要从谷歌大脑和康奈尔大学在 NIPS 2019 年发表的题为“输血:理解医学成像的迁移学习”的论文中提取思想。
作者试图通过以下方式来理解它:
-
性能评估:将从随机初始化训练并直接应用于任务的模型与在 ImageNet 上为相同任务预先训练的模型进行比较。
-
表示分析:使用典型相关分析比较不同模型的隐藏表示。
-
对收敛的影响:由于特征被重用,模型收敛所需的时间显著减少。
在这三个实验的基础上,作者得出了如下结论:
-
迁移学习不会显著影响医学成像任务的表现。这意味着预训练模型是过度参数化的。
-
对于上下文,ImageNet 有 1000 个类别,而医学图像有小得多的预测向量。
-
输入图像的大小差别很大。作为 ImageNet 一部分的自然图像具有 224 X 224 的尺寸,而医学图像可以具有更大的尺寸,例如 512 X 512、1024 X 1024 等。这种图像太大,不能直接输入神经网络。
-
-
发现使用来自网络的最后两层的预先训练的权重对收敛具有最大的影响。
我认为医学图像的迁移学习不如自然图像有效。对于医学图像,我们仍然应该坚持基于图像的形态和手头的任务来训练我们自己的网络的方法。
结论
祝贺你读完了这么长的一章。本章涵盖了与医学图像和人工智能应用相关的广泛信息。您探索了用于捕捉医学图像的不同模态,以及某些模态如何对某些解剖结构有用。然后,我们讨论了存储这些图像的两种不同格式 DICOM 和 NIFTI,以及如何利用与图像相关的元数据来拓宽您对所接收图像的理解。通常,这被更令人着迷的深度学习应用程序忽略了,但元数据包含了很多见解。最后,您了解了二维和三维图像以及可用于解决分类和分割问题的不同架构。
你在本章中学到的东西有巨大的现实价值。分类和定位可用于疾病筛查,也可用于紧急诊断和偶然发现。最近,随着 COVID 病例的大量出现,导致医护人员意外短缺,像分割这样的技术被用于诊断病情危急的患者,因此他们首先得到护理。分割还有助于识别肿瘤轮廓和肿瘤学放射治疗的敏感区域。
最后,深度学习和医疗人工智能领域不仅限于分类或细分。还有许多其他领域有应用高级人工智能技术的巨大潜力,例如从不同序列捕获的图像的配准,从设备到医生可用的图像的重建,以及具有临床背景的图像检索,可以帮助医生跟踪过去对一个病例做了什么。
参考
-
B.H. Menze,A. Jakab,S. Bauer,J. Kalpathy-Cramer,K. Farahani,J. Kirby,等人,“多模态脑肿瘤图像分割基准(BRATS),”IEEE 医学成像汇刊 34(10),1993-2024(2015)DOI:10.1109/TMI . 20143676
-
南 Bakas,H. Akbari,A. Sotiras,M. Bilello,M. Rozycki,J.S. Kirby 等人,“利用专家分割标签和放射特征推进癌症基因组图谱神经胶质瘤 MRI 集合”,自然科学数据,4:170117(2017)DOI:10.1038/sdata . 2017.117。
-
南 Bakas,M. Reyes,A. Jakab,S. Bauer,M. Rempfler,a .克里米等人,“在 BRATS 挑战中确定用于脑肿瘤分割、进展评估和总体生存预测的最佳机器学习算法”,arXiv 预印本 arXiv:1811.02629 (2018)。
-
南 Bakas,H. Akbari,A. Sotiras,M. Bilello,M. Rozycki,J. Kirby 等人,“TCGA-GBM 集合术前扫描的分割标签和放射特征”,癌症成像档案,2017 年。DOI:10.7937/K9/tcia . 2017 . klxwjj 1 q。
-
南 Bakas,H. Akbari,A. Sotiras,M. Bilello,M. Rozycki,J. Kirby 等人,“TCGA-LGG 集合术前扫描的分割标签和放射特征”,癌症成像档案,2017 年。DOI:10.7937/K9/tcia . 2017 . gjq 7 r0ef。