標題

• • 由投影重建圖像
  • • 投影和雷登變換 Johann Radon
    • 反投影
    • 濾波反投影重建

由投影重建圖像

本由投影重建圖像，主要是雷登變換與雷登把變換的應用，所以也沒有太多的研究，只為了保持完整性，而添加到這里。

# 自制旋轉投影圖像# 模擬一個圖像 img = np.zeros((100, 100)) img[45:46] = 0.25 img[46:47] = 0.5 img[47:48] = 0.75 img[48:51, :] = 1. img[51:52, :] = 0.75 img[52:53, :] = 0.5 img[53:54, :] = 0.25total_img = img.copy()# 旋轉并疊加, it's 32 on the book, but here I want to try 64, but degree / 2, see how is work, seems all the same for i in range(64): angle = (5.625) / 2 * (i + 1)C1 = cv2.getRotationMatrix2D((img.shape[1]/2.0, img.shape[0]/2.0), angle, 1)new_img = cv2.warpAffine(img, C1, (img.shape), borderValue=0)total_img = total_img + new_imgtotal_img = normalize(total_img)plt.figure(figsize=(13, 5)) plt.subplot(131), plt.imshow(img, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([]) plt.subplot(132), plt.imshow(new_img, 'gray'), plt.title('new_img'), plt.xticks([]), plt.yticks([]) plt.subplot(133), plt.imshow(total_img, 'gray'), plt.title('total_img'), plt.xticks([]), plt.yticks([]) plt.tight_layout() plt.show()

投影和雷登變換 Johann Radon

$$g(\rho ,\theta )={\int }_{?\infty }^{\infty }{\int }_{?\infty }^{\infty }f(x,y)\delta (x\theta +ysin\theta ?\rho )dxdy$$

離散情況的雷登變換
$$g(\rho ,\theta )=\sum _{x=0}^{M?1}\sum _{y=0}^{N?1}f(x,y)\delta (x\theta +ysin\theta ?\rho )$$

由雷登變換得到的圖稱之為正弦圖(Sinogram)

from scipy import ndimagedef radon_transform(img, steps):"""radon tranform for gray imageparam: img: input Imageparam: steps: step for the transform, same of image height"""channels = img.shape[0]dst = np.zeros((channels, channels), dtype=np.float32)for i in range(steps):res = ndimage.rotate(img, -i * 180 / steps, reshape=False).astype(np.float32)dst[:, i] = sum(res)# dst = ndimage.rotate(dst, 270, reshape=False) # 旋轉后就跟書上的結果一致return dst # 雷登變換 img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0533(a)(circle).tif", 0) m, n = img_ori.shape[1], img_ori.shape[0] img = idea_low_pass_filter(img_ori, (m, n), D0=10) img_radon = radon_transform(img, img.shape[0])plt.figure(figsize=(24, 12))plt.subplot(131), plt.imshow(img, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([]) plt.subplot(132), plt.imshow(img_radon, 'gray'), plt.title('Randon Transform'), plt.xticks([]), plt.yticks([]) # plt.subplot(133), plt.imshow(img_cv2_guass, 'gray'), plt.title('CV2 Guass') plt.tight_layout() plt.show()

# 兩個物體的雷登變換 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0534(a)(ellipse_and_circle).tif', 0) #直接讀為灰度圖像 img_radon = radon_transform(img_ori, img_ori.shape[0])plt.figure(figsize=(24, 12)) plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([]) plt.subplot(132), plt.imshow(img_radon, 'gray'), plt.title('Randon Transform'), plt.xticks([]), plt.yticks([]) # plt.subplot(133), plt.imshow(img_cv2_guass, 'gray'), plt.title('CV2 Guass') plt.tight_layout() plt.show()

