使用自相关和 python 生成空间相关的统计噪声

Question

我正在尝试模仿本文中描述的方法： http://www.ajnr.org/content/34/8/1506 在“模拟剂量减少的噪声添加”标题下 据我了解，作者基本上完成了以下步骤。

1. 首先，从体模数据中测量 CT 噪声的光谱特性
2. 计算噪声自相关函数是根据这些数据计算得出的
3. 围绕自相关峰值生成窗口并保存为卷积滤波器
4. 要生成具有适当功率谱的噪声，请将此滤波器应用于高斯白噪声并缩放至所需的标准偏差。

随后可以将生成的空间相关噪声添加到患者的真实 CT 图像中，以生成具有与体模扫描相同光谱特性的噪声的 CT 图像。

幻像扫描的dicom文件链接：https://drive.google.com/file/d/1ouFfqSxWo7PFV4CXhYI1TEbNZp-vM8Bq/view?usp=sharing

带输出的 Goole colabhttps://colab.research.google.com/drive/1utFiDNbElyeKGuyeHg3_rVG4ZDc5xoHQ?usp=sharing

下面还有python代码。

我能够生成的空间相关统计噪声似乎过于模糊，与 Phantom 扫描中发现的光谱特性不相符。

我想知道是否有人能看出我哪里出错了？

亲切的问候 /////////

#! 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()

Answer 1

我的项目遇到了同样的问题！对于最后一个卷积，请只对测量的 ACF 进行卷积，而不是对 ACF 的绝对值进行卷积。 ACF的标志也有信息。