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 ?????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 ????? 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 ????? fft, shape reconstruction, fourier descriptor MATLAB Answers — New Questions