相关文章推荐
朝气蓬勃的油条  ·  Sql ...·  7 月前    · 
淡定的大脸猫  ·  Visual FoxPro 术语 - ...·  1 年前    · 
冷静的灯泡  ·  asp.net - Can not ...·  1 年前    · 
兴奋的石榴  ·  Android ...·  1 年前    · 
坏坏的春卷  ·  mysql数据库对后缀为.dump的文件的导 ...·  1 年前    · 
Code  ›  第5章 Python 数字图像处理(DIP) - 图像复原与重建17 - 由投影重建图像、雷登变换、投影、反投影、反投影重建
transform 图像滤波 图像复原 数字图像处理
https://blog.csdn.net/jasneik/article/details/115099488
高大的小摩托
11 月前
total_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_img total_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 + y \; sin\theta -\rho) \; dxdy g ( ρ , θ ) = ∫ − ∞ ∞ ​ ∫ − ∞ ∞ ​ f ( x , y ) δ ( x θ + y s i n θ − ρ ) d x d y

离散情况的雷登变换
g(\rho,\theta) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y) \; \delta(x\;\theta + y \; sin\theta -\rho) g ( ρ , θ ) = x = 0 ∑ M − 1 ​ y = 0 ∑ N − 1 ​ f ( x , y ) δ ( x θ + y s i n θ − ρ )

由雷登变换得到的图称之为正弦图(Sinogram) 

from scipy import ndimage
def radon_transform(img, steps):
    radon tranform for gray image
    param: img: input Image
    param: 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(x,y)=θ=0∑π​fθ​(x,y)
 

反投影的图像有时称为层图，我们可把层图理解为一幅由其生成投影的图像的一个近似。




    

 

from scipy import ndimage
def inverse_radon_transform(image, steps):
    inverse radon tranform for radon transform image
    param: image: input Image
    param: steps: step for the transform,  normaly same of image height
    return: 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)=∫0π​∫−∞∞​∣ω∣G(ω,θ)ej2πω(xcosθ+ysinθ)dωdθ
 
         f(x,y) = \int_0^\pi \bigg[\int_{-\infty}^{\infty}|\omega| G(\omega,\theta)e^{j2\pi\omega\rho} d\omega\bigg]_{\rho = xcos\theta+ysin\theta} d\theta 
     f(x,y)=∫0π​[∫−∞∞​∣ω∣G(ω,θ)ej2πωρdω]ρ=xcosθ+ysinθ​dθ
 
         h(\omega) = \begin{cases}c + (c-1)cos\frac{2\pi \omega}{M}, & {0 \leq\omega \leq (M-1)} \\ 0, & others\end{cases} 
     h(ω)={c+(c−1)cosM2πω​,0,​0≤ω≤(M−1)others​
 
         h_{RL}(n\delta)=\begin{cases} -\frac{1}{4\delta^2}, & n =0 \\ 0, & n为偶数 \\ -\frac{1}{(n\pi\delta)^2}, & n 为奇数 \end{cases} 
     hRL​(nδ)=⎩⎪⎨⎪⎧​−4δ21​,0,−(nπδ)21​,​n=0n为偶数n为奇数​
 
     hSL​=π2δ2(4n2−1)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] = 0
    hw = 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.54
hw = 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] = 0
np_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 window
    param: img, input image
    param: hamming: bool, if True return hamming windows, if False return hann window
    return normalize 2d hamming or hann window
    Problem: 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_radon
plt.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] = 0
    filterRL[int(N/2)] = 1 / (4 * np.power(d, 2.0))
    return filterRL
def 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 filterSL
def 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_radon
img_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()
                    标题瑞利噪声爱尔兰（伽马）噪声指数噪声均匀噪声椒盐噪声周期噪声估计噪声参数只存在噪声的复原 - 空间滤波均值滤波器算术平均滤波器几何均值滤波器谐波平均滤波器反（逆）谐波平均滤波器统计排序滤波器中值、最大值、最小值、中点 滤波器修正阿尔法均值滤波器自适应滤波器自适应局部降噪滤波器自适中值滤波器使用频率域滤波降低周期噪声陷波滤波深入介绍最优陷波滤波线性位置不变退化估计退化函数采用观察法估计退化函数采用试验法估计退化函数采用建模法估计退化函数运动模糊函数OpenCV Motion Blur逆滤波逆滤波最小均方误差
				