# 長方形的雷登變換 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0539(a)(vertical_rectangle).tif', 0) #直接讀為灰度圖像img_radon = radon_transform(img_ori, img_ori.shape[0])plt.figure(figsize=(24, 12))plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([]) plt.subplot(132), plt.imshow(img_radon, 'gray'), plt.title('Randon Transform'), plt.xticks([]), plt.yticks([]) # plt.subplot(133), plt.imshow(img_cv2_guass, 'gray'), plt.title('CV2 Guass') plt.tight_layout() plt.show()

# 多個物體雷登變換 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0539(c)(shepp-logan_phantom).tif', 0) #直接讀為灰度圖像img_radon = radon_transform(img_ori, img_ori.shape[0])plt.figure(figsize=(24, 12))plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([]) plt.subplot(132), plt.imshow(img_radon, 'gray'), plt.title('Randon Transform'), plt.xticks([]), plt.yticks([]) # plt.subplot(133), plt.imshow(img_cv2_guass, 'gray'), plt.title('CV2 Guass') plt.tight_layout() plt.show()

反投影

$$f_{\theta }(x,y)=g(xcos\theta +ysin\theta ,\theta )$$
$$f(x,y)={\int }_0^{\pi }{f}_{\theta }(x,y)d\theta$$
$$f(x,y)=\sum _{\theta =0}^{\pi }{f}_{\theta }(x,y)$$

反投影的圖像有時稱為層圖，我們可把層圖理解為一幅由其生成投影的圖像的一個近似。

from scipy import ndimagedef inverse_radon_transform(image, steps):"""inverse radon tranform for radon transform imageparam: image: input Imageparam: steps: step for the transform, normaly same of image heightreturn: inverse radon transform for image, image is un-normalized"""channels = len(image[0])dst = np.zeros((steps, channels, channels))for i in range(steps):# 傳入的圖像中的每一列都對應于一個角度的投影值# 這里用的圖像是上篇博文里得到的Radon變換后的圖像裁剪后得到的temp = image[:, i]# 這里利用維度擴展和重復投影值數組來模擬反向均勻回抹過程temp_expand_dim = np.expand_dims(temp, axis=0)temp_repeat = temp_expand_dim.repeat(channels, axis=0)dst[i] = ndimage.rotate(temp_repeat, i*180 / steps, reshape=False).astype(np.float64)# 各個投影角度的投影值已經都保存在origin數組中，只需要將它們相加即可 iradon = np.sum(dst, axis=0)return iradon # 雷登反變換 # img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0539(c)(shepp-logan_phantom).tif', 0) #直接讀為灰度圖像# img_radon = radon_transform(img_ori, img_ori.shape[0])img_inverse_radon = inverse_radon_transform(img_radon, img_radon.shape[0])plt.figure(figsize=(24, 12))plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original') plt.subplot(132), plt.imshow(img_radon, 'gray'), plt.title('Randon Transform') plt.subplot(133), plt.imshow(img_inverse_radon, 'gray'), plt.title('Inverse Randon Transform') plt.tight_layout() plt.show()

# 兩個物體的雷登反變換 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0534(a)(ellipse_and_circle).tif', 0) #直接讀為灰度圖像img_radon = radon_transform(img_ori, img_ori.shape[0]) img_inverse_radon = inverse_radon_transform(img_radon, img_radon.shape[0])plt.figure(figsize=(24, 12)) plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original') plt.subplot(132), plt.imshow(img_radon, 'gray'), plt.title('Randon Transform') plt.subplot(133), plt.imshow(img_inverse_radon, 'gray'), plt.title('Inverse Randon Transform') plt.tight_layout() plt.show()

濾波反投影重建

$$f(x,y)={\int }_0^{\pi }{\int }_{?\infty }^{\infty }|\omega |G(\omega ,\theta )e^{j2\pi \omega (xcos\theta +ysin\theta )}d\omega d\theta$$

$$f(x,y)={\int }_0^{\pi }[{\int }_{?\infty }^{\infty }|\omega |G(\omega ,\theta )e^{j2\pi \omega \rho }d\omega {]}_{\rho =xcos\theta +ysin\theta }d\theta$$

