U

#### user3220072

##### Guest

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}")
```

<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>