本章涵盖以下内容：
- 加载并处理原始数据文件
- 实现一个用于表示数据的 Python 类
- 将数据转换为 PyTorch 可用的格式
- 对训练数据和验证数据进行可视化
现在我们已经讨论了项目的高层目标，也勾勒了数据在系统中的流动方式，接下来就进入本章要做的具体内容。现在是时候为原始数据实现基础的数据加载和数据处理流程了。这里介绍的技术是基础性的，适用于你今后承担的任何大型项目。对于那些少见的、能够事先拿到已经准备完善的数据的研究者来说：你们真幸运！而我们其余大多数人，则还得忙着写加载和解析数据的代码。图 12.1 展示了第 11 章中给出的项目高层结构图。本章剩余部分将聚焦于第 1 步：数据加载。
图 12.1 我们的端到端肺癌检测项目，本章聚焦于其中的主题：第 1 步，数据加载
我们的目标是：给定原始 CT 扫描数据和一份这些 CT 的标注列表，能够产出一个训练样本。听起来似乎很简单，但在我们真正能够加载、处理并提取出所关心的数据之前，还需要发生很多事情。图 12.2 展示了我们需要完成哪些步骤，才能把原始数据变成训练样本。幸运的是，我们上一章已经提前对数据有了初步认识，不过在这一方面我们还有更多工作要做。
图 12.2 将原始数据转换为 sample tuple 所需的数据变换流程。这些 sample tuple 将作为模型训练流程的输入。
这是一个关键时刻：我们开始把沉重而粗糙的原始数据进行转化。即便还不能说是点石成金，至少也要先把它变成神经网络最终能“炼成金子”的材料。我们最早在第 4 章中讨论过这种把原始数据转成张量的机制。
12.1 原始 CT 数据文件
我们的原始 CT 数据构成了 LUNA 数据集的大部分内容，存放在各个 subset 目录中。每个 subset 目录包含若干个 CT 扫描，而每个扫描由两个文件表示：一个 .mhd 文件和一个 .raw 文件。.mhd 文件包含元数据头信息，.raw 文件则包含构成三维数组的原始字节。每个文件名都以一个唯一标识符开头，这个标识符叫作该 CT 扫描的 series UID（来自医学数字成像与通信标准 DICOM 的命名体系）。例如，对于 series UID 为 1.2.3 的扫描，会对应两个文件：1.2.3.mhd 和 1.2.3.raw。
我们的 Ct 类会读取这两个文件，并生成三维数组张量，以及一个转换矩阵，用于把病人坐标系（我们会在 12.4 节中更详细讨论）转换为三维数组张量所需的 index、row、column 坐标；在图中，这些坐标显示为 (I, R, C)，在代码中则以 _irc 作为变量后缀。现在你不必过度纠结这些细节；只要记住，在把这些坐标真正应用到 CT 数据之前，我们需要做一些坐标系转换。
我们还会加载 LUNA 提供的标注数据。这些数据会给出一组结节坐标，每个结节还带有一个恶性标志，以及相关 CT 扫描的 series UID。通过把结节坐标和坐标系变换信息结合起来，我们就能得到该结节中心所在体素的 index、row 和 column。
有了 (I, R, C) 坐标之后，我们就可以从 CT 数据中裁剪出一个小的三维块，作为模型输入。我们还会构造训练样本 tuple 的其余部分，最终这个 tuple 将包含：样本数组、结节状态标志、series UID，以及这个样本在该 CT 结节候选列表中的索引。这个 sample tuple 正是 PyTorch 对我们的 Dataset 子类所期待的返回内容，也代表着我们从最初的原始数据通向 PyTorch 张量标准结构的桥梁的最后一段。
对数据进行限制或裁剪、以免模型淹没在噪声里，是非常重要的；与此同时，也要确保我们裁剪得不要太激进，以至于把真正的信号也切掉了。我们还希望确保数据范围行为良好，尤其是在归一化之后。通过 clamp（截断）去除异常值会很有帮助，特别是在数据中容易出现极端离群值时。
12.2 解析 LUNA 的标注数据
我们需要做的第一件事，就是开始加载数据。在一个新项目中，这通常是个不错的起点。无论如何，首先确保自己知道如何处理原始输入数据，这是必须要做的；而了解数据在加载之后会呈现什么样子，也能帮助我们决定早期实验的结构。我们当然可以尝试直接加载单个 CT 扫描，但我们认为，先去解析 LUNA 提供的 CSV 文件更有意义，因为这些文件包含了每个 CT 扫描中感兴趣点的信息。正如图 12.3 所示，我们希望从中获取一些坐标信息、一个表明该坐标是否是结节的标记，以及 CT 扫描的唯一标识符。由于 CSV 文件里的信息种类更少，而且更容易解析，我们希望在真正开始加载 CT 之前，它们能给我们一些线索，帮助我们知道之后该关注什么。
图 12.3 candidates.csv 中的 LUNA 标注包含了 CT 的 series、结节候选的位置，以及一个指示该候选是否真的是结节的标志。
candidates.csv 文件包含了所有那些“看起来可能像结节”的肿块信息，不管这些肿块最终是恶性肿瘤、良性肿瘤，还是完全别的什么东西。我们将以此为基础建立一个完整的候选列表，再把它划分为训练集和验证集。下面这段 Bash shell 会话展示了该文件包含的内容：
$ wc -l candidates.csv #1
551066 candidates.csv
$ head data/part2/luna/candidates.csv #2
seriesuid,coordX,coordY,coordZ,class #3
1.3...6860,-56.08,-67.85,-311.92,0
1.3...6860,53.21,-244.41,-245.17,0
1.3...6860,103.66,-121.8,-286.62,0
1.3...6860,-33.66,-72.75,-308.41,0
...
$ grep ',1$' candidates.csv | wc -l #4
1351
#1 统计文件行数
#2 打印文件前几行
#3 .csv 文件第一行定义列名
#4 统计以 1 结尾的行数，也就是表示该候选为结节的记录数
注意：seriesuid 列中的值为了适配纸面排版被截断了。
因此，我们有 551,000 行数据，每一行都包含一个 seriesuid（代码中我们称作 series_uid）、一组三维空间中的 (X, Y, Z) 坐标，用来描述候选结节的位置，以及一个表示结节状态的 class 列。class 列是布尔值：0 表示该候选不是结节（例如血管），1 表示该候选是结节，而这个结节可能是恶性的，也可能是良性的。被标记为结节的候选共有 1,351 个。
annotations.csv 文件则包含了一部分被标记为结节的候选信息。我们尤其关心其中的 diameter_mm 信息：
$ wc -l annotations.csv
1187 annotations.csv #1
$ head data/part2/luna/annotations.csv
seriesuid,coordX,coordY,coordZ,diameter_mm #2
1.3.6...6860,-128.6994211,-175.3192718,-298.3875064,5.651470635
1.3.6...6860,103.7836509,-211.9251487,-227.12125,4.224708481
1.3.6...5208,69.63901724,-140.9445859,876.3744957,5.786347814
1.3.6...0405,-24.0138242,192.1024053,-391.0812764,8.143261683
...
#1 这个数字与 candidates.csv 中的不一样。
#2 最后一列也不一样。
我们拥有大约 1,200 个结节的大小信息。这很有用，因为我们可以利用它，确保训练集和验证集中都包含具有代表性的结节大小分布。否则，就有可能验证集里恰好只剩下一些极端值，从而让人误以为模型表现不佳。这个原则具有广泛适用性：无论你是在预测一部电影会有多少播放量、把文本分类为垃圾信息与否，还是在长音频中检测单词，纳入多样化的数据都能避免结果偏斜，并确保模型在真实世界场景中表现稳定。拥有一个具有代表性的数据集，是构建可靠机器学习模型的关键。
12.2.1 训练集与验证集
对于任何标准的监督学习任务（分类是最典型的例子），我们都会将数据划分为训练集和验证集。我们希望这两个集合都能够代表我们预期在真实世界中看到并正常处理的输入范围。如果其中任意一个集合与实际使用场景存在实质性差异，那么模型的行为很可能会与我们的预期不同——所有训练过程和统计数据，一旦迁移到生产环境中，就不再具有预测性！我们并不是要把这件事变成一门精确科学，但在今后的项目中，你需要始终留意那些迹象：它们可能表明你训练和测试所用的数据，并不适合你的实际运行环境。
再回到我们的结节。我们将按照大小对它们排序，然后每隔 N 个取一个放入验证集。这样应该能给我们带来想要的那种代表性分布。不幸的是，annotations.csv 中提供的位置信息，并不总是能与 candidates.csv 中的坐标精确对齐：
$ grep 100225287222365663678666836860 annotations.csv
1.3.6...6860,-128.6994211,-175.3192718,-298.3875064,5.651470635 #1
1.3.6...6860,103.7836509,-211.9251487,-227.12125,4.224708481
$ grep '100225287222365663678666836860.*,1$' candidates.csv
1.3.6...6860,104.16480444,-211.685591018,-227.011363746,1
1.3.6...6860,-128.94,-175.04,-297.87,1
#1 这两个坐标彼此非常接近。
如果我们把两个文件中对应的坐标截断一下来看，就得到 (-128.70, -175.32, -298.39) 与 (-128.94, -175.04, -297.87)。由于这个结节的直径大约为 5 mm，很显然，这两个点都被用来表示该结节的“中心”，但它们并不完全重合。面对这种数据不一致，一种完全合理的选择是：认为处理这种不匹配不值得，于是直接忽略这个文件。不过，我们还是决定做这部分繁琐工作，把它们对齐起来，因为真实世界的数据集往往就是这么不完美，而这正是你在把多个不同来源的数据组合起来时，必须要做的那种典型工作。
12.2.2 统一标注数据与候选数据
现在我们已经知道原始数据文件长什么样了，接下来就来构建一个 getCandidateInfoList 函数，把这些数据拼接整合起来。我们会在文件开头定义一个 named tuple，用来保存每个结节的信息。
代码清单 12.1 获取候选数据（dsets.py）
from collections import namedtuple
# ... line 27
CandidateInfoTuple = namedtuple(
'CandidateInfoTuple',
'isNodule_bool, diameter_mm, series_uid, center_xyz',
)
这些 tuple 还不是我们的训练样本，因为它们缺少我们所需要的 CT 数据块。相反，它们表示的是一个经过清洗、整理和统一之后的人类标注数据接口。把“处理脏数据”这件事和“模型训练”明确隔离开来非常重要。否则，你的训练循环会很快变得杂乱无章，因为其中不断夹杂着各种特例处理和其他干扰，而训练代码本该专注于训练本身。
提示：要把负责数据清洗的代码与项目其他部分清晰分离。如果有必要，不要害怕把数据重写一遍并保存回磁盘。
我们的候选信息列表将包含：结节状态（也就是模型要学习分类的目标）、直径（有助于在训练时覆盖不同大小的结节，因为大结节和小结节不会具有同样的特征）、series UID（用于定位正确的 CT 扫描），以及候选中心点（用于在更大的 CT 中找到该候选位置）。构建这些 NoduleInfoTuple 实例列表的函数，一开始先使用一个内存缓存装饰器，然后获取磁盘上实际存在的文件列表：
@functools.lru_cache(1) #1
def getCandidateInfoList(require_on_disk=True): #2
data_directory = Path(f"{mhd_data_folder}")
mhd_files = list(data_directory.rglob("subset*/*.mhd"))
if not mhd_files:
print(f"Warning: No .mhd files found under {data_directory}")
present_on_disk_set = {path.stem for path in mhd_files}
#1 标准库提供的内存缓存
#2 require_on_disk 默认会过滤掉那些所属数据子集尚未下载到本机的 series
由于解析某些数据文件可能比较慢，我们会把这个函数调用的结果缓存到内存中。这样做之后会非常有用，因为在后续章节中我们会更频繁地调用这个函数。通过谨慎地在内存或磁盘上应用缓存来加速数据流水线，往往能显著提高训练速度。你在自己的项目中也要留意这类机会。
在获取候选 series UID 之后，我们希望把 annotations.csv 中的直径信息合并进来。首先，我们需要按 series_uid 对标注进行分组，因为这是我们用来交叉引用两个文件中各行数据的第一个关键字段：
diameter_dict = {}
with open('data/part2/luna/annotations.csv', "r") as f:
for row in list(csv.reader(f))[1:]:
series_uid = row[0]
annotationCenter_xyz = tuple([float(x) for x in row[1:4]])
annotationDiameter_mm = float(row[4])
diameter_dict.setdefault(series_uid, []).append(
(annotationCenter_xyz, annotationDiameter_mm)
)
接着，我们再利用 candidates.csv 中的信息构建完整候选列表：
candidateInfo_list = []
with open('data/part2/luna/candidates.csv', "r") as f:
for row in list(csv.reader(f))[1:]:
series_uid = row[0]
if series_uid not in present_on_disk_set and require_on_disk: #1
continue
isNodule_bool = bool(int(row[4]))
candidateCenter_xyz = tuple([float(x) for x in row[1:4]])
candidateDiameter_mm = 0.0
for annotation_tup in diameter_dict.get(series_uid, []):
annotationCenter_xyz, annotationDiameter_mm = annotation_tup
for i in range(3):
delta_mm = abs(candidateCenter_xyz[i] - annotationCenter_xyz[i])
if delta_mm > annotationDiameter_mm / 4: #2
break
else:
candidateDiameter_mm = annotationDiameter_mm
break
candidateInfo_list.append(CandidateInfoTuple(
isNodule_bool,
candidateDiameter_mm,
series_uid,
candidateCenter_xyz,
))
#1 如果某个 series_uid 不存在于磁盘上，说明它属于一个我们没有下载的数据子集，应当跳过。
#2 先用直径除以 2 得到半径，再把半径再除以 2，以确保两个结节中心点相对结节尺寸而言不会相差太远。（这里做的是边界框式判断，而不是真正的欧氏距离判断。）
对于某个给定 series_uid 下的每一条候选记录，我们都会遍历之前收集到的、属于同一 series_uid 的标注，检查这两个坐标是否足够接近，以至于可以认为它们描述的是同一个结节。如果是，那很好——现在我们就有了该结节的直径信息。如果没有找到匹配项，也没关系；我们就把这个结节的直径暂时设为占位值 0.0。由于我们目前使用这些直径信息，仅仅是为了让训练集和验证集中的结节大小分布更合理，因此某些结节直径是占位值，按理说问题不大。不过我们也要记住自己做了这个假设，以防之后发现它不成立。
为了仅仅把结节直径合并进来，就要写这么多略显繁琐的代码。很遗憾，根据你的原始数据具体情况，这种数据操作与模糊匹配的工作其实相当常见。不过，一旦做到这一步，接下来我们就只需要对数据排序并返回即可：
candidateInfo_list.sort(reverse=True) #1
return candidateInfo_list
#1 所有真实结节样本会排在前面，并且从最大的开始；之后才是所有非结节样本（它们没有结节尺寸信息）。
noduleInfo_list 中 tuple 成员的顺序，正是由这个排序决定的。我们使用这种排序方法，是为了帮助确保在对数据做切片时，切出来的那一段能包含一部分具有代表性的真实结节，并且结节直径的分布也比较合理。
12.3 加载单个 CT 扫描
接下来，我们需要能够把磁盘上的 CT 数据——那堆比特——转换成一个 Python 对象，从中可以提取出三维结节密度数据。图 12.4 展示了这条从 .mhd 和 .raw 文件到 Ct 对象的路径。我们的结节标注信息就像一张地图，指向原始数据中那些有意思的部分。而在真正沿着这张地图去找到目标数据之前，我们得先把原始数据变成一种可寻址的形式。
图 12.4 加载一份 CT 扫描会产生一个体素数组，以及一个从病人坐标到数组索引的转换关系。这个过程对于把原始数据转换为结构化、机器可读、适合分析的格式至关重要。.mhd 和 .raw 文件提供原始体数据，而这些数据会被转换为便于索引与处理的三维体素数组。
CT 扫描的原生文件格式是 DICOM（<www.dicomstandard.org>）。DICOM 标准的第一个版本写于 1984 年。正如我们对那个时代的一切计算相关标准大概都会有所预期的那样，它多少有点混乱（例如，某些现已废弃的章节曾经专门讨论应该使用什么数据链路层协议，因为那时候以太网还没有取得最终胜利）。
注意：我们已经替你完成了为这些原始数据文件寻找合适解析库的工作；但对其他那些你可能从未听说过的格式，你还是得自己去找解析器。我们强烈建议你花时间去做这件事！Python 生态里几乎什么格式都有解析器，你的时间几乎肯定更应该花在项目真正有创新性的部分，而不是自己去给某种冷门数据格式写解析器。
好在，LUNA 已经把本章要用的数据转换成了 MetaIO 格式，这比原始 DICOM 好用得多。即便你从未听说过这种格式，也不必担心！我们可以把数据文件格式当成一个黑盒，利用开源的 SimpleITK 库把它们加载成更熟悉的 NumPy 数组：
class Ct:
def __init__(self, series_uid):
import SimpleITK as sitk
mhd_path = glob.glob(
f'{mhd_data_folder}/{series_uid}.mhd' #1
)[0]
ct_mhd = sitk.ReadImage(mhd_path) #2
ct_a = np.array(sitk.GetArrayFromImage(ct_mhd), dtype=np.float32) #3
#1 我们并不关心某个 series_uid 具体属于哪个 subset，因此这里对 subset 使用通配。
#2 sitk.ReadImage 在读取传入的 .mhd 文件时，也会隐式读取对应的 .raw 文件。
#3 重新创建一个 np.array，因为我们希望把数值类型转换为 np.float32。
在真实项目中，你当然应该理解原始数据里包含了哪些信息；但依赖像 SimpleITK 这样的第三方代码去解析磁盘上的比特，完全没有问题。至于到底应该在多大程度上了解输入细节、又在多大程度上可以放心接受数据加载库交给你的结果，这种平衡，大概需要一些经验才能把握。只要记住，我们真正关心的主要是数据，而不是比特。重要的是信息本身，而不是它在磁盘上如何表示。
能够唯一标识某个具体样本，通常很有用。比如，当你需要明确指出到底是哪一个样本导致了问题，或者哪个样本的分类结果很差时，这种唯一标识会极大提升我们隔离与调试问题的能力。具体到不同数据类型，这种唯一标识有时是一个数字或字符串，有时则更复杂，比如一个 tuple。
我们用 CT 扫描创建时分配的 series instance UID（series_uid）来标识具体的 CT 扫描。DICOM 大量使用唯一标识符（UID）来标识单个 DICOM 文件、文件组、治疗过程等等。这些标识符在概念上类似于 UUID（docs.python.org/3/library/u…），但它们的生成过程不同，格式也不同。对我们而言，只需把它们看作不透明的 ASCII 字符串，用作引用不同 CT 扫描的唯一键即可。按照官方定义，DICOM UID 中合法的字符只有数字 0 到 9，以及点号 .。
前面提到的 10 个子集，每个大约包含 90 个 CT 扫描（总共 888 个），而每个 CT 扫描都由两个文件表示：一个扩展名为 .mhd，一个扩展名为 .raw。不过，这种“数据被拆分成多个文件”的实现细节，被 sitk 的这些函数很好地隐藏起来了，我们不需要直接操心。
到这里为止，ct_a 已经是一个三维数组了。它的三个维度全都是空间维度。由于只有单一强度通道，类似灰度图像，所以通道维并不会显式体现在数组中。
12.3.1 Hounsfield 单位
回想一下我们前面说过的话：我们需要理解的是数据，而不是存储数据的比特。这里正好是一个绝佳例子。如果不理解数据数值及其范围的细节，我们就可能把一些会妨碍模型学习目标的值喂给模型。
继续在 Ct 类的 __init__ 方法里，我们需要对 ct_a 的数值做一些清理。CT 扫描的体素值使用的是 Hounsfield 单位（HU），这是一种有点特别的单位：空气是 -1000 HU（对我们的目的而言，这已经足够接近 0 g/cc，也就是每立方厘米 0 克），水是 0 HU（1 g/cc），而骨骼至少在 +1000 HU 以上（约 2–3 g/cc）。
有些 CT 扫描仪会把那些位于扫描仪视野之外的体素，用对应于“负密度”的 HU 值来表示。对我们来说，患者身体外部都应该被看作空气，因此我们会通过把下界设为 -1000 HU，来丢弃这些“视野范围”信息。同样地，骨骼、金属植入物等的精确密度，并不是我们这个任务所关心的内容，所以我们也会把密度上限封顶在大约 2 g/cc（1000 HU），尽管从生物学上说，这在多数情况下并不精确：
ct_a.clip(-1000, 1000, ct_a)
高于 0 HU 的数值与密度之间并不完全成比例，但我们关心的肿瘤通常大约在 1 g/cc（0 HU）附近，因此我们就忽略“HU 与常见密度单位并不完美对应”这个事实。这样做没有问题，因为我们的模型会直接以 HU 作为输入进行训练。
我们希望把这些离群值从数据中去掉：它们与我们的目标并不直接相关，而且这些异常值会让模型学习变得更困难。这种情况可能以很多方式出现，但一个常见例子是：batch normalization（第 8 章讨论过）如果吃进这些异常值，那么它用来估计“如何最好地归一化数据”的统计量就会被扭曲。要始终留意是否有机会对数据进行清洗。
到这里，我们构建出的这些值都会赋给 self：
self.series_uid = series_uid
self.hu_a = ct_a
知道我们的数据范围落在 -1000 到 +1000 之间，这一点非常重要。如果我们不考虑 HU 数值与其他附加数据之间的量纲差异，那么新加入的通道很容易被原始 HU 值淹没。有时候，构建一个好模型确实需要领域知识；理解 Hounsfield 单位是什么，以及它与被扫描物质密度之间的关系，就是一个很好的例子。
12.4 使用病人坐标系定位结节
深度学习模型通常需要固定大小的输入（当然也有例外，但这里暂时不相关），因为输入神经元的数量是固定的。我们需要能够生成一个包含候选区域的固定大小数组，以便把它作为分类器输入。我们希望用一个以候选为中心的 CT 局部裁剪块来训练模型，因为这样一来，模型就不必额外学习如何去注意那些躲在输入角落里的结节。通过减少预期输入的变化范围，我们让模型的任务变得更容易。
12.4.1 病人坐标系
不幸的是，我们在 12.2 节加载的所有候选中心数据，单位都是毫米，而不是体素！我们不能直接把“毫米位置”拿去当作数组索引，然后还指望一切按预期工作。正如图 12.5 所示，我们需要把当前以毫米为单位的坐标系 (X, Y, Z)，转换成 CT 数据数组切片所使用的体素寻址坐标系 (I, R, C)。这正是“必须一致地处理单位”这一原则的典型例子。
图 12.5 利用变换信息把候选坐标系中的结节中心坐标 (X, Y, Z) 转换为数组索引 (I, R, C)
正如我们前面提到过的，在处理 CT 扫描时，我们把数组维度称作 index、row、column，而不是 X、Y、Z，因为 X、Y、Z 在这里另有专门含义，如图 12.6 所示。病人坐标系规定：正 X 指向病人的左侧（left），正 Y 指向病人的后方（posterior），正 Z 指向病人的头部方向（superior）。Left-posterior-superior 有时缩写为 LPS。
图 12.6 我们这位穿着不合场景服装的病人，用来示意病人坐标系的各个轴向
病人坐标系以毫米为单位，其原点位置是任意定义的，并不对应于 CT 体素数组的原点，如图 12.7 所示。
图 12.7 数组坐标（左）和病人坐标（右）有着不同的原点和尺度
病人坐标系常被用于以一种不依赖于具体某一次扫描的方式，来描述感兴趣的解剖位置。定义 CT 数组与病人坐标系之间关系的元数据，存储在 DICOM 文件头中，而 MetaIO 格式也在其头信息里保留了这些数据。正是这些元数据，使我们能够构建图 12.5 中所示的从 (X, Y, Z) 到 (I, R, C) 的变换。
12.4.2 CT 扫描形状与体素大小
CT 扫描之间最常见的差异之一，就是体素大小；它们通常并不是立方体。相反，它们可能是 1.125 mm × 1.125 mm × 2.5 mm 这样的尺寸。通常 row 和 column 两个维度上的体素大小相同，而 index 维度上的值会更大，但也可能存在其他比例关系。
如果直接用方形像素来绘制这些非立方体体素，图像看起来就会有些变形，有点类似于使用墨卡托投影地图时，靠近南北极区域出现的那种变形（mng.bz/pZqE）。当然，这个类比并不完全恰当，因为这里的变形是均匀且线性的——图 12.8 中的病人从上到下看起来比现实中更加矮胖、胸腔更鼓。要想让图像呈现出真实比例，我们就需要施加一个缩放因子。
图 12.8 一张在 index 轴方向上使用非立方体体素的 CT 扫描。注意肺部从上到下被压缩得多么明显。
了解这些细节，在试图从视觉上解释结果时会很有帮助。如果没有这些知识，我们很容易误以为数据加载出了问题，可能会觉得图像之所以看起来那么“矮胖”，是不是因为我们不小心跳过了一半切片之类的。把时间浪费在调试一个其实一直正常工作的环节上，是很容易发生的事；而熟悉自己的数据，就能帮助你避免这种情况。
CT 图像通常是 512 行 × 512 列，而 index 维度大约从 100 张切片到 250 张切片不等（250 张切片乘以每张 2.5 毫米，通常已经足够覆盖感兴趣的解剖区域）。每个 CT 文件都会在元数据中指定体素大小（以毫米计）。
12.4.3 在毫米坐标与体素地址之间转换
我们将定义一些工具代码，用来在“以毫米表示的病人坐标”（在代码里我们会用 _xyz 后缀来标注变量）和“数组坐标 (I, R, C)”（代码里用 _irc 后缀表示）之间做转换。
你可能会问：SimpleITK 库不是已经提供了这类转换函数吗？确实，Image 实例提供了两个方法——TransformIndexToPhysicalPoint 和 TransformPhysicalPointToIndex——就是用来做这个转换的（只不过涉及从 CRI〔column,row,index〕到 IRC 的顺序调换）。不过，我们希望在不必一直保留 Image 对象的前提下完成这一计算，因此这里我们会手动进行数学变换。
翻转坐标轴（以及可能存在的旋转或其他变换）被编码在一个 3 × 3 矩阵中，该矩阵以 tuple 形式由 ct_mhd.GetDirections() 返回。要从体素索引走到物理坐标，我们需要按顺序执行以下四步：
- 把坐标顺序从 IRC 调整为 CRI，以便与 XYZ 对齐。
- 使用体素尺寸对索引进行缩放。
- 使用 Python 里的
@运算符，与 directions 矩阵做矩阵乘法。 - 加上原点偏移量。
若要从 XYZ 反向回到 IRC，则需要按相反顺序执行这些步骤的逆操作。
我们把体素大小保存在 named tuple 里，因此这里会先把它们转换成数组。
代码清单 12.2 在病人坐标（XYZ）与数组索引（IRC）之间进行转换（utils.py）
IrcTuple = collections.namedtuple('IrcTuple', ['index', 'row', 'col'])
XyzTuple = collections.namedtuple('XyzTuple', ['x', 'y', 'z'])
def irc2xyz(coord_irc, origin_xyz, vxSize_xyz, direction_a):
cri_a = np.array(coord_irc)[::-1] #1
origin_a = np.array(origin_xyz)
vxSize_a = np.array(vxSize_xyz)
coords_xyz = (direction_a @ (cri_a * vxSize_a)) + origin_a #2
return XyzTuple(*coords_xyz)
def xyz2irc(coord_xyz, origin_xyz, vxSize_xyz, direction_a):
origin_a = np.array(origin_xyz)
vxSize_a = np.array(vxSize_xyz)
coord_a = np.array(coord_xyz)
cri_a = ((coord_a - origin_a) @ np.linalg.inv(direction_a)) / vxSize_a #3
cri_a = np.round(cri_a) #4
return IrcTuple(int(cri_a[2]), int(cri_a[1]), int(cri_a[0])) #5
#1 在转换为 NumPy 数组时顺便交换顺序
#2 我们计划中的后面三步，在这一行里一次完成
#3 最后三步的逆操作
#4 在转成整数之前先做规范的四舍五入
#5 调整顺序并转为整数
呼——如果这一段看起来有点重，也别担心。只要记住，我们需要做坐标转换，而这些函数可以先当成黑盒来使用。把病人坐标（_xyz）转换成数组坐标（_irc）所需的元数据，就存放在 MetaIO 文件中，与 CT 数据本身放在一起。我们在读取 ct_a 的同时，也会从 .mhd 文件中提取体素尺寸和位置信息：
class Ct:
def __init__(self, series_uid):
import SimpleITK as sitk
mhd_path = glob.glob(
f'{mhd_data_folder}/{series_uid}.mhd'
)[0]
ct_mhd = sitk.ReadImage(mhd_path)
# ...
self.origin_xyz = XyzTuple(*ct_mhd.GetOrigin())
self.vxSize_xyz = XyzTuple(*ct_mhd.GetSpacing())
self.direction_a = np.array(ct_mhd.GetDirection()).reshape(3, 3) #1
#1 先把 directions 转成数组，再把这个九元素数组 reshape 成正确的 3 × 3 矩阵形状
这些属性就是我们在调用 xyz2irc 转换函数时，除了待转换点本身之外所需传入的参数。有了这些属性之后，我们的 Ct 对象实现就已经具备了把候选中心从病人坐标转换为数组坐标所需的全部数据。
12.4.4 从 CT 扫描中提取结节
正如我们在第 11 章中提到的，对于一个包含肺结节的病人 CT 扫描而言，其中高达 99.9999% 的体素都不会属于实际结节（更不用说癌变部分了）。再说一遍，这样的比例就相当于一台高清电视屏幕上某个地方有两个颜色不对的像素，或者一整架小说中只有一个拼错的单词。强迫模型在如此巨量的数据中寻找那些我们真正关心的结节线索，效果大概不会比让你在一堆你根本不认识的外文小说中找出那个拼写错误的单词更好！顺便问一句，你到现在在这本书里发现过拼错的单词了吗？;)
相反，正如图 12.9 所示，我们会从每个候选位置周围提取一个局部区域，让模型一次只关注一个候选。这就像是允许你只阅读那门外语中的单独一段——当然仍然不轻松，但至少没那么令人绝望！为模型寻找各种“缩小问题范围”的方法，通常会很有帮助，尤其是在项目初期，我们还在努力让第一个可运行实现落地的时候。
图 12.9 利用候选中心的数组坐标 (I, R, C) 信息，从更大的 CT 体素数组中裁剪出候选样本
getRawCandidate 函数接收一个用病人坐标系 (X, Y, Z) 表示的中心点——这正是 LUNA CSV 数据中的表示方式——以及一个以体素为单位的宽度。它会返回一个立方体形状的 CT 数据块，而候选中心则会被转换为数组坐标：
def getRawCandidate(self, center_xyz, width_irc):
center_irc = xyz2irc(
center_xyz,
self.origin_xyz,
self.vxSize_xyz,
self.direction_a,
)
slice_list = []
for axis, center_val in enumerate(center_irc):
start_ndx = int(round(center_val - width_irc[axis]/2))
end_ndx = int(start_ndx + width_irc[axis])
slice_list.append(slice(start_ndx, end_ndx))
ct_chunk = self.hu_a[tuple(slice_list)]
return ct_chunk, center_irc
真实的实现中还需要处理这样一些情况：候选中心与裁剪宽度的组合，会让裁剪区域的边缘落到数组之外。不过正如前面所说，我们会略过那些会掩盖函数主要意图的复杂细节。
12.5 一个直接明了的 Dataset 实现
我们第一次接触 PyTorch 的 Dataset 实例是在第 7 章，但这是我们第一次自己动手实现一个。通过继承 Dataset，我们就能把任意形式的数据接入到 PyTorch 生态系统的其余部分中。每一个 Ct 实例都代表了数百个不同的样本，而这些样本都可以用于训练模型，或者用于验证模型效果。
我们的 LunaDataset 类会对这些样本做归一化处理，把每个 CT 中的结节“摊平”成一个统一集合，使我们在取样本时不必关心它原本来自哪一个 Ct 实例。这种“摊平”方式往往正是我们希望处理数据的方式，不过我们很快也会看到：当类别不平衡或存在其他采样约束时，仅仅简单地摊平往往还不够。
在实现上，我们会先从 Dataset 子类必须满足的要求出发，然后再反推实现内容。这与之前我们使用的数据集不同；此前用到的都是外部库已经提供好的类，而这里我们需要自己实现并实例化这个类。做完之后，我们就能像使用之前那些数据集一样来使用它。幸运的是，这个自定义子类的实现并不难，因为 PyTorch API 只要求任何 Dataset 子类必须提供以下两个函数：
- 一个
__len__实现，它必须返回一个单一且固定的值，用来表示数据集长度（在某些场景中，这个值还会被缓存） - 一个
__getitem__方法，它接收一个索引，并返回一个 tuple，里面包含供训练使用的样本数据（或用于验证的数据）
先来看看这两个函数的签名和返回值长什么样：
def __len__(self):
return len(self.candidateInfo_list)
def __getitem__(self, ndx):
# ... line 200
return (
candidate_t, #1
pos_t, #1
candidateInfo_tup.series_uid, #1
torch.tensor(center_irc), #1
)
#1 这就是我们的训练样本。
我们的 __len__ 实现非常直接：我们有一个候选列表，每个候选就是一个样本，因此数据集大小就是样本数量。实现不一定非得像这里这样简单；后面几章你会看到它甚至会变得“更简单”一些。但重点在于，我们是有设计空间的。唯一规则是：如果 __len__ 返回 N，那么 __getitem__ 必须能对从 0 到 N-1 的所有输入都返回合法结果。
至于 __getitem__，它接收 ndx（根据前面的规则，通常就是整数），并返回一个由四项组成的 sample tuple，正如图 12.2 所示。这个四元组是当前我们这个数据集特有的；对于其他问题，它的内容则会不同。不过，构建这个 tuple 要比单纯求数据集长度复杂得多，所以我们来具体看看。
该方法的第一部分意味着：我们需要先构造 self.candidateInfo_list，并提供 getCtRawNodule 函数（这两部分我们很快会在 12.5.1 和 12.5.2 中讲到）：
def __getitem__(self, ndx):
candidateInfo_tup = self.candidateInfo_list[ndx]
width_irc = (32, 48, 48)
candidate_a, center_irc = getCtRawCandidate( #1
candidateInfo_tup.series_uid,
candidateInfo_tup.center_xyz,
width_irc,
)
#1 返回值 candidate_a 的形状是 (32, 48, 48)；三个轴分别对应 depth、height 和 width。
接下来，在 __getitem__ 方法中，我们还需要把数据整理成下游代码所期待的正确数据类型和数组维度：
candidate_t = torch.from_numpy(candidate_a)
candidate_t = candidate_t.to(torch.float32)
candidate_t = candidate_t.unsqueeze(0) #1
#1 .unsqueeze(0) 会增加一个 Channel 维。
现在你不必太纠结为什么这里要对维度做这些处理；下一章里会出现真正消费这个输出的代码，也正是那些代码预先施加了这些约束。对你自己今后实现的任何自定义 Dataset 来说，这种转换都应该是意料之中的一部分。这些变换正是把“蛮荒西部”般的现实数据，变成井井有条张量的关键环节。
最后，我们还需要构建分类张量：
pos_t = torch.tensor([
not candidateInfo_tup.isNodule_bool,
candidateInfo_tup.isNodule_bool
],
dtype=torch.long,
)
这个张量有两个元素，分别对应两种候选类别（非结节与结节）。我们当然也可以只用一个输出值来表示“结节状态”（阳性或阴性），但这里我们把它视为两个类别，而 nn.CrossEntropyLoss 期望每个类别对应一个输出值，因此我们就按这种方式来构造。你最终要构造什么样的张量，当然还是取决于你所做项目的具体类型。
我们来看一下最终的 sample tuple 长什么样（其中较大的 nodule_t 输出并不太适合直接阅读，因此清单中省略了大部分内容）。
代码清单 12.3 检查来自 LunaDataset 的训练样本（p2ch12_explore_data.ipynb）
# In[10]:
LunaDataset()[0]
# Out[10]:
(tensor([[[[-899., -903., -825., ..., -901., -898., -893.], #1
..., #1
[ -92., -63., 4., ..., 63., 70., 52.]]]]), #1
tensor([0, 1]), #2
'1.3.6...287966244644280690737019247886', #3
tensor([ 91, 360, 341])) #4
#1 candidate_t（形状为 [1, 32, 48, 48] 的张量，存储 CT 扫描数据，单位为 HU）
#2 pos_t（结节类、非结节类）
#3 candidate_tup.series_uid（该候选的截断唯一 ID）
#4 center_irc（候选中心的 IRC 坐标）
这里我们就看到了 __getitem__ 返回的四项内容。
12.5.1 使用 getCtRawCandidate 函数缓存候选数组
为了让 LunaDataset 获得还不错的性能，我们需要投入一点精力做磁盘缓存。这样就能避免每次取样本时都从磁盘读取整份 CT 扫描。如果每个样本都这样做，速度会慢到无法接受！你在自己的项目中也要注意性能瓶颈，并在它们真正拖慢系统时尽可能优化。严格说来，我们这里有点“抢跑”了，因为我们还没正式证明自己一定需要缓存。但没有缓存的话，LunaDataset 会慢大约 50 倍！我们会在本章练习中重新回到这个话题。
这个函数本身倒是很简单。它本质上是对前面见过的 Ct.getRawCandidate 方法包了一层基于文件缓存（pypi.python.org/pypi/diskca…）的封装：
@functools.lru_cache(1, typed=True)
def getCt(series_uid):
return Ct(series_uid)
@raw_cache.memoize(typed=True)
def getCtRawCandidate(series_uid, center_xyz, width_irc):
ct = getCt(series_uid)
ct_chunk, center_irc = ct.getRawCandidate(center_xyz, width_irc)
return ct_chunk, center_irc
这里我们用了几种不同的缓存方式。首先，我们把 getCt 的返回值缓存在内存里，这样如果多次请求同一个 Ct 实例，就不必一次次从磁盘重新加载全部数据了。对于重复请求而言，这能带来巨大的速度提升；不过，我们只在内存里保留一个 CT，因此如果访问顺序安排不好，cache miss 依然会相当频繁。
调用 getCt 的 getCtRawCandidate 函数，其输出也被缓存了，因此当缓存填充好之后，getCt 就几乎不会再被调用。这些值通过 Python 库 diskcache 缓存在磁盘上。至于为什么要设计成这种缓存结构，我们会在第 13 章讨论。眼下你只需要知道：从磁盘直接读取 2^15 个 float32 数值，远比读取 2^25 个 int16 数值、再把它们转换成 float32、然后从中选出一个长度为 2^15 的子集要快得多。自第二次遍历数据开始，输入阶段的 I/O 时间应该就会降到几乎可以忽略。
注意：如果这些函数的定义将来发生了实质性变化，我们就必须把磁盘中的缓存值删除。否则，即便新函数已经不再把给定输入映射到旧输出，缓存仍然会继续把旧结果返回给我们。这些数据存放在 data-unversioned/cache 目录下。
12.5.2 在 LunaDataset.__init__ 中构建数据集
几乎每一个项目都需要把样本分成训练集和验证集。这里我们会通过一个设计来实现：每隔 10 个样本取 1 个作为验证集成员，这由 val_stride 参数指定。我们还会接受一个 isValSet_bool 参数，用它来决定当前这个数据集实例到底保留训练数据、验证数据，还是全部数据：
class LunaDataset(Dataset):
def __init__(self,
val_stride=0,
isValSet_bool=None,
series_uid=None,
):
self.candidateInfo_list = copy.copy(getCandidateInfoList()) #1
if series_uid:
self.candidateInfo_list = [
x for x in self.candidateInfo_list if x.series_uid == series_uid
]
#1 对返回值做一次拷贝，以免修改 self.candidateInfo_list 时影响缓存中的原始副本
如果传入了一个真值的 series_uid，那么该实例就只会包含这个 series 中的结节。这对于可视化或调试很有帮助，因为这样更容易只盯着一个有问题的 CT 扫描看。
12.5.3 训练/验证划分
我们允许 Dataset 按照 1/N 的比例，把数据划分出一个子集专门用于验证模型。具体如何处理这个子集，取决于 isValSet_bool 参数的值：
if isValSet_bool:
assert val_stride > 0, val_stride
self.candidateInfo_list = self.candidateInfo_list[::val_stride]
assert self.candidateInfo_list
elif val_stride > 0:
del self.candidateInfo_list[::val_stride] #1
assert self.candidateInfo_list
#1 从 self.candidateInfo_list 中删除验证图像（也就是列表中每隔 val_stride 的元素）。因为前面做过拷贝，所以不会修改原始列表。
这样一来，我们就可以创建两个 Dataset 实例，并且确信训练数据与验证数据之间是严格隔离的。当然，这依赖于 self.candidateInfo_list 存在一个一致的排序顺序；而这一点，我们通过候选信息 tuple 的稳定排序，再加上 getCandidateInfoList 在返回前对列表排序来保证。
关于训练数据与验证数据的分离，还有另一个需要注意的点是：具体到某些任务时，我们可能需要确保同一个病人的数据只能出现在训练集中或测试集中，而不能两边同时出现。但在这里，这不是问题；否则，我们就需要在进入“结节级别”之前，先对病人和 CT 扫描进行划分。
我们来看一下 p2ch12_explore_data.ipynb 中的数据：
# In[2]:
from p2ch12.dsets import getCandidateInfoList, getCt, LunaDataset
candidateInfo_list = getCandidateInfoList(require_on_disk=False)
positiveInfo_list = [x for x in candidateInfo_list if x[0]]
diameter_list = [x[1] for x in positiveInfo_list]
# In[4]:
for i in range(0, len(diameter_list), 100):
print('{:4} {:4.1f} mm'.format(i, diameter_list[i]))
# Out[4]:
0 32.3 mm
100 17.7 mm
200 13.0 mm
300 10.0 mm
400 8.2 mm
500 7.0 mm
600 6.3 mm
700 5.7 mm
800 5.1 mm
900 4.7 mm
1000 4.0 mm
1100 0.0 mm
1200 0.0 mm
1300 0.0 mm
我们有少数几个非常大的候选，最大可达 32 mm，但很快就下降到它的一半左右。大多数候选的尺寸集中在 4 到 10 mm 之间，还有数百个候选根本没有尺寸信息。这和预期是一致的；你可能还记得，真实结节的数量比带直径标注的结节更多。对数据做这种快速的 sanity check 往往非常有帮助；及早发现问题或错误假设，可能为你节省数小时的时间！
更大的启示是：若要让训练集和验证集真正有效，它们应当具备以下几个性质：
- 两个集合都应当包含预期输入各种变化形式的样本。
- 除非是出于特定目的（例如训练模型对离群值具有鲁棒性），否则两个集合都不应包含那些不具代表性的样本。
- 训练集不应当向验证集泄露在真实世界数据中本不可能获得的“不公平提示”。（例如，把同一个样本，或同一个结点的不同样本，同时放进训练集和验证集，这种情况被称为训练集泄漏。）
12.5.4 渲染数据
再次提醒，你可以直接使用 p2ch12_explore_data.ipynb，或者启动 Jupyter Notebook 并输入：
# In[7]:
from p2ch12.vis import findPositiveSamples, showCandidate
positiveSample_list = findPositiveSamples()
# In[8]:
series_uid = positiveSample_list[11][2]
showCandidate(series_uid)
提示：关于 Jupyter 的 Matplotlib inline magic（这是他们的叫法，不是我们的！），可参考：mng.bz/rrmD。
这会生成类似于本章前面那些展示 CT 与结节切片的图像。
如果你有兴趣，我们也鼓励你根据自己的需求和审美，修改 p2ch12/vis.py 中渲染代码的实现。这些渲染代码大量使用了 Matplotlib（matplotlib.org），这个库本身过于庞大复杂，我们不打算在这里详细展开。
请记住，对数据进行渲染，并不仅仅是为了生成好看的图片。真正的目的是让你对输入数据长什么样建立直觉。能够一眼看出“这个有问题的样本比我的其他数据噪声大得多”，或者“奇怪，这个样本看起来其实很正常”，在排查问题时非常有用。有效的渲染还会激发出这样的洞见：“也许如果我按某种方式修改它，就能解决现在遇到的问题。” 当你开始处理越来越难的项目时，这种程度的熟悉感将会变得必不可少。
注意：由于各个子集的划分方式，再加上构建 LunaDataset.candidateInfo_list 时所采用的排序方式，noduleSample_list 中条目的顺序会高度依赖于代码执行时磁盘上到底有哪些子集存在。当你试图第二次找到某个特定样本时，尤其是在你后来又解压了更多子集的情况下，请记住这一点。
12.6 结论
在第 11 章中，我们让自己的脑子理解了数据。而在本章中，我们则让 PyTorch 的“脑子”也理解了数据！通过把基于 DICOM / Meta-image 的原始数据转换成张量，我们已经为实现模型和训练循环搭好了舞台，而这些内容会在下一章中出现。
不要低估我们目前为止已经做出的这些设计决策所带来的影响：输入尺寸的选择、缓存结构的设计、以及训练集和验证集的划分方式，都会影响整个项目最终是成功还是失败。以后，尤其是在你自己的项目里，一旦发现有必要，请不要犹豫，回头重新审视这些决定。
12.7 练习
实现一个程序，对某个 LunaDataset 实例进行迭代，并计时整个过程所需的时间。为了节省时间，或许可以加一个选项，把迭代限制在前 N=1000 个样本上：
- 第一次运行需要多长时间？
- 第二次运行需要多长时间？
- 清空缓存对运行时间有什么影响？
- 如果改为使用最后
N=1000个样本，第一次 / 第二次运行时间会怎样变化？
修改 LunaDataset 的实现，使其在 __init__ 中将样本列表随机化。清空缓存后再运行修改版。这样会对第一次和第二次运行时间造成什么影响？
撤销随机化，并注释掉 getCt 上的 @functools.lru_cache(1, typed=True) 装饰器。清空缓存并运行修改版。现在运行时间又会怎样变化？
小结
- 解析和加载原始数据所需的代码，往往并不简单。在这个项目中，我们实现了一个
Ct类，用来从磁盘加载数据，并提供围绕感兴趣点裁剪区域的访问能力。 - 如果解析和加载过程代价较高，那么缓存会很有用。要记住，有些缓存适合放在内存中，有些则更适合放在磁盘上。在数据加载流水线中，这两者各有其用武之地。
- PyTorch 的
Dataset子类用于把原始形态的数据转换为适合送入模型的张量。我们可以利用这一机制，把现实世界中的数据接入到 PyTorch API 中。 Dataset的子类需要提供两个方法的实现：__len__和__getitem__。其他辅助方法是允许的，但不是必须的。- 将数据划分为合理的训练集和验证集时，必须确保没有样本同时出现在两个集合中。这里我们通过保持一致的排序顺序，并每隔 10 个样本取 1 个放入验证集来实现这一点。
- 数据可视化非常重要；能够从视觉上检查数据，往往能为发现错误或问题提供关键线索。我们这里使用的是 Jupyter Notebook 和 Matplotlib 来渲染数据。