$$h(\omega )=\{\begin{array}{ll}c+(c?1)cos\frac{2\pi \omega }{M},& 0\le \omega \le (M?1)\\ 0,& others\end{array}$$

當$c=0.54$時，該函數稱為漢明窗(Richard Hamming)；當$c=0.5$時，該函數稱為韓窗(Julius von Hann)。漢明窗和韓窗的主要區別是，韓窗末尾的一些點為零。通常兩者的差別在圖像處理應用中是覺察不到的。

反投影圖像$f(x,y)$是按如下步驟得到的：

計算每個投影的一維傅里葉變換。

將每個全里葉變換乘以濾波器傳弟函數$|\omega |$，這個傳遞函數已乘以一個合適的窗（如漢明窗）。

得到的每個濾波后的變換的一維傅里葉反變換。

對步驟3得到的所有一維反變換進行積分（求和）。

顧名思義，濾波反投影就是先對投影值進行濾波，然后利用的得到的值進行反投影，簡單來說濾波的主要目的就是使邊緣（對應于圖像的高頻部分）更加明顯。理論上來說，濾波函數應該是：
$$h(\omega )=|\omega |$$

但是這個是一個理想濾波器，沒辦法實現，因此需要考慮其他能夠實現并且能夠使得到的圖像更加接近原始圖像的濾波器。這里我僅介紹兩種R—L濾波器和S—L濾波器，下面是這兩種濾波器的具體形式：
$$h_{RL}(n\delta )=?\begin{array}{ll}?\frac{1}{4{\delta }^2},& n=0\\ 0,& n為偶數\\ ?\frac{1}{(n\pi \delta {)}^2},& n為奇數\end{array}$$

$$h_{SL}=\frac{1}{{\pi }^2{\delta }^2(4n^2?1)}$$

利用以上兩個濾波器和投影值進行卷積，然后再進行反投影，就可以得到得到原始的圖像，大致上來說，就是在上面的代碼中增加濾波器的構建和投影與濾波器的卷積過程，具體的代碼如下：

def box_filter(M, c):hw = np.zeros((M, ))for w in range(M):if 0 <= w <= M - 1:hw[w] = c + (c - 1) * np.cos(2 * np.pi * w / M)else:hw[w] = 0hw = np.fft.fft(hw)hw = np.fft.fftshift(hw)hw = np.real(hw) # hw = normalize(np.real(hw))return hw # 1維漢明窗 x = np.linspace(-5, 5, 200)M = x.shape[0]y = abs(x)c = 0.54hw = np.zeros_like(x) for w in range(M):if 0 <= w <= M - 1:hw[w] = c + (c - 1) * np.cos(2 * np.pi * w / M)else:hw[w] = 0np_hamming = np.hamming(M) np_hann= np.hanning(M)# plt.plot(x, y) plt.plot(x, hw) plt.plot(x, np_hamming) plt.plot(x, np_hann) plt.show()

def hamming_window(img, hamming=True):"""2d hamming and hann windowparam: img, input imageparam: hamming: bool, if True return hamming windows, if False return hann windowreturn normalize 2d hamming or hann windowProblem: I still not very sure if this right, since result is not very good."""M, N = img.shape[1], img.shape[0]if hamming:u = np.hamming(M)v = np.hamming(N)else:u = np.hanning(M)v = np.hanning(N)u, v = np.meshgrid(u, v)high_pass = np.sqrt(u**2 + v**2)# high_pass = np.sqrt((u - M//2)**2 + (v - N//2)**2)kernel = high_pass # 1 / (1 + ((high_pass * W) / (high_pass ** 2 - D0 **2 + 1e-5))**(2*order))return kernel # 二維漢明窗 img_hamming = hamming_window(img_radon, hamming=True)# img_radon_hamming = normalize(img_hamming * img_radon) img_radon_hamming = img_hamming * img_radonplt.figure(figsize=(24, 12))plt.subplot(131), plt.imshow(img_radon, 'gray'), plt.title('img radon'), plt.xticks([]), plt.yticks([]) plt.subplot(132), plt.imshow(img_hamming, 'gray'), plt.title('img_hamming'), plt.xticks([]), plt.yticks([]) plt.subplot(133), plt.imshow(img_radon_hamming, 'gray'), plt.title('img_radon_hamming'), plt.xticks([]), plt.yticks([]) plt.tight_layout() plt.show()

#兩種濾波器的實現 def rl_filter(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) + 1e-5) # 1e-5 加上一個不為零的小數，防止出現除0的問題if np.mod(i - N / 2, 2) == 0:filterRL[i] = 0filterRL[int(N/2)] = 1 / (4 * np.power(d, 2.0))return filterRLdef sl_filter(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 inverse_filter_radon_transform(image, steps):#定義用于存儲重建后的圖像的數組channels = len(image[0])origin = np.zeros((steps, channels, channels)) # filter = box_filtechannelschannels, 0.48) # filter = rl_filter(channels, 1)filter = sl_filter(channels, 1)for i in range(steps):projectionValue = image[:, i]projectionValueFiltered = np.convolve(filter, projectionValue, "same")projectionValueExpandDim = np.expand_dims(projectionValueFiltered, axis=0)projectionValueRepeat = projectionValueExpandDim.repeat(channels, axis=0)origin[i] = ndimage.rotate(projectionValueRepeat, i*180/steps, reshape=False).astype(np.float64)iradon = np.sum(origin, axis=0)return iradon # 兩種濾波器的實現 hx = box_filter(200, 0.5) sl = sl_filter(200, 1)plt.figure(figsize=(10, 4)) plt.subplot(1, 2, 1), plt.plot(x, hx) plt.subplot(1, 2, 2), plt.plot(x, sl) plt.show()

# 多個物體的反投影 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0539(c)(shepp-logan_phantom).tif', 0) #直接讀為灰度圖像img_radon = radon_transform(img_ori, img_ori.shape[0]) img_inverse_radon = inverse_radon_transform(img_radon, img_radon.shape[0]) img_filter_inverse_radon = inverse_filter_radon_transform(img_radon, img_radon.shape[0])plt.figure(figsize=(18, 6))plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([]) plt.subplot(132), plt.imshow(img_inverse_radon, 'gray'), plt.title('Inverse Randon Transform'), plt.xticks([]), plt.yticks([]) plt.subplot(133), plt.imshow(img_filter_inverse_radon, 'gray'), plt.title('Inverse Filter Randon Transform') plt.xticks([]), plt.yticks([]) plt.tight_layout() plt.show()

# 長方形的反投影 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0539(a)(vertical_rectangle).tif', 0) #直接讀為灰度圖像img_radon = radon_transform(img_ori, img_ori.shape[0]) img_inverse_radon = inverse_radon_transform(img_radon, img_radon.shape[0]) img_hamming = hamming_window(img_radon, hamming=True) img_radon_hamming = img_hamming * img_radonimg_filter_inverse_radon = inverse_filter_radon_transform(img_radon_hamming, img_radon_hamming.shape[0])plt.figure(figsize=(18, 6))plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([]) plt.subplot(132), plt.imshow(img_inverse_radon, 'gray'), plt.title('Inverse Randon Transform'), plt.xticks([]), plt.yticks([]) plt.subplot(133), plt.imshow(img_filter_inverse_radon, 'gray'), plt.title('Inverse Filter Randon Transform') plt.xticks([]), plt.yticks([]) plt.tight_layout() plt.show()

總結

以上是生活随笔為你收集整理的第5章 Python 数字图像处理(DIP) - 图像复原与重建17 - 由投影重建图像、雷登变换、投影、反投影、反投影重建的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。