Pytorch实现CT图像正投影（FP）与反投影（FBP）的模块

FP/FBP Modules

有关CT图像重建或图像处理的训练任务有时需要数据在投影域和图像域上进行变换，为了能使梯度在投影域和图像域之间进行传播，需要实现Forward Projection与Back Projection模块。
参考文献中提到的平行束正反投方法，可以进行网络中正反投模块的Pytorch实现。

[1] Zhao, J. , et al. “Unsupervised Learnable Sinogram Inpainting Network (SIN) for Limited Angle CT reconstruction.” (2018).

正反投模块的原理示意图下图1所示：radon变换即正投的步骤可以看成是旋转+累加，radon反变换即反投的过程可以看成是滤波+旋转（反投）+累加的过程。radon反变换的滤波可以在时域进行也可以在频域进行，因为时域的卷积等于频域的乘积，本博文是在频率上进行得Ramp滤波。

图1 Radon变换与反变换示意图

Pytorch实现

正投模块（FP）的实现：

class FP(nn.Module):def __init__(self, viewNum, chanNum, batchSize):super(FP, self).__init__()self.viewNum = viewNumself.chanNum = chanNumself.batchSize = batchSizedef forward(self, x):'''x: image x is a tensor (batchSize*netChanNum*imgSize*imgSize)'''sino = torch.from_numpy( np.zeros((self.batchSize, 1, self.chanNum, self.viewNum))).type(torch.FloatTensor) # batchSize*channel*512*360sino = sino.cuda()''' rotate'''for i in range(self.viewNum):angle = - 180/self.viewNum*(i+1) * math.pi / 180 - math.piA = np.array([[np.cos(angle), -np.sin(angle)],[np.sin(angle), np.cos(angle)]])  theta = np.array([[A[0, 0], A[0, 1], 0], [A[1, 0], A[1, 1], 0]])                                    theta = torch.from_numpy(theta).type(torch.FloatTensor)theta = theta.unsqueeze(0)theta = theta.repeat(self.batchSize,1,1)theta = theta.cuda()''' interpolation'''grid = F.affine_grid(theta, x.size())x_rotate = F.grid_sample(x, grid) # 4*1*512*512''' accumulation'''sino[:,:,:,i] = torch.sum(x_rotate, dim=2) sino = sino*0.5sino = sino.cuda()    return sino

反投模块（FBP）的实现：

class FBP(nn.Module):def __init__(self, viewNum, chanNum, batchSize, netChanNum, chanSpacing):super(FBP, self).__init__()self.viewNum = viewNum # projection的投影角度数self.chanNum = chanNum # projection的通道数self.batchSize = batchSizeself.netChanNum = netChanNum # 输入FBP网络数据的通道数self.chanSpacing = chanSpacingdef forward(self, x):'''x:  projection (batchSize*netChanNum*chanNum*viewNum) 4*1*512*360type(x) is a tensor''''''频域滤波'''projectionValue = convolution(x,self.batchSize,self.netChanNum,self.chanNum,self.viewNum,self.chanSpacing) # 2*1*512*360projectionValue = projectionValue.cuda()sino_rotate = np.zeros((self.batchSize, self.netChanNum, self.viewNum, self.chanNum, self.chanNum)) # batchSize*netChanNum*viewNum*chanNum*chanNumsino_rotate = torch.from_numpy(sino_rotate).type(torch.FloatTensor)sino_rotate = sino_rotate.cuda()AglPerView = math.pi/self.viewNum'''设置FOV,生成mask将FOV以外的区域置零'''FOV = torch.ones((self.batchSize,self.netChanNum,self.chanNum,self.chanNum))x_linespace = np.arange(1,self.chanNum+1,1)  # (512,)y_linespace = np.arange(1,self.chanNum+1,1)  # (512,)x_mesh,y_mesh = np.meshgrid(x_linespace,y_linespace) # 512*512XPos = (x_mesh-256.5) * self.chanSpacing # 512*512YPos = (y_mesh-256.5) * self.chanSpacing # 512*512R = np.sqrt(XPos**2 + YPos**2) # 512*512R = torch.from_numpy(R).type(torch.FloatTensor) # 512*512R = R.repeat(self.batchSize,self.netChanNum,1,1) # 2*1*512*512FOV[R>=self.chanSpacing*self.chanNum/2] = 0 # 2*1*512*512FOV = FOV.cuda()''' rotate interpolation'''for i in range(self.viewNum):projectionValueFiltered = torch.unsqueeze(projectionValue[:,:,:,i],3) # 2*1*512*1projectionValueRepeat = projectionValueFiltered.repeat(1,1,1,512) # 2*1*512*512projectionValueRepeat = projectionValueRepeat * FOV  # 2*1*512*512angle = -math.pi/2 + 180/self.viewNum*(i+1) * math.pi / 180A = np.array([[np.cos(angle), -np.sin(angle)],[np.sin(angle), np.cos(angle)]])theta = np.array([[A[0, 0], A[0, 1], 0], [A[1, 0], A[1, 1], 0]])theta = torch.from_numpy(theta).type(torch.FloatTensor)theta = theta.unsqueeze(0)theta = theta.cuda()theta = theta.repeat(self.batchSize,1,1)grid = F.affine_grid(theta, torch.Size((self.batchSize, self.netChanNum, 512, 512)))sino_rotate[:,:,i,:,:] = F.grid_sample(projectionValueRepeat, grid) # 512*512''' accumulation'''iradon = torch.sum(sino_rotate, dim=2)  iradon = iradon*AglPerViewreturn iradon

频域滤波的实现：其中调用到的Ramp()函数参照博文

https://blog.csdn.net/kouwang9779/article/details/115961582