​本博客将实现Python版本的双目三维重建系统，项目代码实现包含：`双目标定`，`立体校正（含消除畸变)`，`立体匹配`，`视差计算`和`深度距离计算/3D坐标计算` 的知识点。限于篇幅，本博客不会过多赘述算法原理，而是手把手教你，如果搭建一套属于自己的双目三维重建的系统。该系统包含：
支持双USB连接线的双目摄像头
支持单USB连接线的双目摄像头(左右摄像头被拼接在同一个视频中显示)
支持单目相机标定:mono_camera_calibration.py
支持双目相机标定:stereo_camera




    


				
关于投影的基础知识：
假设我们要用一束细细的，平行的X射线从左到右穿过（通过一个图像平面），这里我们假设物体吸收的射线束能量 比背景吸收的射线束能量多。我们利用放在放在另一端的X射线吸收检测器来检测射线通过这个物体的能量
所以当左侧有一排X射线（光束带），右侧有一排X射线检测器（检测器带/吸收剖面）的时候，就会形成一副函数图像（通过接收器的位置和应该接收点接收能量的离散值形成），我们的目的是为了复原X射线穿过物体的图像，但是只是知道这个函数有什么用呢？
反投影法：
我们需要基于上述函数的信息来重建一

				
参考博文进行了平行束滤波反投影的修改，将时域滤波修改为频域滤波，重建后消除原博文中图像的竖条状伪影。
https://blog.csdn.net/hsyxxyg/article/details/106433940
频域平行束滤波反投影（反radon变换）
产生频域滤波信号：
nextpow2函数的实现可参见github：
https://github.com/freenowill/Denoise-/blob/ba99babb685f6c80f07a23275c8c2127de601815/Denois

				
图像重建之由投影重建图像简述部分引言
最初接触由投影重建图像这块内容的时候是在车牌识别中。上图是在90°的投影下的结果。
下面我们开一个简单图像的特定角度下的投影 当我们收集到各个角度的投影后，并希望通过这些投影的图像重建原图像。上图的最后一张图片就是直接重建的图像。可见，其有非常明显的“晕环”现象。
雷登变换
雷登变换阐述了一幅图像与其在各个角度下投影的具体表示。 下面我们在matlab上展示一副

import numpy as np
import matplotlib.pyplot as plt
from skimage.transform import radon, iradon
# 生成一个圆形模板
N = 128
R = 40
x, y = np.indices((N, N))
mask = ((x-N//2)**2 + (y-N//2)**2) < R**2
# 生成雷登变换图
theta = np.linspace(0, 180, 180, endpoint=False)
sinogram = radon(mask.astype(float), theta=theta, circle=True)
# 显示正弦图
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))
ax1.imshow(mask, cmap='gray')
ax1.set_title('Original')
ax2.imshow(sinogram, cmap='gray', extent=(-90, 90, 0, sinogram.shape[0]), aspect='auto')
ax2.set_title('Sinogram')
plt.show()
# 重建图像
reconstruction = iradon(sinogram, theta=theta, circle=True)
reconstruction *= np.pi / (2. * len(theta))  # 标准化
error = np.sqrt(np.mean((mask - reconstruction) ** 2))
# 显示重建图像和误差
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))
ax1.imshow(mask, cmap='gray')
ax1.set_title('Original')
ax2.imshow(reconstruction, cmap='gray')
ax2.set_title('Reconstruction\nRMSE: %.3g' % error)
plt.show()
这个示例代码生成一个直径为80的圆形模板，并使用Radon函数生成180个角度下的正弦图像，然后使用iradon函数进行重建。你可以根据自己的需求修改代码中的参数。
 
推荐文章
朝气蓬勃的油条  ·  Sql Server中用变量名作为表名的解决方案_preparedstatement sql變量為表名怎麼處理-CSDN博客
7 月前
淡定的大脸猫  ·  Visual FoxPro 术语 - Open Database Connectivity (ODBC) | Microsoft Learn
1 年前
冷静的灯泡  ·  asp.net - Can not find system.web.http - Stack Overflow
1 年前
兴奋的石榴  ·  Android JPEG编解码(YUV,RGB,JPEG格式转换) - 简书
1 年前
坏坏的春卷  ·  mysql数据库对后缀为.dump的文件的导出和导入_daiqinge的博客-CSDN博客
1 年前
今天看啥   ·   Py中国   ·   codingpro   ·   藏经阁   ·   小百科   ·   link之家   ·   卧龙AI搜索
删除内容请联系邮箱 2879853325@qq.com
Code - 代码工具平台
© 2024 ~ 沪ICP备11025650号