image thumbnail

Discrete Curvature, Normals and Evolute

version 6.0.2 (3.14 KB) by Moreno, M.
Fast calculation of the discrete 1D curvature, evolute and normals of a 1D, 2D or 3D curve.

55 Downloads

Updated 4 Apr 2022

View License

Rapid calculation of discrete 1D curvature 'k', normal 'n' and evolute 'e' for a curve 'X' in the two or three-dimensional Euclidean space. Optional parameter 'f' for rotating frame corrections for preventing normal inversion at the inflexion points of X. Calling [k, n, e] = kappa(X, f) returns:
  • Discrete 1-D curvature vector 'k' calculated as the inverse of the radius of the circumscribing circle for every triplet of points in X. The end-values of the curvature are corrected with linear mid-point extrapolation.
  • Normals 'n' of the curve X calculated as the normalised difference between X and its evolute.
  • Evolute 'e' of the curve X calculated as the locus of the centres of the previously-defined circumscribing circles.
Use of matrix form for maximising the speed of calculations and avoiding unnecessary pre-allocations. Some examples of the use of this function are:
% 2D Cardioid
t = 2 * pi * (0 : 1 / (1000 - 1) : 1)';
X = [(0.5 + cos(t)) .* sin(t), (0.5 + cos(t)) .* cos(t)];
[k, n, e] = kappa(X);
K = X + k .* n / max(k);
figure
hold on
plot(X(:, 1), X(:, 2))
plot(K(:, 1), K(:, 2))
plot(e(:, 1), e(:, 2))
axis equal
xlim([-1, 1])
ylim([-0.75, 1.75])
figure
hold on
h = plot(X(:, 1), X(:, 2), 'LineWidth', 1.2, 'Color', '#0072BD');
h.Color(4) = 0.8;
h = plot(K(:, 1), K(:, 2), 'LineWidth', 1.0, 'Color', '#0072BD');
h.Color(4) = 0.6;
for i = 1 : 1000 - 1
h = plot([X(i, 1), K(i, 1)], [X(i, 2), K(i, 2)], 'Color', '#0072BD');
h.Color(4) = 0.1;
end
axis equal
xlim([-1, 1])
ylim([-0.75, 1.75])
figure
hold on
h = plot(X(:, 1), X(:, 2), 'LineWidth', 1.2, 'Color', '#D95319');
h.Color(4) = 0.8;
h = plot(e(:, 1), e(:, 2), 'LineWidth', 1.0, 'Color', '#D95319');
h.Color(4) = 0.6;
for i = 1 : 1000 - 1
h = plot([X(i, 1), e(i, 1)], [X(i, 2), e(i, 2)], 'Color', '#D95319');
h.Color(4) = 0.1;
end
axis equal
xlim([-1, 1])
ylim([-0.75, 1.75])
% 3D Derived Lissajous
t = 2 * pi * (0 : 1 / (1000 - 1) : 1)';
X = [cos(t), sin(t), 1.2 * cos(2 * t)];
[k, n, ~] = kappa(X);
K = X + k .* n / 8;
E = 1 ./ k;
E = X + E .* n / 16;
figure
hold on
plot3(X(:, 1), X(:, 2), X(:, 3))
plot3(K(:, 1), K(:, 2), K(:, 3))
plot3(E(:, 1), E(:, 2), E(:, 3))
view (3)
axis equal
xlim([-1, 1])
ylim([-1, 1])
figure
hold on
for i = 1 : 1000 - 1
plot3([X(i, 1), K(i, 1)], [X(i, 2), K(i, 2)], [X(i, 3), K(i, 3)], ...
'Color', [0.9 0.9447 0.9741]);
end
plot3(X(:, 1), X(:, 2), X(:, 3), 'LineWidth', 1.2, 'Color', ...
[0.2 0.5576 0.7928]);
plot3(K(:, 1), K(:, 2), K(:, 3), 'LineWidth', 1.0, 'Color', ...
[0.4 0.6682 0.8446]);
view(3)
axis equal
xlim([-1, 1])
ylim([-1, 1])
figure
hold on
for i = 1 : 1000 - 1
plot3([X(i, 1), E(i, 1)], [X(i, 2), E(i, 2)], [X(i, 3), E(i, 3)], ...
'Color', [0.985 0.9325 0.9098]);
end
plot3(X(:, 1), X(:, 2), X(:, 3), 'LineWidth', 1.2, 'Color', ...
[0.88 0.46 0.2784]);
plot3(E(:, 1), E(:, 2), E(:, 3), 'LineWidth', 1.0, 'Color', ...
[.91 0.595 0.4588]);
view(3)
axis equal
xlim([-1, 1])
ylim([-1, 1])
% Correction of Frenet's frame of reference with rotating frame motion
b = (-1 : 2 / (1000 - 1) : 1)';
b = [b, b .^ 2 .* (b - 1)];
[k, n] = kappa(b, 0);
[K, N] = kappa(b, 1);
k = b + k .* n / max(k) / 10;
K = b + K .* N / max(K) / 10;
figure
hold on
plot(b(:, 1), b(:, 2), 'Color', 'k')
plot(k(:, 1), k(:, 2), 'Color', '#0072BD')
plot(K(:, 1), K(:, 2), 'Color', '#D95319')
title('Frenet (blue) and Rotated Frame (orange) curvatures')
axis equal
figure
hold on
plot(n, 'Color', '#0072BD')
plot(N, 'Color', '#D95319')
title('Frenet (blue) and Rotated Frame (orange) normals')

Cite As

Moreno, M. (2022). Discrete Curvature, Normals and Evolute (https://www.mathworks.com/matlabcentral/fileexchange/107929-discrete-curvature-normals-and-evolute), MATLAB Central File Exchange. Retrieved .

MATLAB Release Compatibility
Created with R2022a
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

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

Start Hunting!