I’m trying to estimate the SD of pixel values in the spatial domain from the 2D NPS in Python.

I would expect to need the sum of the NPS values divided by the total number of pixels. However I only reach to correct estimate by taking the mean divided by the total number of pixels.

Can anyone point me in the direction as to why that is?

Please see code examples below.

“SUM:”

import numpy as np
from scipy.fftpack import fft2, fftshift

# Generate a random 128x128 image
image = np.random.rand(128, 128)*100

# Subtract the mean of the image
image = image - np.mean(image)

# Compute the 2D Fourier transform of the image
F = fft2(image)

# Shift the zero frequency component to the center
F_shifted = fftshift(F)

# Calculate the Noise Power Spectrum (NPS)
NPS = np.abs(F_shifted)**2

# Calculate the total power in the NPS
total_power = np.sum(NPS)

# The total power in the NPS corresponds to the sum of the squared deviations
# from the mean in the spatial domain. For an image of size NxN:
N = image.shape[0]
variance_from_nps = total_power / (N * N)

# Compute the standard deviation from the variance
sigma_from_nps = np.sqrt(variance_from_nps)

# Compute the standard deviation directly from the image
sd_direct = np.std(image)

print(f"Standard Deviation from NPS: {sigma_from_nps}")
print(f"Standard Deviation directly from the image: {sd_direct}")

Standard Deviation from NPS: 3692.2142651255785
Standard Deviation directly from the image: 28.84542394629358

“Mean:”

import numpy as np
from scipy.fftpack import fft2, fftshift

# Generate a random 128x128 image
image = np.random.rand(128, 128)*100

# Subtract the mean of the image
image = image - np.mean(image)

# Compute the 2D Fourier transform of the image
F = fft2(image)

# Shift the zero frequency component to the center
F_shifted = fftshift(F)

# Calculate the Noise Power Spectrum (NPS)
NPS = np.abs(F_shifted)**2

# Calculate the total power in the NPS
total_power = np.mean(NPS)

# The total power in the NPS corresponds to the sum of the squared deviations
# from the mean in the spatial domain. For an image of size NxN:
N = image.shape[0]
variance_from_nps = total_power / (N * N)

# Compute the standard deviation from the variance
sigma_from_nps = np.sqrt(variance_from_nps)

# Compute the standard deviation directly from the image
sd_direct = np.std(image)

print(f"Standard Deviation from NPS: {sigma_from_nps}")
print(f"Standard Deviation directly from the image: {sd_direct}")

Standard Deviation from NPS: 28.894608971414826
Standard Deviation directly from the image: 28.894608971414826