Python实现平行束滤波反投影——Inverse Radon Transformation

参考博文进行了平行束滤波反投影的修改，将时域滤波修改为频域滤波，重建后消除原博文中图像的竖条状伪影。

https://blog.csdn.net/hsyxxyg/article/details/106433940

频域平行束滤波反投影（反radon变换）

产生频域滤波信号：
nextpow2函数的实现可参见github：

https://github.com/freenowill/Denoise-/blob/ba99babb685f6c80f07a23275c8c2127de601815/Denoise/nextpow2/nextpow2.py#L20

def Ramp(channNum,channSpacing): N = 2**nextpow2(2*channNum-1)ramp = np.zeros(N)for i1 in range(-(channNum-1),(channNum-1)):if i1==0:ramp[i1+channNum] = 1/(4*channSpacing*channSpacing)elif i1 % 2 ==0:ramp[i1+channNum] = 0elif i1 % 2 ==1:ramp[i1+channNum] = -1/((i1*np.pi*channSpacing)**2)# ramp = channSpacing*np.abs(np.fft.fft(ramp))ramp = channSpacing*np.abs(fft(ramp))# ramp = fftshift(ramp)return ramp

进行滤波反投影：

def IRandonTransform2(image, steps):AglPerView = np.pi/stepschannels = 512origin = np.zeros((steps, channels, channels))step = list(np.arange(0,1,1/100))step2 = step.copy()step2.reverse()# 防止truncation，在投影通道方向进行过渡# 如果sinogram通道方向的边界为0可以不进行过渡，否则重建图像的边缘可能会有一个亮环step_temp = image[0,:] # (360,)step_temp = np.expand_dims(step_temp,axis=0) # (1,360)step_temp = step_temp.repeat(100,axis=0) # (100,360)step = np.array(step) # (100,)step = np.expand_dims(step,axis=1) # (100,1)step = step.repeat(360,axis=1) # (100,360)step = step*step_temp # (100,360)step2_temp = image[-1,:] # (360,)step2_temp = np.expand_dims(step2_temp,axis=0) # (1,360)step2_temp = step2_temp.repeat(100,axis=0) # (100,360)step2 = np.array(step2) # (100,)step2 = np.expand_dims(step2,axis=1) # (100,1)step2 = step2.repeat(360,axis=1) # (100,360)step2 = step2*step2_temp # (100,360)    proj = np.concatenate((step,image,step2),axis=0) # (712,360)proj = np.concatenate((proj,np.zeros((2048-712,360)))) # (2048,360)sino_fft = fft(proj,axis=0) # (2048,360)'''卷积滤波(频域)'''filterData = Ramp(2*len(step)+512,0.5) # (2048,)iLen = len(filterData)filterData = np.expand_dims(filterData,axis=1) #(2048,1)filterData = filterData.repeat(360,axis=1) # (2048,360)image_filter = filterData*sino_fft # (2048,360)image_filter_ = ifft(image_filter,axis=0) # (2048,360)image_filter2 = np.real(image_filter_) # (2048,360)image_filter_final = image_filter2[100:512+100,:] # (512,360)image_filter3 = np.fliplr(image_filter_final) # (512,360)for i in range(steps): # viewNumprojectionValueFiltered = image_filter3[:,i]projectionValueExpandDim = np.expand_dims(projectionValueFiltered, axis=0) # 1*512projectionValueRepeat = projectionValueExpandDim.repeat(channels, axis=0) # 512*512origin[i] = ndimage.rotate(projectionValueRepeat, 180+i*180/steps, reshape=False).astype(np.float64)iradon = np.sum(origin, axis=0)iradon = iradon*AglPerViewreturn iradon，image_filter_final

主函数中用到的读图小函数：

def readdat(file,shape):fid = open(file,'rb')nefid = fid.read()fid.close()neimage = np.reshape(np.frombuffer(nefid,dtype=np.float32),shape,'F')return  neimage def readdicom(filename):img = pydicom.dcmread(filename).pixel_array.astype('int64')-24return img

主函数：

