Calculating Power of Specific Spatial Frequencies of pictures

Hi experts,

I want to calculate the power of the spatial frequency range 20-30 cycles per millimeter in an picture. I am not sure how I can accomplish this using MATLAB code.

I downloaded following code online and made some modifications. Could you check the code for me and tell me whether I am correct?

FYI. "Cycles per millimeter" (c/mm) is a measure of spatial frequency, indicating how many repeating structures (cycles) there are per millimeter in an image. This metric is commonly used in image processing to describe the resolution or detail level of an image.

% Load the image
image = imread('image.png');
image = rgb2gray(image);

% Get the size of the image
[rows, cols] = size(image);

% Compute the Fourier transform of the image
fft_image = fft2(double(image));

% Shift zero frequency component to the center of the spectrum
fft_image_shifted = fftshift(fft_image);

% Compute the magnitude spectrum
magnitude_spectrum = abs(fft_image_shifted);

% Create frequency axis
u = (-cols/2:cols/2-1) * (1 / (cols * 1e-3)); % Frequency in cycles per millimeter
v = (-rows/2:rows/2-1) * (1 / (rows * 1e-3)); % Frequency in cycles per millimeter

% Create meshgrid for frequencies
[U, V] = meshgrid(u, v);

% Calculate the radial frequency
radial_freq = sqrt(U.^2 + V.^2);

% Define the frequency range
freq_min = 20; % 20 cycles per millimeter
freq_max = 30; % 30 cycles per millimeter

% Create a mask for the frequency range
frequency_mask = (radial_freq >= freq_min) & (radial_freq <= freq_max);

% Compute the power in the specified frequency range
energy = sum(magnitude_spectrum(frequency_mask).^2, 'all');

% Display the result
disp(['Power in the frequency range 20-30 cycles per millimeter: ', num2str(energy)]);

[Umar]

Hi @Olivia,

To address your query, “I want to calculate the power of the spatial frequency range 20-30 cycles per millimeter in a picture. I am not sure how I can accomplish this using MATLAB code”, please see updated code snippet below. So, basically this step is required “Convert to Grayscale, if it's a color image for frequency analysis” in which your case it is. Then, apply Fourier Transform by computing the 2D Fourier Transform of the grayscale image. Afterwards, use Shift Zero Frequency to the center for better visualization then Calculate Frequency Range by determining the frequency range in cycles per millimeter based on the image resolution.Finally, Calculate Power by counting the number of pixels in the specified frequency range to calculate the power.

% Read the image

image = imread('/MATLAB Drive/IMG_7236.jpg');

% Convert the image to grayscale

gray_image = rgb2gray(image);

% Perform Fourier Transform on the image

ft_image = fft2(double(gray_image));

% Shift the zero frequency component to the center

ft_image_shifted = fftshift(ft_image);

% Calculate the frequency range in cycles per millimeter

% Assuming the image resolution in pixels per millimeter

resolution = 10; % For example, 10 pixels per millimeter

freq_range = [20 30] / resolution;

% Calculate the power in the specified frequency range

power = sum(sum(abs(ft_image_shifted) >= freq_range(1) & abs(ft_image_shifted) <= freq_range(2)));

disp(['Power of spatial frequency range 20-30 cycles per millimeter: ', num2str(power)]);

This is basically how I will approach to solve the problem. Please see attached results and feel free to customize this code based on your preference.

Hope, this helps resolve your problem. Please let me know if you have any further questions.

[Olivia]

Hi Umar,

Thanks a lot for your answer.

But the line confuses me a lot:

power = sum(sum(abs(ft_image_shifted) >= freq_range(1) & abs(ft_image_shifted) <= freq_range(2)));

Because this code:

abs(ft_image_shifted) 

calculates the power.

Does it mean anything to compare the calculating results with our frequence cutoff?

[Umar]

Hi @Olivia,

To address your query, The line of code provided calculates the power within the spatial frequency range of 20-30 cycles per millimeter in an image. The expression abs(ft_image_shifted) computes the magnitude of the shifted Fourier transform of the image. The code segment abs(ft_image_shifted) >= freq_range(1) & abs(ft_image_shifted) <= freq_range(2) filters out the frequencies falling within the specified range. By summing the logical array resulting from this comparison, the total power within the desired frequency range is obtained. Comparing this calculated power with the frequency cutoff allows for assessing the significance of the spatial frequencies in the image within the specified range. Please let me know if you have any further questions.

Answer 1

[Raghava S N]

Hi Olivia,

Please refer to this answer provided by a MathWorks staff member that describes in detail the workflow that uses Discrete Fourier Transforms(DFTs) to determine the spatial frequency of patterns in images - https://www.mathworks.com/matlabcentral/answers/13896-fftshift-of-an-image#:~:text=In%20order%20to,of%20physical%20distance.

Concretely, the answer describes the following workflow:

The frequency domain for the cut-off spatial frequencies must be created first. This domain should have the units “cycles per millimeter” in this case.

Following which, the functions “fft2” and “fftshift” are used to compute the 2D DFT of the image. Please refer to the documentation of these functions for more information-

“fft2” - https://www.mathworks.com/help/matlab/ref/fft2.html

“fftshift” - https://www.mathworks.com/help/matlab/ref/fftshift.html

After this, displaying the resulting spectral image for the cut-off spatial frequencies will show the peak values in the spectrum.

Finally, the energy can be calculated from the image as the code does in your MATLAB Answers post.

Please refer to another MATLAB Answers post that discusses the same workflow for more details - https://www.mathworks.com/matlabcentral/answers/14587-how-to-find-wavelength-of-a-wave-in-an-image.

Hope this helps!

Raghava

Answer 2

[Paul]

Hi Olivia,

The equations for u and v are correct if rows and cols are both even. Otherwise, they'd have to be adjusted.

Regarding this line

% Compute the power in the specified frequency range
energy = sum(magnitude_spectrum(frequency_mask).^2, 'all');

The comment should change "power" to "energy", and the computation needs to be scaled

% Compute the energy in the specified frequency range
energy = sum(magnitude_spectrum(frequency_mask).^2, 'all')/(rows*cols);

Keep in mind that's only an approximation to the energy constained within the specified frequency region.