通过NPS估算CT体模像素值标准差结果偏差问题咨询

你的结果偏差100倍，核心是NPS计算的两个关键步骤缺失：

1. 未去除图像均值（直流分量） #

NPS是噪声的功率谱，需要先从图像中减去均值，只保留噪声部分。如果直接对原始图像做FFT，图像的直流分量（均值对应的频率分量）会被计算进去，这部分功率远大于噪声，直接拉高了NPS的平均值。

2. FFT结果未做归一化 #

np.fft.fft2的输出幅度和图像的像素总数（128×128）成正比，计算功率谱时需要除以像素总数的平方（或在FFT后直接除以像素总数），否则功率谱数值会被放大128²倍，最终开平方后结果会放大128倍（和你看到的100倍偏差吻合）。

import pydicom
import numpy as np
import matplotlib.pyplot as plt

def load_dicom(file_path):
    """ Load DICOM image from the file path. """
    dicom_file = pydicom.dcmread(file_path)
    return dicom_file.pixel_array, dicom_file.RescaleIntercept

def extract_center_pixels(image, size=128):
    """ Extract the center size x size pixels from the image. """
    center_x, center_y = image.shape[1] // 2, image.shape[0] // 2
    return image[center_y - size // 2 : center_y + size // 2, center_x - size // 2 : center_x + size // 2]

def calculate_nps(image_section):
    """ Calculate the 2D Noise Power Spectrum (NPS). """
    # 去除均值，只保留噪声分量
    noise = image_section - np.mean(image_section)
    # 做FFT并归一化，消除像素数量对幅度的影响
    f_transform = np.fft.fft2(noise) / noise.size
    f_shifted = np.fft.fftshift(f_transform)
    # 计算功率谱（幅度平方）
    magnitude_squared = np.abs(f_shifted) ** 2
    return magnitude_squared

def plot_nps(nps):
    """ Plot the 2D NPS. """
    plt.figure(figsize=(10, 8))
    plt.imshow(np.log(nps + 1), cmap='gray')
    plt.colorbar()
    plt.title('2D Noise Power Spectrum')
    plt.show()

# 主流程
image, intercept = load_dicom('/content/drive/MyDrive/Colab Notebooks/PhantomDII/B4017.dcm')
image = image.astype(np.float64) + intercept
center_image = extract_center_pixels(image)

# 直接计算标准差
print(f"直接计算的像素值标准差: {np.std(center_image)}")

# 计算NPS并估算标准差
nps = calculate_nps(center_image)
std_from_nps = np.sqrt(np.mean(nps))
print(f"通过NPS估算的像素值标准差: {std_from_nps}")

# 绘图部分
plt.figure(figsize=(10, 8))
plt.imshow(center_image, cmap='gray')
plt.colorbar()
plt.title('提取的中心区域')
plt.show()

plt.figure(figsize=(10, 8))
plt.imshow(image, cmap='gray')
plt.colorbar()
plt.title('原始DICOM图像')
plt.show()

plot_nps(nps)

• 去均值：noise = image_section - np.mean(image_section)，确保仅对噪声做FFT，排除直流分量的干扰。
• FFT归一化：np.fft.fft2(noise) / noise.size，消除像素数量对FFT幅度的放大作用，得到正确的功率谱数值。
• 修正了原代码中错误的图像标题（将原始图像标题误写为2D NPS）。

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