def convolution(proj,batchSize,netChann,channNum,viewnum,channSpacing):AglPerView = np.pi/viewnumchannels = 512origin = np.zeros((batchSize,netChann,viewnum, channels, channels))# avoid truncationstep = list(np.arange(0,1,1/100))step2 = step.copy()step2.reverse()step = np.array(step) # (100,)step = np.expand_dims(step,axis=1) # 100*1step = torch.from_numpy(step).type(torch.FloatTensor) # (100,1)step = step.repeat(batchSize,1,1,360) # 2*1*100*360step_temp = proj[:,:,0,:].unsqueeze(2) # 2*1*1*360step_temp = step_temp.repeat(1,1,100,1) # 2*1*100*360step = step.cuda()step = step*step_temp # 2*1*100*360step2 = np.array(step2) # (100,)step2 = np.expand_dims(step2,axis=1) # 100*1step2 = torch.from_numpy(step2).type(torch.FloatTensor) # (100,1)step2 = step2.repeat(batchSize,1,1,360) # 2*1*100*360step2_temp = proj[:,:,-1,:].unsqueeze(2) # 2*1*1*360step2_temp = step2_temp.repeat(1,1,100,1) # 2*1*100*360step2 = step2.cuda()step2 = step2*step2_temp # 2*1*100*360filterData = Ramp(batchSize,netChann,2*100+channNum,channSpacing) # 2*1*2048*360iLen = filterData.shape[2] # 2048proj = torch.cat((step,proj,step2),2) # 2*1*712*360proj = torch.cat((proj,torch.zeros(batchSize,netChann,iLen-proj.shape[2],viewnum).cuda()),2) # 2*1*2048*360sino_fft = fft(proj.detach().cpu().numpy(),axis=2) # 2*1*2048*360image_filter = filterData*sino_fft # 2*1*2048*360image_filter_ = ifft(image_filter,axis=2) # 2*1*2048*360image_filter_ = np.real(image_filter_)image_filter_ = torch.from_numpy(image_filter_).type(torch.FloatTensor)image_filter_final = image_filter_[:,:,100:512+100] # 2*1*512*360return image_filter_final

结果展示

图2 原图像（左），通过FP模块后的正投结果（中），正投结果通过FBP模块后的反投结果（右）

讨论

以这种方法进行正反投之后的误差即原图像与反投结果之差如图3所示，这个误差范围应该可以通过网络的训练进行弥补。

图3 原图像与反投影图像的误差图

Pytorch实现CT图像正投影（FP）与反投影（FBP）的模块相关推荐

【转】CT图像重构方法详解——傅里叶逆变换法、直接反投影法、滤波反投影法
转自:CT图像重构方法详解--傅里叶逆变换法.直接反投影法.滤波反投影法_Absolute Zero-CSDN博客_反投影法绪在做CT图像处理的时候遇到很多问题,对于滤波反变换有许多细 ...
CT图像重构方法详解——傅里叶逆变换法、直接反投影法、滤波反投影法
绪在做CT图像处理的时候遇到很多问题,对于滤波反变换有许多细节存在疑问,经过多天查找资料和利用MATLAB程序一步步实现后终于豁然开朗,于是想要总结成文,作为笔记方便今后查看.文中若有错误欢迎指出! ...
研究 | CT图像迭代重建算法研究进展
上次讲到我实现了一下代数迭代重建(ART),到周六开会的时候才大概了解了我的研究方向应该是统计迭代重建,一下子就把我给搞懵了.按照书上的说法,统计迭代法是在发射型CT(SPECT和PET)中应用广泛, ...
acdsee扫描没有图像_详解CT图像常见伪影成因及解决方法
CT 图像伪影 (artifact) 是指重建图像上与实际解剖结构不相符的密度异常变化.CT 图像比传统平扫X线更容易出现伪影,这是因为 CT 图像是由成千上万独立的原始测量数据重建而得,而计算机重建 ...
matlab 画图断层显示,MATLAB编程实现连续断层工业CT图像的三维重建_张爱东
第26卷第4期核电子学与探测技术 V ol.26 N o.4 2006年 7月Nuclear Electr onics &Detection T echnolo gy July 2006 M ...
CT图像常见伪影及解决方法
- 前言 - CT的伪影理论上可被定义为图像中被重建数值与物体真实衰减系数之间的差异,简单来说,对于图像重建过程中不该出现在图像上的影像,可认为其是伪影(antifacts). - 01 伪影的分类 ...
ct图像中的金属伪影
ct图像中的金属伪影校正方法背景随着现代医学的发展,通过守住在病人体内植入带有金属物质的假体的情况越来越普遍.最常见的是假牙植入.心脏起搏器,以及越来越多的各种关节和假肢等.这些金属物体相比较于人 ...
【深度学习】一个应用—肝脏CT图像自动分割（术前评估）
[深度学习]一个应用-肝脏CT图像自动分割(术前评估) 文章目录 1 目标 2 数据集 3 LITS20173.1 LiTS数据的预处理3.2 LiTS数据的读取3.3 数据增强3.4 数据存储 4 ...
【项目实战课】基于Pytorch的SRGAN图像超分辨实战
欢迎大家来到我们的项目实战课,本期内容是<基于Pytorch的SRGAN图像超分辨实战>.所谓项目实战课,就是以简单的原理回顾+详细的项目实战的模式,针对具体的某一个主题,进行代码级的实战 ...

Pytorch实现CT图像正投影（FP）与反投影（FBP）的模块

FP/FBP Modules

Pytorch实现

结果展示

讨论

Pytorch实现CT图像正投影（FP）与反投影（FBP）的模块相关推荐

最新文章

热门文章