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

Posted 2023-02-16

【中文标题】使用自相关和 python 生成空间相关的统计噪声

【英文标题】：Generating spatially correlated statistical noise using autocorrelation and python

【发布时间】：2021-11-01 18:31:30

【问题描述】：

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

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

随后可以将生成的空间相关噪声添加到患者的真实 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()

【问题讨论】：

会不会是像素比例的原因？即您使用的是像素，但纸张中的比例尺是毫米，反之亦然？

嗨，理查德，感谢您回复我的信息。包含噪声数据的原始 DICOM 文件为 128x128，与生成的自相关函数卷积的高斯噪声也是 128x128，因此应该预期图像纹理相似。因此，我看不出不同的像素比例应该如何解释原始噪声数据和生成的噪声数据之间的纹理差异。你同意吗？

我不确定。我快速浏览了这篇论文，但乍一看并不太明显他们在做什么/如何做事。我知道大多数 CT 是 256x256，而不是 128。我知道 CT 中的噪声不是高斯的，需要通过考虑比例的幂函数（或卷积）类型模型来建模。因此，如果您只是在像素单位中工作，并且没有正确地将所有内容缩放到物理单位，那么我认为您会遇到问题。

嗨，理查德，再次感谢您提供建议。我应该纠正自己并说我从 512x512 像素的幻像 DICOM 文件的统一部分中提取了 128x128 像素的 ROI。 ACF 是由上述噪声产生的。随后，ACF 与具有高斯噪声的 128x128 图像进行卷积，我预计输出图像的 (1D) 功率谱类似于上述提取的 128x128 噪声中的功率谱。

【参考方案1】：

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

【讨论】：

您好，Taka，非常感谢您与我们联系，对于延迟回复您的回答，我们深表歉意。我不确定你所说的测量 ACF 到底是什么意思？ abs 命令从 extractAC 2D 数组的复数中获取幅​​度。复数的大小有符号吗？亲切的问候