网站开发介绍ppt,网站数据分析指标,苏州网页开发公司,wordpress柒比贰主题医学图像处理 利用pytorch实现的可用于反传的Radon变换和逆变换 前言代码实现思路实验结果 前言
Computed Tomography#xff08;CT#xff0c;计算机断层成像#xff09;技术作为如今医学中重要的辅助诊断手段#xff0c;也是医学图像研究的重要主题。如今#xff0c;随… 医学图像处理 利用pytorch实现的可用于反传的Radon变换和逆变换 前言代码实现思路实验结果 前言
Computed TomographyCT计算机断层成像技术作为如今医学中重要的辅助诊断手段也是医学图像研究的重要主题。如今随着深度学习对CT重建、CT分割的研究逐渐深入论文开始越来越多的利用CT的每一个环节来扩充Feature或构造损失函数。
但是每到这个时候一个问题就出现了如果我要构造损失函数我势必要保证这个运算中梯度不会断掉否则起不到优化效果。而Radon变换目前好像没有人直接用其当作损失函数的一部分很奇怪故在此实现pytorch版本的Radon变换已经验证可以反传但是反传的对不对不敢保证只保证能计算出反传的值。希望能帮到需要的同学。
参考了这两篇博文在此十分感谢这位前辈。 Python实现离散Radon变换 Python实现逆Radon变换——直接反投影和滤波反投影
代码
from typing import Optional
import numpy as np
import matplotlib.pyplot as plt
import math
import torch as th
import torch.nn as nn
import torch.nn.functional as F
import SimpleITK as sitk#两种滤波器的实现
def RLFilter(N, d):filterRL np.zeros((N,))for i in range(N):filterRL[i] - 1.0 / np.power((i - N / 2) * np.pi * d, 2.0)if np.mod(i - N / 2, 2) 0:filterRL[i] 0filterRL[int(N/2)] 1 / (4 * np.power(d, 2.0))return filterRLdef SLFilter(N, d):filterSL np.zeros((N,))for i in range(N):filterSL[i] - 2 / (np.pi**2.0 * d**2.0 * (4 * (i - N / 2)**2.0 - 1))return filterSLdef IRandonTransform(image:th.Tensor|np.array, steps:Optional[int]None):Inverse Radon Transform(with Filter, FBP)Parameters:image: (B, C, W, H)# 定义用于存储重建后的图像的数组channels image.shape[-1]B, C, W, H image.shapeif steps None:steps channelsorigin th.zeros((B, C, steps, channels, channels), dtypeth.float32)filter_kernal th.tensor(SLFilter(channels, 1)).unsqueeze(0).unsqueeze(0).float()Filter nn.Conv1d(C, C, (channels), paddingsame,biasFalse)Filter.weight nn.Parameter(filter_kernal) for i in range(steps):# 传入的图像中的每一列都对应于一个角度的投影值# 这里用的图像是上篇博文里得到的Radon变换后的图像裁剪后得到的projectionValue image[:, :, :, i]projectionValue Filter(projectionValue)# 这里利用维度扩展和重复投影值数组来模拟反向均匀回抹过程projectionValueExpandDim projectionValue.unsqueeze(2)projectionValueRepeat projectionValueExpandDim.repeat((1, 1, channels, 1))origin[:,:, i] rotate(projectionValueRepeat, (i / steps) * math.pi)#各个投影角度的投影值已经都保存在origin数组中只需要将它们相加即可iradon th.sum(origin, axis2)return iradondef rotate(image:th.Tensor, angle):Rotate the image in any angles(include negtive).angle should be pi 180B image.shape[0]transform_matrix th.tensor([[math.cos(angle),math.sin(-angle),0],[math.sin(angle),math.cos(angle),0]], dtypeth.float32).unsqueeze(0).repeat(B,1,1)grid F.affine_grid(transform_matrix, # 旋转变换矩阵image.shape).float() # 变换后的tensor的shape(与输入tensor相同)rotation F.grid_sample(image, # 输入tensorshape为[B,C,W,H]grid, # 上一步输出的gird,shape为[B,C,W,H]modenearest) # 一些图像填充方法这里我用的是最近邻return rotationdef DiscreteRadonTransform(image:th.Tensor|np.array, steps:Optional[int]None):Radon TransformParameters:image : (B, C, W, H)channels image.shape[-1] # img_sizeB, C, W, H image.shaperes th.zeros((B, channels, channels), dtypeth.float32)if steps None:steps channelsfor s in range(steps):angle (s / steps) * math.pirotation rotate(image, -angle)res[:, :,s] th.sum(rotation, dim2)return res.unsqueeze(1)if __name__ __main__:origin sitk.ReadImage(/hy-tmp/data/LIDC/CT/0001.nii.gz)t_origin sitk.GetArrayFromImage(origin)t_origin th.tensor(t_origin)t_origin t_origin[40].unsqueeze(0).unsqueeze(0)a nn.Parameter(th.ones_like(t_origin))t_origin t_origin * aret DiscreteRadonTransform(t_origin) # (B, 1, H, W)b th.ones_like(ret)lf nn.MSELoss()loss lf(b, ret)loss.backward() # pytorch不会报错并能返回gradrec IRandonTransform(ret)ret ret.squeeze(0)rec rec.squeeze(0)plt.imsave(test.png, (t_origin.squeeze(0).squeeze(0)).data.numpy(), cmapgray)plt.imsave(test2.png, (ret.squeeze(0)).data.numpy(), cmapgray)plt.imsave(test3.png, (rec.squeeze(0)).data.numpy(), cmapgray)
实现思路
这份代码实际上是参考前文提到的前辈的代码修改而来具体而言就是把各种numpy实现的地方修改为pytorch的对应实现其中pytorch没有对应的API来实现矩阵的Rotate因此还参考了网上其它人实现的手写旋转的pytorch版本。并将其写作Rotate函数在其它任务中也可以调用这里需要注意调用时矩阵需要是方阵否则会出现旋转后偏移中心的问题。
实验结果
原始图像 Radon变换的结果 重建结果 这些图像可以下载下来自己试试。
