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

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

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] = 0
        elif 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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

进行滤波反投影：

def IRandonTransform2(image, steps):
    AglPerView = np.pi/steps
    channels = 512
    origin = 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): # viewNum
        projectionValueFiltered = image_filter3[:,i]
        projectionValueExpandDim = np.expand_dims(projectionValueFiltered, axis=0) # 1*512
        projectionValueRepeat = projectionValueExpandDim.repeat(channels, axis=0) # 512*512
        origin[i] = ndimage.rotate(projectionValueRepeat, 180+i*180/steps, reshape=False).astype(np.float64)
    iradon = np.sum(origin, axis=0)
    iradon = iradon*AglPerView

    return iradon，image_filter_final 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

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

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')-24
    return img
1
2
3
4
5
6
7
8
9
10

主函数：

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 mpimg

ori_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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

滤波反投影结果