Calculating Power of Specific Spatial Frequencies of pictures

50 views (last 30 days)
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)]);
  3 Comments
Olivia
Olivia on 31 Jul 2024
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
Umar on 31 Jul 2024
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.

Sign in to comment.

Answers (2)

Raghava S N
Raghava S N on 31 Jul 2024
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-
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

Paul
Paul on 2 Aug 2024
Edited: Paul on 2 Aug 2024
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.

Categories

Find more on Fourier Analysis and Filtering in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!