使用 PyTorch 和 TIAToolbox 进行全玻片图像分类¶
创建时间： 2023年12月19日 |上次更新时间：2024 年 8 月 27 日 |上次验证： Nov 05， 2024
提示
为了充分利用本教程，我们建议使用此 Colab 版本。这将允许您试验下面提供的信息。
介绍¶
在本教程中，我们将展示如何对整个幻灯片图像 （WSI） 进行分类 在 TIAToolbox 的帮助下使用 PyTorch 深度学习模型。一个 WSI 是通过手术或活检采集的人体组织样本的图像，并且 使用专用扫描仪进行扫描。它们被病理学家和 计算病理学研究人员在微观上研究癌症等疾病 级别输入 例如，为了了解肿瘤生长并帮助改善治疗 为患者。
WSI 难以处理的原因是它们的巨大规模。为 例如，典型的幻灯片图像的大小约为 100,000x100,000 像素，其中每个像素都可以 对应于载玻片上的约 0.25x0.25 微米。这引入了 加载和处理此类图像的挑战，更不用说 单个研究中有数百甚至数千个 WSI（大型研究 产生更好的结果）！
传统的图像处理管道不适用于 WSI 处理，因此我们需要更好的工具。这就是 TIAToolbox 可以的地方 帮助，因为它带来了一组有用的工具来导入和处理组织 幻灯片。通常，WSI 保存在金字塔结构中，其中包含同一图像的多个副本 在各种放大倍率级别下，针对可视化进行了优化。0 级 （或最底层）包含最高层的图像 放大或缩放级别，而金字塔中的较高级别 具有较低分辨率的基础映像副本。金字塔结构是 草图如下。
WSI 金字塔堆栈
（来源)
TIAToolbox 允许我们自动执行常见的下游分析任务，例如
作为 组织
分类。在这个
教程中，我们将展示如何做到这一点：1. 使用
TIAToolbox;和 2.使用不同的 PyTorch 模型对幻灯片进行分类
补丁级别。在本教程中，我们将提供一个使用
TorchVision 模型和自定义 HistoEncoder <https://github.com/jopo666/HistoEncoder>'__ 模型。ResNet18
让我们开始吧！
设置环境¶
要运行本教程中提供的示例，以下软件包 是必需的先决条件。
OpenJpeg 格式
打开幻灯片
皮克斯曼
TIAToolbox
HistoEncoder （用于自定义模型示例）
请在您的终端中运行以下命令来安装这些 包：
apt-get -y -qq install libopenjp2-7-dev libopenjp2-tools openslide-tools libpixman-1-dev pip install -q 'tiatoolbox<1.5' histoencoder && echo “安装完成。”
或者，您可以运行
在 MacOS 上安装必备软件包，而不是 。
有关安装的更多信息，请参见
这里。brew install openjpeg openslideapt-get
运行前清理¶
为确保正确清理（例如在异常终止时），所有
在此运行中下载或创建的文件保存在一个目录中 ，我们将其设置为等于“./tmp/”。简化
maintenance，则目录的名称只出现在这一个位置，因此
如果需要，可以很容易地更改它。global_save_dir
warnings.filterwarnings("ignore")
global_save_dir = Path("./tmp/")
def rmdir(dir_path: str | Path) -> None:
"""Helper function to delete directory."""
if Path(dir_path).is_dir():
shutil.rmtree(dir_path)
logger.info("Removing directory %s", dir_path)
rmdir(global_save_dir) # remove directory if it exists from previous runs
global_save_dir.mkdir()
logger.info("Creating new directory %s", global_save_dir)
下载数据¶
对于我们的示例数据，我们将使用一张全玻片图像，以及来自 Kather 的验证子集 100k 数据集。
wsi_path = global_save_dir / "sample_wsi.svs"
patches_path = global_save_dir / "kather100k-validation-sample.zip"
weights_path = global_save_dir / "resnet18-kather100k.pth"
logger.info("Download has started. Please wait...")
# Downloading and unzip a sample whole-slide image
download_data(
"https://tiatoolbox.dcs.warwick.ac.uk/sample_wsis/TCGA-3L-AA1B-01Z-00-DX1.8923A151-A690-40B7-9E5A-FCBEDFC2394F.svs",
wsi_path,
)
# Download and unzip a sample of the validation set used to train the Kather 100K dataset
download_data(
"https://tiatoolbox.dcs.warwick.ac.uk/datasets/kather100k-validation-sample.zip",
patches_path,
)
with ZipFile(patches_path, "r") as zipfile:
zipfile.extractall(path=global_save_dir)
# Download pretrained model weights for WSI classification using ResNet18 architecture
download_data(
"https://tiatoolbox.dcs.warwick.ac.uk/models/pc/resnet18-kather100k.pth",
weights_path,
)
logger.info("Download is complete.")
读取数据¶
我们创建一个补丁列表和相应的标签列表。为
example 中，中的第一个标签将指示
中的第一个图像块。label_listpatch_list
# Read the patch data and create a list of patches and a list of corresponding labels
dataset_path = global_save_dir / "kather100k-validation-sample"
# Set the path to the dataset
image_ext = ".tif" # file extension of each image
# Obtain the mapping between the label ID and the class name
label_dict = {
"BACK": 0, # Background (empty glass region)
"NORM": 1, # Normal colon mucosa
"DEB": 2, # Debris
"TUM": 3, # Colorectal adenocarcinoma epithelium
"ADI": 4, # Adipose
"MUC": 5, # Mucus
"MUS": 6, # Smooth muscle
"STR": 7, # Cancer-associated stroma
"LYM": 8, # Lymphocytes
}
class_names = list(label_dict.keys())
class_labels = list(label_dict.values())
# Generate a list of patches and generate the label from the filename
patch_list = []
label_list = []
for class_name, label in label_dict.items():
dataset_class_path = dataset_path / class_name
patch_list_single_class = grab_files_from_dir(
dataset_class_path,
file_types="*" + image_ext,
)
patch_list.extend(patch_list_single_class)
label_list.extend([label] * len(patch_list_single_class))
# Show some dataset statistics
plt.bar(class_names, [label_list.count(label) for label in class_labels])
plt.xlabel("Patch types")
plt.ylabel("Number of patches")
# Count the number of examples per class
for class_name, label in label_dict.items():
logger.info(
"Class ID: %d -- Class Name: %s -- Number of images: %d",
label,
class_name,
label_list.count(label),
)
# Overall dataset statistics
logger.info("Total number of patches: %d", (len(patch_list)))
|2023-11-14|13:15:59.299| [INFO] Class ID: 0 -- Class Name: BACK -- Number of images: 211
|2023-11-14|13:15:59.299| [INFO] Class ID: 1 -- Class Name: NORM -- Number of images: 176
|2023-11-14|13:15:59.299| [INFO] Class ID: 2 -- Class Name: DEB -- Number of images: 230
|2023-11-14|13:15:59.299| [INFO] Class ID: 3 -- Class Name: TUM -- Number of images: 286
|2023-11-14|13:15:59.299| [INFO] Class ID: 4 -- Class Name: ADI -- Number of images: 208
|2023-11-14|13:15:59.299| [INFO] Class ID: 5 -- Class Name: MUC -- Number of images: 178
|2023-11-14|13:15:59.299| [INFO] Class ID: 6 -- Class Name: MUS -- Number of images: 270
|2023-11-14|13:15:59.299| [INFO] Class ID: 7 -- Class Name: STR -- Number of images: 209
|2023-11-14|13:15:59.299| [INFO] Class ID: 8 -- Class Name: LYM -- Number of images: 232
|2023-11-14|13:15:59.299| [INFO] Total number of patches: 2000
正如你所看到的，对于这个补丁数据集，我们有 9 个具有 ID 的类/标签 0-8 和关联的类名。描述 补丁：
返回 ⟶ 背景（空玻璃区域）
LYM ⟶ 淋巴细胞
NORM ⟶ 正常结肠粘膜
DEB ⟶ 碎片
MUS ⟶ 平滑肌
STR ⟶ 癌症相关基质
ADI ⟶ 脂肪
粘液 ⟶ 粘液
TUM ⟶ 结直肠腺癌上皮
对图像补丁进行分类¶
我们演示了如何获取
数字滑轨首先使用模式，然后是大滑轨
using 模式。patchwsi
定义模型PatchPredictor¶
PatchPredictor 类运行用 PyTorch 编写的基于 CNN 的分类器。
model可以是任何经过训练的 PyTorch 模型，其约束是 它应该遵循 （docs） <https://tia-toolbox.readthedocs.io/en/latest/_autosummary/tiatoolbox.models.models_abc.ModelABC.html>'__ 类结构。有关此问题的更多信息，请参阅我们的 Advanced Model 示例 Notebook 技术。 为了加载自定义模型，您需要编写一个小型的 预处理函数，如 ，它确保 输入张量的格式适合加载的网络。tiatoolbox.models.abc.ModelABCpreproc_func(img)或者，您也可以以字符串形式传递 论点。这指定了执行预测的 CNN 模型 并且它必须是此处列出的型号之一。 该命令将如下所示：.
pretrained_modelpredictor = PatchPredictor(pretrained_model='resnet18-kather100k', pretrained_weights=weights_path, batch_size=32)pretrained_weights：使用 时， 默认情况下，还将下载相应的预训练权重。 您可以通过参数使用自己的一组权重覆盖默认值。pretrained_modelpretrained_weightbatch_size：每次输入到模型中的图像数。高等 此参数的值需要更大的 （GPU） 内存容量。
# Importing a pretrained PyTorch model from TIAToolbox
predictor = PatchPredictor(pretrained_model='resnet18-kather100k', batch_size=32)
# Users can load any PyTorch model architecture instead using the following script
model = vanilla.CNNModel(backbone="resnet18", num_classes=9) # Importing model from torchvision.models.resnet18
model.load_state_dict(torch.load(weights_path, map_location="cpu", weights_only=True), strict=True)
def preproc_func(img):
img = PIL.Image.fromarray(img)
img = transforms.ToTensor()(img)
return img.permute(1, 2, 0)
model.preproc_func = preproc_func
predictor = PatchPredictor(model=model, batch_size=32)
预测补丁标签¶
我们创建一个预测器对象，然后使用
模式。然后我们计算分类准确率和
混淆矩阵。predictpatch
with suppress_console_output():
output = predictor.predict(imgs=patch_list, mode="patch", on_gpu=ON_GPU)
acc = accuracy_score(label_list, output["predictions"])
logger.info("Classification accuracy: %f", acc)
# Creating and visualizing the confusion matrix for patch classification results
conf = confusion_matrix(label_list, output["predictions"], normalize="true")
df_cm = pd.DataFrame(conf, index=class_names, columns=class_names)
df_cm
|2023-11-14|13:16:03.215| [INFO] Classification accuracy: 0.993000
预测整个幻灯片的补丁标签¶
现在，我们引入 ，一个指定
为模型配置图像读取和预测写入
预测引擎。这是通知分类器哪个级别所必需的
分类器应读取、处理数据并生成
输出。IOPatchPredictorConfig
的参数定义为：IOPatchPredictorConfig
input_resolutions： 一个列表，以字典的形式， 指定每个输入的分辨率。List 元素必须位于 与 Target 中的顺序相同。如果您的模型 只接受一个输入，你只需要放一个字典 指定 和 .请注意，TIAToolbox 支持具有多个输入的模型。有关更多信息 单位和分辨率，请参阅 TIAToolbox 文档。model.forward()'units''resolution'patch_input_shape：最大输入的形状 （height， width） 格式。stride_shape：两个之间的步幅（步长）的大小 Consecutive Patchs，用于补丁提取过程。如果 user 设置等于 、 patches 将被提取和处理，没有任何重叠。stride_shapepatch_input_shape
wsi_ioconfig = IOPatchPredictorConfig(
input_resolutions=[{"units": "mpp", "resolution": 0.5}],
patch_input_shape=[224, 224],
stride_shape=[224, 224],
)
该方法将 CNN 应用于输入补丁并获取
结果。以下是参数及其说明：predict
mode：要处理的输入类型。从 中选择，或根据您的应用程序进行选择。patchtilewsiimgs：输入列表，应该是指向 输入磁贴或 WSI。return_probabilities：设置为 True 可按类获取 probabilities 以及 input patches 的预测标签。如果你 希望合并预测以生成 OR 模式的预测图，则可以设置 。tilewsireturn_probabilities=Trueioconfig：使用类设置 IO 配置信息。IOPatchPredictorConfigresolution和 （未在下面显示）：这些参数 指定 WSI 级别的级别或每像素微米分辨率 我们计划从中提取补丁，并且可以代替 .在这里，我们将 WSI 级别指定为 ， 相当于 0 级。一般来说，这是 最高分辨率。在这种特殊情况下，图像只有一个 水平。有关更多信息，请参阅文档。unitioconfig'baseline'masks：与列表中 WSI 的掩码相对应的路径列表。这些掩码指定原始 WSI 中的区域 我们想从中提取补丁。如果特定 WSI 指定为 ，则该 的所有补丁的标签 WSI（甚至背景区域）将被预测。这可能会导致 不必要的计算。imgsNonemerge_predictions：您可以将此参数设置为如果 生成补丁分类结果的 2D 地图所需的。 但是，对于大型 WSI，这将需要较大的可用内存。一 替代（默认）解决方案是设置 ， 然后使用该函数生成 2D 预测图，稍后您将看到。Truemerge_predictions=Falsemerge_predictions
由于我们使用的是大型 WSI，因此补丁提取和预测
进程可能需要一些时间（请确保将 if
您可以访问启用 Cuda 的 GPU 和 PyTorch+Cuda）。ON_GPU=True
with suppress_console_output():
wsi_output = predictor.predict(
imgs=[wsi_path],
masks=None,
mode="wsi",
merge_predictions=False,
ioconfig=wsi_ioconfig,
return_probabilities=True,
save_dir=global_save_dir / "wsi_predictions",
on_gpu=ON_GPU,
)
我们通过以下方式了解预测模型如何对整个幻灯片图像进行处理
可视化 .我们首先需要合并补丁预测
输出，然后将它们可视化为原始图像上的叠加层。如
before，该方法用于合并补丁
预测。这里我们设置参数以生成预测地图
1.25 倍放大倍率。如果您想要更高/更低的分辨率
（更大/更小）预测地图，您需要更改这些参数
因此。合并预测后，使用该函数将预测地图叠加在
WSI 缩略图，应以用于
预测合并。wsi_outputmerge_predictionsresolution=1.25, units='power'overlay_patch_prediction
overview_resolution = (
4 # the resolution in which we desire to merge and visualize the patch predictions
)
# the unit of the `resolution` parameter. Can be "power", "level", "mpp", or "baseline"
overview_unit = "mpp"
wsi = WSIReader.open(wsi_path)
wsi_overview = wsi.slide_thumbnail(resolution=overview_resolution, units=overview_unit)
plt.figure(), plt.imshow(wsi_overview)
plt.axis("off")
将预测图叠加在此图像上，如下所示：
# Visualization of whole-slide image patch-level prediction
# first set up a label to color mapping
label_color_dict = {}
label_color_dict[0] = ("empty", (0, 0, 0))
colors = cm.get_cmap("Set1").colors
for class_name, label in label_dict.items():
label_color_dict[label + 1] = (class_name, 255 * np.array(colors[label]))
pred_map = predictor.merge_predictions(
wsi_path,
wsi_output[0],
resolution=overview_resolution,
units=overview_unit,
)
overlay = overlay_prediction_mask(
wsi_overview,
pred_map,
alpha=0.5,
label_info=label_color_dict,
return_ax=True,
)
plt.show()
使用病理特异性模型进行特征提取¶
在本节中，我们将展示如何从预训练的 存在于 TIAToolbox 外部的 PyTorch 模型，使用 WSI 推理 由 TIAToolbox 提供的引擎。为了说明这一点，我们将使用 HistoEncoder，一种特定于计算病理学的模型，已经 以自我监督的方式进行训练，以从组织学中提取特征 图像。该模型已在此处提供：
“HistoEncoder：数字病理学的基础模型” （https://github.com/jopo666/HistoEncoder） 作者：Pohjonen、Joona 和 赫尔辛基大学。
我们将 umap 缩减绘制为特征图的 3D （RGB） 以 可视化特征如何捕获某些 上述组织类型。
# Import some extra modules
import histoencoder.functional as F
import torch.nn as nn
from tiatoolbox.models.engine.semantic_segmentor import DeepFeatureExtractor, IOSegmentorConfig
from tiatoolbox.models.models_abc import ModelABC
import umap
TIAToolbox 定义了一个 ModelABC，它是一个继承 PyTorch nn.Module 并指定模型的外观，以便在 TIAToolbox 推理引擎。histoencoder 模型不遵循此 结构，因此我们需要将其包装在一个 output 和 methods 为 TIAToolbox 引擎所期望的。
class HistoEncWrapper(ModelABC):
"""Wrapper for HistoEnc model that conforms to tiatoolbox ModelABC interface."""
def __init__(self: HistoEncWrapper, encoder) -> None:
super().__init__()
self.feat_extract = encoder
def forward(self: HistoEncWrapper, imgs: torch.Tensor) -> torch.Tensor:
"""Pass input data through the model.
Args:
imgs (torch.Tensor):
Model input.
"""
out = F.extract_features(self.feat_extract, imgs, num_blocks=2, avg_pool=True)
return out
@staticmethod
def infer_batch(
model: nn.Module,
batch_data: torch.Tensor,
*,
on_gpu: bool,
) -> list[np.ndarray]:
"""Run inference on an input batch.
Contains logic for forward operation as well as i/o aggregation.
Args:
model (nn.Module):
PyTorch defined model.
batch_data (torch.Tensor):
A batch of data generated by
`torch.utils.data.DataLoader`.
on_gpu (bool):
Whether to run inference on a GPU.
"""
img_patches_device = batch_data.to('cuda') if on_gpu else batch_data
model.eval()
# Do not compute the gradient (not training)
with torch.inference_mode():
output = model(img_patches_device)
return [output.cpu().numpy()]
现在我们有了包装器，我们将创建我们的特征提取 建模并实例化一个 DeepFeatureExtractor，以允许我们在 WSI 上使用这个模型。我们将使用与 上图，但这次我们将从 WSI 使用 HistoEncoder 模型，而不是预测 每个补丁。
# create the model
encoder = F.create_encoder("prostate_medium")
model = HistoEncWrapper(encoder)
# set the pre-processing function
norm=transforms.Normalize(mean=[0.662, 0.446, 0.605],std=[0.169, 0.190, 0.155])
trans = [
transforms.ToTensor(),
norm,
]
model.preproc_func = transforms.Compose(trans)
wsi_ioconfig = IOSegmentorConfig(
input_resolutions=[{"units": "mpp", "resolution": 0.5}],
patch_input_shape=[224, 224],
output_resolutions=[{"units": "mpp", "resolution": 0.5}],
patch_output_shape=[224, 224],
stride_shape=[224, 224],
)
当我们创建 时，我们将传递参数。这将自动创建一个
使用 OTSU 阈值对组织区域进行掩膜，以便提取器
仅处理那些包含组织的补丁。DeepFeatureExtractorauto_generate_mask=True
# create the feature extractor and run it on the WSI
extractor = DeepFeatureExtractor(model=model, auto_generate_mask=True, batch_size=32, num_loader_workers=4, num_postproc_workers=4)
with suppress_console_output():
out = extractor.predict(imgs=[wsi_path], mode="wsi", ioconfig=wsi_ioconfig, save_dir=global_save_dir / "wsi_features",)
这些功能可用于训练下游模型，但在 为了获得对这些功能所代表的含义的一些直觉，我们将使用 UMAP 缩减，用于在 RGB 空间中可视化特征。要点 以相似颜色标记的 label 应该具有相似的特征，因此我们可以检查 如果特征自然分离到不同的组织区域 当我们在 WSI 缩略图上叠加 UMAP 缩减时。我们将绘制它 以及上面的补丁级预测图，以查看 特征与以下单元格中的补丁级预测进行了比较。
# First we define a function to calculate the umap reduction
def umap_reducer(x, dims=3, nns=10):
"""UMAP reduction of the input data."""
reducer = umap.UMAP(n_neighbors=nns, n_components=dims, metric="manhattan", spread=0.5, random_state=2)
reduced = reducer.fit_transform(x)
reduced -= reduced.min(axis=0)
reduced /= reduced.max(axis=0)
return reduced
# load the features output by our feature extractor
pos = np.load(global_save_dir / "wsi_features" / "0.position.npy")
feats = np.load(global_save_dir / "wsi_features" / "0.features.0.npy")
pos = pos / 8 # as we extracted at 0.5mpp, and we are overlaying on a thumbnail at 4mpp
# reduce the features into 3 dimensional (rgb) space
reduced = umap_reducer(feats)
# plot the prediction map the classifier again
overlay = overlay_prediction_mask(
wsi_overview,
pred_map,
alpha=0.5,
label_info=label_color_dict,
return_ax=True,
)
# plot the feature map reduction
plt.figure()
plt.imshow(wsi_overview)
plt.scatter(pos[:,0], pos[:,1], c=reduced, s=1, alpha=0.5)
plt.axis("off")
plt.title("UMAP reduction of HistoEnc features")
plt.show()
我们看到，来自补丁级预测器的预测图，以及 来自我们的自我监督特征编码器的特征映射，捕获类似的 有关 WSI 中组织类型的信息。这是一个很好的理智 检查我们的模型是否按预期工作。它还表明， HistoEncoder 模型提取的特征正在捕获 组织类型之间的差异，以便它们编码 组织学相关信息。
从这里去哪里¶
在本笔记本中，我们将展示如何使用 and 类及其方法来预测
大图块和 WSI 的 label 或 extract 特征。我们
introduce 和 helper
合并补丁预测输出并可视化
生成的预测图作为输入图像/WSI 上的叠加层。PatchPredictorDeepFeatureExtractorpredictmerge_predictionsoverlay_prediction_mask
所有过程都在 TIAToolbox 中进行，我们可以轻松地将
按照我们的示例代码拼凑在一起。请务必设置
inputs 和 options 正确。我们鼓励您进一步调查
更改函数对预测输出的影响
参数。我们已经演示了如何使用您自己的预训练模型或
一个由研究社区为
TIAToolbox 框架对大型 WSI 进行推理，即使模型
结构未在 TIAToolbox 模型类中定义。predict
您可以通过以下资源了解更多信息：