Using FFT Coefficients/Descriptors to reconstruct particle shape

Hi, I am trying to use Fourier descriptors to reconstruct the shape of a particle. I have written most of the code (with a lot of help) and only the last step which involves using Fourier descriptors to reconstruct the shape of the original image remains.

I need help to know which functions to use to achieve this. I would really appreciate if anyone could help.

Here is the code:

I have also attached the original image and another file which demonstrates what I need to do.

scale_factor = 100/90;           % (units: microns/pixel)
                                 % to convert pixels to microns,
                                 % multiply # of pixels by scale_factor
%1. Open and show image of particle
img = imread('1_50.JPG');        %Read image
subplot (3,3,1)
imshow(img)                      %show image

%2. Binarize and show image
BW = im2bw(img,0.45);            %binarize image with a threshold value of 0.45
img1= bwareaopen (BW, 1000);     %Remove small objects
img2= imfill(img1, 'holes');     %fill holes
subplot(3,3,2)                   
imshow(img2)                     %show binarized image

%3. Setting up the axes 
% for plotting into multiple axes in one loop, it's convenient to set up
% all the axes first, before anything is plotted into them

% 3.1. set up the 2nd axes:
ax3 = subplot(3,3,3);
hold on
box on                        %draws a box that encompasses the entire axes area
axis image                    %sets the aspect ratio of the axes (ax2) to be equal and 
                              %adjusts the limits of the axes so that the data units are 
                              %of equal size along both the x-axis and y-axis.
set(ax3,'YDir','reverse')     %reverses the y-axis to begin from the top
                              %since the origin of the image is at the top left-corner
% 3.2. set up the 3rd axes:
ax4 = subplot(3,3,4);
hold on

%4. Compute centroid of binarized image 
centriod_value= regionprops(img2, 'Centroid');  %Finds centroid of image 
centroid = cat(1,centriod_value.Centroid);      %Stores centroid coordinates in a two-column matrix

%5. Tracing the boundary of image
p_boundary= bwboundaries(img2);                % trace the exterior boundaries of particle
                                               %(in row, column order, not x,y)
number_of_boundaries = size(p_boundary,1);    %specify only one boundary
centroid = centroid*scale_factor;

% 6. Setting up equally spaced angles before initing the loop
%use 129 points to have an angle spacing of 360/128
N_angles = 129;
interp_angles = linspace(0,360,N_angles).';

%7. Initiate for loop to plot boundary, centroid, & a graph of radius vs angle

for k = 1 : number_of_boundaries              % initiate loop to go round the selected boundary
    thisBoundary = p_boundary{k}*scale_factor;
    y = thisBoundary(:,1);                    % rows
    x = thisBoundary(:,2);                    % columns
    
    plot(ax3, x, y, 'g', 'LineWidth', 2);
    
    plot(ax3, centroid(:,1),centroid(:,2),'b.')
    
    % 7.1. Calculate the angles in degrees
    deltaY = thisBoundary(:,1) - centroid(k,2); % boundary(:,1) is y, but centroid(k,2) is y
    deltaX = thisBoundary(:,2) - centroid(k,1); % boundary(:,2) is x, but centroid(k,1) is x
    % angles = atand(deltaY ./ deltaX); % atand gives angles in range [-90,90]
    angles = atan2d(deltaY,deltaX); % use atan2d to get angles in range [-180,180]
    
    % 7.2. Calculate the radii
    radius = sqrt(deltaY.^2 + deltaX.^2); % use deltaX and deltaY, which were just calculated
    
    % 7.3. Add 360 degrees to negative angles so that range is [0,360] instead of [-180,180]
    idx = angles < 0;                  %idx reps angles less than 0 degreees
    angles(idx) = angles(idx)+360;     %360 degrees is added to idx angles
    
    %7.4. Sort angles in ascending order
    % sort the angles, so the first one is smallest (closest to zero)
    % and the last one is largest (closest to 360)
    % this is useful for getting a smooth line that goes from 0 to 360 in the plot
    % I returns the corresponding indices of the angles
    [angles,I] = sort(angles);
    
    %7.5. Reorder distances the same way angles was reordered using I
 
    radius = radius(I);
    
    %7.6. Use only unique angles in case an angle is repeated
    % in case any angles are repeated, use only the unique set and
    % corresponding distances in interp1 (in this data, one angle is
    % repeated twice, with the same distance each time)
    %The ~ (tilde) here is a placeholder indicating that we're not interested 
    %in capturing the actual unique elements, only their indices
    [~,I] = unique(angles);
    
    % 7.7. perform linear interpolation/extrapolation:
    %angles(I)= sample points, radius(I)= correspoinding points
    %and interp_angles = coordinates of the querry points (QP) 
    %In general, the functions returns interpolated values at specific QP
    interp_distances = interp1(angles(I),radius(I),interp_angles,'linear','extrap');
    
    % 7.8 Plot distance vs. angle.
    plot(ax4,interp_angles,interp_distances)
    
