OiO.lk Community platform!

Oio.lk is an excellent forum for developers, providing a wide range of resources, discussions, and support for those in the developer community. Join oio.lk today to connect with like-minded professionals, share insights, and stay updated on the latest trends and technologies in the development field.
  You need to log in or register to access the solved answers to this problem.
  • You have reached the maximum number of guest views allowed
  • Please register below to remove this limitation

Estimating noise (std) of pixel values from NPS?

  • Thread starter Thread starter user3220072
  • Start date Start date
U

user3220072

Guest
I have written python code to open a DICOM file of a CT phantom, extract the center (homogenous) 128x128 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?

Code:
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}")
<p>I have written python code to open a DICOM file of a CT phantom, extract the center (homogenous) 128x128 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?</p>
<pre><code>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}")
</code></pre>
 
Top