基于自相关和Python生成CT空间相关噪声的问题求助

复现背景 #

尝试复现该论文中添加噪声模拟剂量降低章节的方法，对作者实现步骤的理解如下：

• 从体模数据测量CT噪声的频谱特性
• 基于上述数据计算噪声自相关函数
• 在自相关峰值周围取窗口，保存为卷积滤波器
• 将该滤波器作用于高斯白噪声，缩放至所需标准差，得到具有对应功率谱的噪声
  生成的空间相关噪声可添加到实际患者CT图像中，得到噪声频谱特性与体模扫描一致的CT图像。

原实现代码 #

#! pip install pydicom
import matplotlib.pyplot as plt
import pydicom
import pydicom.data
import numpy as np
from numpy.fft import fft, ifft
from numpy import zeros
from scipy import signal

base = ""
pass_dicom1 = "Catphan36A.dcm"  # Phantom noise data

filename = pydicom.data.data_manager.get_files(base, pass_dicom1)[0]
ds = pydicom.dcmread(filename)

print("# show CT of phantom")
plt.imshow(ds.pixel_array, cmap=plt.cm.bone)  
plt.show()

n=512 # get center 128x128 pixels of phantom  scan, i.e. the uniform noise
dataNoise= ds.pixel_array
dataNoise = dataNoise[int(n*(3/8)):int(n*(3/8))+int(n/4), int(n*(3/8)):int(n*(3/8))+int(n/4)]

print("Show 12x128 uniform noise from Phantom")
plt.imshow(dataNoise, cmap="gray")  # show 12x128 uniform noise from Phantom
plt.show()

# do 2d DT of the phantom noise 
dataNoiseFT = np.fft.ifft2(dataNoise)

# compute the autocorrelation function of the phantom noise and shift the data to center to obtain kernel
dataAC = np.fft.ifft2(dataNoiseFT*np.conjugate(dataNoiseFT))
shiftedAC = np.fft.fftshift(dataAC)

print("Show 128x128 autocorrelation kernel")
plt.imshow(abs(shiftedAC), cmap="gray")  # show 128x128 kernel 
plt.show()

print("Show 32x32 autocorrelation kernel")
n = 128 # downsize kernel to 32x32
extractedAC = abs(shiftedAC)[int(n*(3/8)):int(n*(3/8))+int(n/4), int(n*(3/8)):int(n*(3/8))+int(n/4)]
extractedAC = extractedAC
plt.imshow(abs(extractedAC), cmap="gray")  # show 32x32 kernel
plt.show()

print("Generate gaussian noise 128x128 with SD of 90")
gaussNoise = np.random.normal(0, 90,(128,128)) # genereate Gaussian noise 128x128
plt.imshow(gaussNoise, cmap="gray")  # set the color map to bone
plt.show()

print("Convolve the Gaussian noise with the 32x32 autocorrelation kernel to obtain noise pattern spatially correlated with the noise in the phantom scan")
# convolve the Gaussian noise with the 32x32 autocorrelation kernel
spatialNoise = signal.convolve2d(gaussNoise, abs(extractedAC))

plt.imshow(spatialNoise, cmap="gray")  # set the color map to bone
plt.show()

问题原因与修复方案 #

目前生成的空间相关统计噪声过于模糊，与体模扫描的频谱特性不匹配，由以下几个代码错误导致：

1. 傅里叶变换逻辑错误 #

计算噪声频谱时误用了逆傅里叶变换ifft2，需替换为正傅里叶变换fft2：

# 原错误代码
dataNoiseFT = np.fft.ifft2(dataNoise)
# 修改为
dataNoiseFT = np.fft.fft2(dataNoise)

2. 自相关核未做归一化 #

提取的32x32自相关核没有归一化，卷积时相当于对高斯噪声做了过度平滑，直接导致噪声模糊。需要对核做归一化处理，保证核的总权重为1，避免噪声能量异常缩放：

# 提取核后新增归一化代码
extractedAC = abs(shiftedAC)[int(n*(3/8)):int(n*(3/8))+int(n/4), int(n*(3/8)):int(n*(3/8))+int(n/4)]
extractedAC = extractedAC / np.sum(extractedAC)

3. 卷积后噪声未匹配目标标准差 #

卷积操作会改变噪声的标准差，未做后续缩放调整无法匹配体模噪声的实际强度。卷积完成后需将噪声缩放到目标标准差：

# 可根据需求修改目标标准差，示例为匹配体模噪声的标准差
target_std = np.std(dataNoise)
# 建议使用same模式保持输出尺寸与输入噪声一致
spatialNoise = signal.convolve2d(gaussNoise, extractedAC, mode='same')
spatialNoise = spatialNoise / np.std(spatialNoise) * target_std

4. 自相关核尺寸调整 #

如果修正以上三点后频谱仍不匹配，可尝试保留更大尺寸的自相关核，32x32可能截断了有效相关范围，导致频谱特性偏差。

内容的提问来源于stack exchange，提问作者user3220072