end

%8. Apply fft algorithm to values of radius 
L = length(interp_distances);    %Define the number of the radiuses (points)
f = 0:L-1;             %Define the sampling frequencies
A = fft(interp_distances);       %Transforming the radius function into frequency domain using FFT
B = (abs(A)/2535.2);      % Normalize the amplitude of fourier descriptors using the first descriptors

%9. Plot fft descriptor vs. frequency
subplot(3,3,5)
bar(f,B,3)       %plot the frequencies vs. the nornalised descriptors
xlim([0 128])

%10. Plot phase angle vs. frequency
subplot(3,3,6)
bar(f, angle(A),3)   % plot the frequencies vs phase angles 
xlim([0 128])
ylim([-4 4])

%11. Create a table of values 
table1 = table(A, f', abs(A), B, angle(A));
table1.Properties.VariableNames = {'FFT_coeffs', 'Frequency', 'Amplitude', 'Norm_Amplitude' 'Phase'};
disp(table1);

%12. Reconstruct particle boundary using FFT descriptors ?????

Answer 1

[Suraj Kumar]

Hi Johnson,

As I interpret it, you aim to use Fourier descriptors for reconstructing a particle's boundary. To accomplish this, you can consider following these steps and refer to the attached code snippets:

1. Use the “ifft” function to perform the inverse Fast Fourier Transform on the coefficients A. This step converts the frequency domain data back into the spatial domain, giving reconstructed radii of the particle boundary.

reconstructed_radius = ifft(A);

2. Convert these radii into x and y coordinates using the centroid of the particle and the interpolated angles. The "cosd" and "sind" functions compute the cosine and sine of the angles respectively and are multiplied by the reconstructed radii to get the x and y coordinates.

reconstructed_x = centroid(1,1) + reconstructed_radius .*    cosd(interp_angles);
reconstructed_y = centroid(1,2) + reconstructed_radius .* sind(interp_angles);

3. Plot the reconstructed particle boundary in a 3x3 grid of subplots. The boundary is visualized in red, with the centroid marked in blue.

subplot(3,3,7)
plot(reconstructed_x, reconstructed_y, 'r', 'LineWidth', 2);
hold on
plot(centroid(:,1), centroid(:,2), 'b.')
axis image
set(gca, 'YDir', 'reverse')
title('Reconstructed Boundary')
xlabel('X (microns)')
ylabel('Y (microns)')

Refer to the output for better understanding:

Please refer to the documentation links below:

• ifft- https://www.mathworks.com/help/matlab/ref/ifft.html
• sind- https://www.mathworks.com/help/matlab/ref/sind.html
• cosd- https://www.mathworks.com/help/matlab/ref/cosd.html

Happy Coding!

[Johnson]

Hi Suraj. This is exactly what I need. Thank you so much!

[Suraj Kumar]

You are welcome.