import numpy as np
from scipy import ndimage
from scipy.signal import convolve
import matplotlib.pyplot as plt
import pydicom
from scipy.fftpack import fft,ifft,fftshift,ifftshift
import ctypes as ct
import scipy.io as scio
import matplotlib.image as mpimgori_image = readdicom(r'E:\Metal Artifact Reduction\DLrelatedMAC_512_360_ratio1_noNorm\7_SimSource_img\noMetal\1_output_noM.dcm')
proj = readdat(r'E:\Metal Artifact Reduction\DLrelatedMAC_512_360_ratio1_noNorm\7_SimSource_img\Label_Train\1_label_train.dat',[512,360])
iradon,image_filter_final = IRandonTransform2(proj, 360)
plt.subplot(2, 2, 1)
plt.imshow(ori_image, cmap='gray')
plt.title('original image')
plt.subplot(2, 2, 2)
plt.imshow(iradon, cmap='gray')
plt.title('iradon image')
plt.subplot(2,2,3)
plt.imshow(proj,'gray')
plt.title('original projection')
plt.subplot(2,2,4)
plt.imshow(image_filter_final,'gray')
plt.title('filtered projection')
plt.show()

滤波反投影结果

Python实现平行束滤波反投影——Inverse Radon Transformation相关推荐

利用python进行平行束FBP重建结果
所涉及的知识点: FBP重建的步骤: 1.正投影 2.对每一个角度的投影数据补零,是至少补N-1个零:N是一个角度投影数据的长度. 3.构建和补零后数据长度一样的空间滤波器. 4.分别将空间滤波器和每 ...
pythonsl火车加字_荐Python实现Radon变换——直接反投影和滤波反投影
前几天我学习了Radon变换并用Python做了一个简单的程序(见上一篇博文),昨天看了一下逆Radon变换,尝试了一下简单的实现.我们可以通过对Sinogram图使用逆Radon变换来还原原始图像, ...
对投影值进行线性插值之后再进行滤波反投影的Python实现
前面一篇文章中我介绍了滤波反投影,实际中我们的扫描都是分立而非连续的,因此我们通常需要对投影值进行插值之后再进行滤波反投影,这样能够获得更好的效果.我现在先把代码贴上来,具体的数学过程过几天再详细讲. ...
【转】由投影重建图像：滤波反投影、FDK、TFDK三维重建算法理论基础
转自:由投影重建图像:滤波反投影.FDK.TFDK三维重建算法理论基础_m0_37357063的博客-CSDN博客_fdk算法 1. 基础理论从: [1] RafaelC.Gonzalez, Rich ...
【转】CT图像重构方法详解——傅里叶逆变换法、直接反投影法、滤波反投影法
转自:CT图像重构方法详解--傅里叶逆变换法.直接反投影法.滤波反投影法_Absolute Zero-CSDN博客_反投影法绪在做CT图像处理的时候遇到很多问题,对于滤波反变换有许多细 ...
【youcans 的 OpenCV 例程 200 篇】112. 滤波反投影重建图像
欢迎关注『youcans 的 OpenCV 例程 200 篇』系列,持续更新中欢迎关注『youcans 的 OpenCV学习课』系列,持续更新中 [youcans 的 OpenCV 例程 2 ...
matlab fbp fan arc,滤波反投影重建算法（FBP）实现及应用（matlab）
滤波反投影重建算法实现及应用(matlab) 1. 滤波反投影重建算法原理滤波反投影重建算法常用在CT成像重建中,背后的数学原理是傅立叶变换:对投影的一维傅立叶变换等效于对原图像进行二维的傅立叶变换 ...
CT图像重构方法详解——傅里叶逆变换法、直接反投影法、滤波反投影法
绪在做CT图像处理的时候遇到很多问题,对于滤波反变换有许多细节存在疑问,经过多天查找资料和利用MATLAB程序一步步实现后终于豁然开朗,于是想要总结成文,作为笔记方便今后查看.文中若有错误欢迎指出! ...
滤波反投影重建算法（FBP）实现及应用（matlab）
滤波反投影重建算法实现及应用(matlab) 1. 滤波反投影重建算法原理滤波反投影重建算法常用在CT成像重建中,背后的数学原理是傅立叶变换:对投影的一维傅立叶变换等效于对原图像进行二维的傅立叶变换 ...

Python实现平行束滤波反投影——Inverse Radon Transformation

频域平行束滤波反投影（反radon变换）

滤波反投影结果

Python实现平行束滤波反投影——Inverse Radon Transformation相关推荐

最新文章

热门文章