File Exchange

## Elliptic fourier for shape analysis

version 1.3.0.0 (6.65 KB) by Auralius Manurung

### Auralius Manurung (view profile)

Implementation of elliptic fourier for shape analysis.

Updated 19 May 2016

1) plot_chain_code(ai, color, line_width)
This function will plot the given chain code. The chain code (ai) should be in
column vector.
Example:
>> ai = [5 4 1 2 3 4 3 0 0 1 0 1 0 0 0 7 7 1 1 0 7 5 4 5 4 5 0 6 5 4 1 3 4 4 4 4 6];
>> plot_chain_code(ai)

2) plot_fourier_approx(ai, n, m, normalized, color, line_width)

This function will plot the Fourier approximation, given a chain code (ai),
number of harmonic elements (n), and number of points for reconstruction (m).
Normalization can be applied by setting "normalized = 1".

3) output = calc_traversal_dist(ai, n, m, normalized)

This function will generate position coordinates of chain code (ai). Number of
harmonic elements (n), and number of points for reconstruction (m) must be
specified.
The output is a matrix of [x1, y1; x2, y2; ...; xm, ym].

3) output = fourier_approx(ai, n, m, normalized)

This function will generate position coordinates of Fourier approximation of
chain code (ai). Number of harmonic elements (n), and number of points for
reconstruction (m) must be specified.
The output is a matrix of [x1, y1; x2, y2; ...; xm, ym].

4) output = calc_harmonic_coefficients(ai, n)

This function will calculate the n-th set of four harmonic coefficients.
The output is [an bn cn dn]

5) [A0, C0] = calc_dc_components(ai)

This function will calculate the bias coefficients A0 and C0.

6) output = calc_traversal_dist(ai)

Traversal distance is defined as accumulated distance travelled by every
component of the chain code assuming [0 0] is the starting position.
Example:
>> x = calc_traversal_dist([1 2 3])
x =
1 1
1 2
0 3

7) output = calc_traversal_time(ai)

Traversal time is defined as accumulated time consumed by every
component of the chain code.
Example:
>> x = calc_traversal_time([1 2 3])
x =

1.4142
2.4142
3.8284

### Cite As

Auralius Manurung (2020). Elliptic fourier for shape analysis (https://www.mathworks.com/matlabcentral/fileexchange/32800-elliptic-fourier-for-shape-analysis), MATLAB Central File Exchange. Retrieved .

bentozer anass

Braiki Marwa

### Braiki Marwa (view profile)

Hello Auralius, thanks for sharing this code. I have a problem when I use the code for the images that have two or more objects. How can change this code to resolve this problem. I tried with this modification, but it still the same problem!!

function output = fourier_approx(ai, n, m, normalized,xi,yi)

% This function will generate position coordinates of fourier approximation of
% chain code (ai).Number of harmonic elements (n), and number of points for
% reconstruction (m) must be specified.

for i = 1 : n
harmonic_coeff = calc_harmonic_coefficients(ai, i);
a(i) = harmonic_coeff(1, 1);
b(i) = harmonic_coeff(1, 2);
c(i) = harmonic_coeff(1, 3);
d(i) = harmonic_coeff(1, 4);
end

[A0, C0] = calc_dc_components(ai);

% Normalization procedure
if (normalized == 1)
% Remove DC components
A0 = 0;
C0 = 0;

% Compute theta1
theta1 = 0.5 * atan(2 * (a(1) * b(1) + c(1) * d(1)) / ...
(a(1)^2 + c(1)^2 - b(1)^2 - d(1)^2));

costh1 = cos(theta1);
sinth1 = sin(theta1);

a_star_1 = costh1 * a(1) + sinth1 * b(1);
b_star_1 = -sinth1 * a(1) + costh1 * b(1);
c_star_1 = costh1 * c(1) + sinth1 * d(1);
d_star_1 = -sinth1 * c(1) + costh1 * d(1);

% Compute psi1
psi1 = atan(c_star_1 / a_star_1) ;

% Compute E
E = sqrt(a_star_1^2 + c_star_1^2);

cospsi1 = cos(psi1);
sinpsi1 = sin(psi1);

for (i = 1 : n)
normalized = [cospsi1 sinpsi1; -sinpsi1 cospsi1] * [a(i) b(i); c(i) d(i)] * ...
[cos(theta1 * i) -sin(theta1 * i); sin(theta1 * i) cos(theta1 * i)];

a(i) = normalized(1,1) / E;
b(i) = normalized(1,2) / E;
c(i) = normalized(2,1) / E;
d(i) = normalized(2,2) / E;
end

end % end if normalized

for t = 1 : m
x_ = 0.0;
y_ = 0.0;

for i = 1 : n
x_ = x_ + (a(i) * cos(2 * i * pi * t / m) + b(i) * sin(2 * i * pi * t / m));
y_ = y_ + (c(i) * cos(2 * i * pi * t / m) + d(i) * sin(2 * i * pi * t / m));
end

output(t,1) = xi+A0 + x_;
output(t,2) = yi+C0 + y_;
end

end

Auralius Manurung

### Auralius Manurung (view profile)

Hi Richard,
Yes, though there is a limit to it.
Yo can play with the 'n' and the 'm' parameters.
It should becoming closer and closer to the original image.

Richard

### Richard (view profile)

Hi Auralius, Thanks for making this. I see a problem when I run example 1. If I change n increasing it from 10 the resulting approximation plot does not continue to get closer to the original image. For example n of 100 is pretty much the same as n = 1000. Worse yet there is quite a bit of visible error around the figure.

Shouldn't this approach the original image as n is greatly increased?

thanks

Auralius Manurung

### Auralius Manurung (view profile)

Auralius Manurung

### Auralius Manurung (view profile)

Hi Chrstopher Cramer,

Yes! You are absolutely right!

Chrstopher Cramer

### Chrstopher Cramer (view profile)

Hi Auralius, nice work. I think I've found a couple of additional errors. In both calc_harmonic_coefficients.m and calc_dc_components.m, you have a conditional:

if (p > 2)

in all of these cases, this should either be:

if (p > 1)

or maybe

if (p >= 2)

Auralius Manurung

### Auralius Manurung (view profile)

Hi João Neves,

Sorry for the late reply and thank you for your valuable inputs.

You are definitely right. I will fix this in the next update.

Cheers,
Auralius

João Neves

### João Neves (view profile)

Fantastic work, however I would like to know if there is an error in the calc_dc_components function.
Is this line correct?

sum_c0 = sum_c0 + delta_x / (2 * delta_t) * (t(p))^2 + delta * t(p);

or it should be

sum_c0 = sum_c0 + delta_y / (2 * delta_t) * (t(p))^2 + delta * t(p);

Thanks again for this great work

 19 May 2016 1.3.0.0 Updated as per feedback from Chrstopher Cramer. 24 Nov 2012 1.3.0.0 Corrected mistake on the equation based on João Neves' suggestion. 8 Sep 2011 1.1.0.0 Added function descriptions for better understanding
##### MATLAB Release Compatibility
Created with R2009a
Compatible with any release
##### Platform Compatibility
Windows macOS Linux