2021年认证杯SPSSPRO杯数学建模
A题 医学图像的配准
原题再现：
图像的配准是图像处理领域中的一个典型问题和技术难点，其目的在于比较或融合同一对象在不同条件下获取的图像。例如为了更好地综合多种信息来辨识不同组织或病变，医生可能使用多种仪器对患者的同一部位进行成像。在综合多幅图像时，首先需要将它们严格对齐，使得图上同一个坐标的位置对应的是真实对象的同一个点，这个过程称之为配准。现在的许多医学成像技术，包括 CT、MRI、PET 等，最终生成的是人体的断层影像。在这里，我们主要关心的是断层成像的配准问题。
我们考虑对一个患者的腹部进行断层成像。由于人体组织是柔软的，所以即使使用同一台成像设备，两次成像的结果也并不完全一致。最终输出时还会对图像进行自动放缩，所以输出图片的大小也并不完全相同。想要精确配准，需要将其中一次的成像结果进行某种仿射变换（或非线性变换），以尽可能地匹配另一次的结果（或将两次结果都映射到同一个标准模板中）。求得合适的变换就是图像配准的核心任务。
第二阶段问题：多模态的配准是指对来自不同设备的图像进行配准，例如对CT和MRI图像进行配准。有的组织或病变部位在单一的成像技术下与周边组织的区分不明显，所以我们可以对多种成像技术得到的图片进行融合处理，让每个像素点的成像结果表现为一个多维向量（每个分量都是一种成像技术的成像结果），这样可以更好地识别组织或病变的细节。多模态的配准则是图像融合的第一步。
现在我们有对患者同一身体部位（同一时间）的CT、MRI和PET成像结果。但在图像融合处理时遇到了两个问题：首先，每一种成像技术对不同组织的区分能力是不同的，例如有些不同的组织在CT下看起来区别不大，但在MRI下区分却十分明显；另外一些组织在MRI下区别不大，但在CT下区分却十分明显。所以对不同的成像设备而言，即使是同一个位置的成像结果，也并非完全相似，这给配准带来了难度。第二，在进行断层成像时，虽然对每个设备而言，我们能够确切地知道每个断层的位置，但不同设备扫描的断层位置并不完全相同。请你设计一个有效的方法，对这样的成像结果进行图像的融合。
整体求解过程概述(摘要)
一种基于多模态医学图像的图像融合方法可以显著提高融合图像的质量。一种有效的图像融合技术通过保留从源图像中收集到的所有可行的和显著的信息而不引入任何缺陷或不必要的扭曲来产生输出图像。
大多数深度学习方法采用所谓的单流高到低、低到高的网络结构，能够获得满意的整体配准结果。然而，一些严重变形的局部区域的精确定位，这是精确定位手术目标的关键，却经常被忽视。因此，这些方法对一些难以对齐的区域不敏感，例如畸形肝叶的图像配准融合。针对这一问题，我们提出一种新的无监督配准网络，即全分辨率残差配准网络(F3RNet)，用于严重变形器官的变形配准。该方法以残差学习的方式结合了两个并行处理流。一个流利用全分辨率信息，促进准确的体素级注册。另一个流学习深度多尺度残差表示以获得鲁棒特征识别与提取。此外，我们还使用了提升小波变换，混合融合等方法对图像特征进行精确分类。最后对三维卷积进行因式分解，得到图像配准与融合结果。我们选用了腹部和肺部的CT-MRI数据集对所提方法进行验证，实验结果表明所提方法可以获取更高质量的配准融合图像，同时又能显著提高配准融合效率。
为了验证我们提出方法的有效性，我们选用自备的医学CT图像和MRI图像数据集进行了评估，本文配准模型在肺部图像与脑部图像的配准值有一定程度的提升。从模型的运行的结果可以看出，当AUC值为0.883时，得到Dice稀疏为0.575，准确率为0.884，灵敏度为0.647，特异性为0.929，F1得分为0.640。模型中还进行了病变程度与匹配程度的相关性分析，在保证快捷的同时方便简洁地解决了问题，具有一定的推广意义。
问题分析：
考虑到在进行疾病诊断时，不同的医学图像在进行疾病诊断时所发挥的作用时不同的，为更全面的了解病人的病情，往往需要多种医学图像提供不同的信息，通过这些信息的综合分析，能够为临床诊断治疗提供全面的信息。深度学习的发展为多模态医学配准和图像融合提供了新的思路，利用深度学习的模型寻找待配准图像对像素点间的空间对应关系，使固定图像与浮动图像的像素点在空间解剖结构上对齐，根据对齐关系进行利用融合特性对这些提取的特征进行融合，能够更清晰的反应病变的情况。该方案的具体解决措施如下：
首先，提出一种全分辨率残差配准网络用于配准不同医学图像中的相关特征，通过设计全分辨率和多尺度残差块这两个并行的网络，全分辨率网络在密集网络表现较好，能够通过常规残差块提取不同医学图像中的低阶特征。多尺度残差块利用连续池化和卷积操作来增加识别范围，善于捕捉高级特征，从而提高识别性能。这两个并行网络能够全面提取不同医学图像的特征，不会因为医学图像分辨率不同导致特征提取不全面的情况。
其次，在提取完全特征的基础上，使用小波变换融合和混合图像融合相结合的方法进行医学图像的融合。小波变换融合利用小波变换对源参考的医学图像进行分解，利用定量度量度量技术的性能，并比较技术的效率，以找到合适的融合规则。再使用混合图像融合技术可以从有噪声、失真的图像中提取原始图像特征，从而获得改进和增强的图像质量。
模型假设：
在数学建模的过程中，为了使模型简单明确，简化运算过程，在不影响模型意义与计算精度的前提下，建立了如下假设：
(1)假设同一组的医学图像来自于同一位患者。
(2)假设用来采集每一张医学图像的设备都能正常采集医学图像。
(3)假设采集到的医学图像没有人为修改与损坏，并且能够直接使用。
(4)假设同一组诊断图像为患者同一部位的医学图像。
(5)假设患者在两次成像之间身体没有发生其他部分的病变，或者存在导致成像偏差的身体异状。
(6)假设每台计算机辅助医疗系统的参数设定与配置都保持一致。
论文缩略图：
程序代码：
from typing import Union, Optional, List, Tuple, Text, BinaryIO
import pathlib
import torch
import math
import warnings
import numpy as np
from PIL import Image, ImageDraw, ImageFont, ImageColor
__all__=["make_grid","save_image", "draw_bounding_boxes", "draw_segmentation_masks"]
def make_grid(
tensor: Union[torch.Tensor, List[torch.Tensor]],
nrow: int = 8,
padding: int = 2,
normalize: bool = False,
value_range: Optional[Tuple[int, int]] = None,
scale_each: bool = False,
pad_value: int = 0,
**kwargs
) -> torch.Tensor:
if not (torch.is_tensor(tensor) or
(isinstance(tensor, list) and all(torch.is_tensor(t) for t in tensor))):
raise TypeError(f'tensor or list of tensors expected, got {type(tensor)}')
if "range" in kwargs.keys():
warning = "range will be deprecated, please use value_range instead."
warnings.warn(warning)
value_range = kwargs["range"]
# if list of tensors, convert to a 4D mini-batch Tensor
if isinstance(tensor, list):
tensor = torch.stack(tensor, dim=0)
if tensor.dim() == 2: # single image H x W
tensor = tensor.unsqueeze(0)
if tensor.dim() == 3: # single image
if tensor.size(0) == 1: # if single-channel, convert to 3-channel
tensor = torch.cat((tensor, tensor, tensor), 0)
tensor = tensor.unsqueeze(0)
if tensor.dim() == 4 and tensor.size(1) == 1: # single-channel images
tensor = torch.cat((tensor, tensor, tensor), 1)
if normalize is True:
tensor = tensor.clone() # avoid modifying tensor in-place
if value_range is not None:
assert isinstance(value_range, tuple), \
"value_range has to be a tuple (min, max) if specified. min and max are numbers"
def norm_ip(img, low, high):
img.clamp_(min=low, max=high)
img.sub_(low).div_(max(high - low, 1e-5))
def norm_range(t, value_range):
if value_range is not None:
norm_ip(t, value_range[0], value_range[1])
else:
norm_ip(t, float(t.min()), float(t.max()))
if scale_each is True:
for t in tensor: # loop over mini-batch dimension
norm_range(t, value_range)
else:
norm_range(tensor, value_range)
if tensor.size(0) == 1:
return tensor.squeeze(0)
# make the mini-batch of images into a grid
nmaps = tensor.size(0)
xmaps = min(nrow, nmaps)
ymaps = int(math.ceil(float(nmaps) / xmaps))
height, width = int(tensor.size(2) + padding), int(tensor.size(3) + padding)
num_channels = tensor.size(1)
grid = tensor.new_full((num_channels, height * ymaps + padding, width * xmaps + padding), pad_value)
k = 0
for y in range(ymaps):
for x in range(xmaps):
if k >= nmaps:
break
# Tensor.copy_() is a valid method but seems to be missing from the stubs
# https://pytorch.org/docs/stable/tensors.html#torch.Tensor.copy_
grid.narrow(1, y * height + padding, height - padding).narrow( # type: ignore[attr-defined]
2, x * width + padding, width - padding
).copy_(tensor[k])
k = k + 1
return grid
@torch.no_grad()
def save_image(
tensor: Union[torch.Tensor, List[torch.Tensor]],
fp: Union[Text, pathlib.Path, BinaryIO],
format: Optional[str] = None,
**kwargs
) -> None:
grid = make_grid(tensor, **kwargs)
# Add 0.5 after unnormalizing to [0, 255] to round to nearest integer
ndarr = grid.mul(255).add_(0.5).clamp_(0, 255).permute(1, 2, 0).to('cpu', torch.uint8).numpy()
im = Image.fromarray(ndarr)
im.save(fp, format=format)
def draw_bounding_boxes(
image: torch.Tensor,
boxes: torch.Tensor,
labels: Optional[List[str]] = None,
colors: Optional[List[Union[str, Tuple[int, int, int]]]] = None,
fill: Optional[bool] = False,
width: int = 1,
font: Optional[str] = None,
font_size: int = 10
) -> torch.Tensor:
if not isinstance(image, torch.Tensor):
raise TypeError(f"Tensor expected, got {type(image)}")
elif image.dtype != torch.uint8:
raise ValueError(f"Tensor uint8 expected, got {image.dtype}")
elif image.dim() != 3:
raise ValueError("Pass individual images, not batches")
ndarr = image.permute(1, 2, 0).numpy()
img_to_draw = Image.fromarray(ndarr)
img_boxes = boxes.to(torch.int64).tolist()
if fill:
draw = ImageDraw.Draw(img_to_draw, "RGBA")
else:
draw = ImageDraw.Draw(img_to_draw)
txt_font = ImageFont.load_default() if font is None else ImageFont.truetype(font=font, size=font_size)
for i, bbox in enumerate(img_boxes):
if colors is None:
color = None
else:
color = colors[i]
if fill:
if color is None:
fill_color = (255, 255, 255, 100)
elif isinstance(color, str):
# This will automatically raise Error if rgb cannot be parsed.
fill_color = ImageColor.getrgb(color) + (100,)
elif isinstance(color, tuple):
fill_color = color + (100,)
draw.rectangle(bbox, width=width, outline=color, fill=fill_color)
else:
draw.rectangle(bbox, width=width, outline=color)
if labels is not None:
draw.text((bbox[0], bbox[1]), labels[i], fill=color, font=txt_font)
return torch.from_numpy(np.array(img_to_draw)).permute(2, 0, 1).to(dtype=torch.uint8)
def draw_segmentation_masks(
image: torch.Tensor,
masks: torch.Tensor,
alpha: float = 0.2,
colors: Optional[List[Union[str, Tuple[int, int, int]]]] = None,
) -> torch.Tensor:
if not isinstance(image, torch.Tensor):
raise TypeError(f"Tensor expected, got {type(image)}")
elif image.dtype != torch.uint8:
raise ValueError(f"Tensor uint8 expected, got {image.dtype}")
elif image.dim() != 3:
raise ValueError("Pass individual images, not batches")
elif image.size()[0] != 3:
raise ValueError("Pass an RGB image. Other Image formats are not supported")
num_masks = masks.size()[0]
masks = masks.argmax(0)
if colors is None:
palette = torch.tensor([2 ** 25 - 1, 2 ** 15 - 1, 2 ** 21 - 1])
colors_t = torch.as_tensor([i for i in range(num_masks)])[:, None] * palette
color_arr = (colors_t % 255).numpy().astype("uint8")
else:
color_list = []
for color in colors:
if isinstance(color, str):
# This will automatically raise Error if rgb cannot be parsed.
fill_color = ImageColor.getrgb(color)
color_list.append(fill_color)
elif isinstance(color, tuple):
color_list.append(color)
color_arr = np.array(color_list).astype("uint8")
_, h, w = image.size()
img_to_draw = Image.fromarray(masks.byte().cpu().numpy()).resize((w, h))
img_to_draw.putpalette(color_arr)
img_to_draw = torch.from_numpy(np.array(img_to_draw.convert('RGB')))
img_to_draw = img_to_draw.permute((2, 0, 1))
return (image.float() * alpha + img_to_draw.float() * (1.0 - alpha)).to(dtype=torch.uint8)
