I have written python code to open a DICOM file of a CT phantom, extract the center (homogenous) 128×128 pixels. Subsequently calculate the std of these pixels. Finally I’m calculating 2D NPS and using the squared 2D NPS mean to estimate the std of the original pixel values. However, the estimate I get with the code below is x 100 larger than the direct std computation from the pixel values. What am I doing wrong?

import pydicom
import numpy as np
from scipy.fftpack import fftshift, fft2
import matplotlib.pyplot as plt
t_transform = 0
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). """
    f_transform = np.fft.fft2(image_section)

    
    f_shifted = np.fft.fftshift(f_transform)  # Shift zero frequency component to the center
    
    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"Standard Deviation of the Pixel Values by direct calculation:")

print(np.mean(center_image))

print(np.std(center_image))

nps = calculate_nps(center_image)

plt.figure(figsize=(10, 8))
plt.imshow(center_image, cmap='gray')
plt.colorbar()
plt.title('Extracted')
plt.show()

plt.figure(figsize=(10, 8))
plt.imshow(image, cmap='gray')
plt.colorbar()
plt.title('2D Noise Power Spectrum')
plt.show()

plot_nps(nps)
# Calculate the standard deviation from NPS
std_deviation = np.sqrt(np.mean(nps))
print(f"Standard Deviation of the Pixel Values calculated from the NPS: {std_deviation}")