[Math] Curvature of a 2D discrete curve

differential-geometrydiscrete geometry

Suppose one has a discrete plane curve in 2D, such that the curve is composed entirely of connected straight line segments. Clearly, the usual definition of curvature does not apply, as this curve is not smooth.

However, is there a curvature definition one could use, that would converge to the usual definition as more and more line segments were used to describe a smooth curve? That is, for a discrete plane curve, at each vertex point, our discrete curvature should be defined such that at the limit of the smooth curve, the new curvature at the vertex, should be the curvature of the curve at that point.

Is such a definition possible?
If it is, are there standard definitions for this discrete curvature?
Is more than a single definition even possible?

$\qquad\qquad\qquad\qquad$ enter image description here

Best Answer

Heuristically what I would do is, for each polygonal approximation:

Decide on a value for both the direction and a curvature at each vertex, by fitting a circle through it and its neighbor vertices on each side.
Along each of the straight segments, assign a curvature as a quadratic function of arc length, such that (a) it agrees with the values a the endpoints that you have already decided on, and (b) its integral over the entire segment is exactly the different in assumed direction between the vertices.

If your vertices always lie on the original smooth curve, as is the case in your illustration, this ought to converge nicely (that is, locally uniformly) towards the true curvature.

For some applications it may work just as well simply to interpolate linearly in step 2. Using a quadratic correction has the advantage that you can integrate it to find a continuously varying direction at each point on the segment, which can be used to simulate shading, bouncing off the curve, and the like.

Related Solutions

[Math] Discrete points curvature analysis

I am not sure if I understand your question correctly, but I think you can use the angle, i.e., approach 3. You can make use of the inner product.

Consider point $\textbf{z}(k)=(x(k), y(k))$ and its neighbors $\textbf{z}(k-1)=(x(k-1), y(k-1))$ and $\textbf{z}(k+1)=(x(k+1), y(k+1))$. Then, the angle between the line through $\textbf{z}(k-1)$ and $\textbf{z}(k)$ and the line through $\textbf{z}(k+1)$ and $\textbf{z}(k)$ has the following value: $$ \alpha=\arccos \left( \frac{(\textbf{z}(k-1)-\textbf{z}(k))\cdot(\textbf{z}(k+1)-\textbf{z}(k))}{\| \textbf{z}(k-1)-\textbf{z}(k) \| \|\textbf{z}(k+1)-\textbf{z}(k)\|} \right). $$ If you write this out, you would get $$ \cos{\alpha}=\frac{(x(k-1)-x(k))(x(k+1)-x(k))+(y(k-1)-y(k)(y(k+1)-y(k))}{\sqrt{(x(k-1)-x(k))^2+(y(k-1)-y(k))^2}\sqrt{(x(k+1)-x(k))^2+(y(k+1)-y(k))^2}} $$

In Matlab, it may look like this:

function corners(x, y, threshold)

% Compute angle between consecutive lines
cosangle = ((x(1:end-2) - x(2:end-1)) .* (x(3:end) - x(2:end-1)) + ...
    (y(1:end-2) - y(2:end-1)) .* (y(3:end) - y(2:end-1))) ./ ...
    hypot(x(1:end-2) - x(2:end-1), y(1:end-2) - y(2:end-1)) ./ ...
    hypot(x(3:end) - x(2:end-1), y(3:end) - y(2:end-1))
angle = acos(cosangle);

% Logical vector with emphasize on sharp corners (or endpoints)
emph = [true, angle < threshold, true];

% Make plot
figure();
hold on;
plot(x, y, 'g', 'LineWidth', 3);
plot(x(emph), y(emph), 'k.', 'MarkerSize', 30)
plot(x, y, 'r.', 'MarkerSize', 20);
axis equal;

Then the code

% Case 1
x_coords = [153, 168, 186, 202, 214, 234, 248, 270, 286, 311, 327, 350, 363, 386, 398, 413, 428, 415, 394, 378, 353, 335, 307, 290, 264, 248, 226, 211, 188, 174, 157, 145, 129, 116, 100, 90, 98, 114, 139, 152, 167, 178, 195, 207, 221, 234, 246, 266, 280, 290, 278, 266, 250, 232, 214, 196, 180, 166];
y_coords = [120, 121, 125, 128, 129, 131, 131, 131, 129, 128, 127, 128, 129, 131, 133, 137, 139, 137, 135, 133, 133, 133, 135, 137, 140, 142, 145, 147, 148, 147, 146, 144, 141, 138, 132, 125, 115, 111, 107, 103, 100, 96, 91, 87, 81, 75, 69, 58, 50, 41, 34, 30, 28, 27, 29, 34, 38, 41];
corners(x_coords, y_coords, 2*pi/3);
saveas(gcf, 'example1', 'png');

% Case 2
x_coords = [387, 376, 355, 335, 321, 301, 288, 268, 255, 235, 222, 200, 187, 168, 156, 139, 127, 111, 99, 117, 130, 145, 158, 179, 193, 217, 235, 263, 283, 315, 336, 368, 389, 420, 440, 467, 483, 506, 520, 537, 549, 564, 566, 564, 562, 560, 558, 556, 554, 551, 540, 536, 532, 530, 528, 527, 526, 526, 527, 528, 531, 532, 534, 534, 534, 530, 524, 514, 503, 494, 486, 474, 465, 452, 444, 433, 423, 409, 396];
y_coords = [45, 41, 43, 51, 58, 71, 78, 89, 94, 101, 105, 109, 110, 112, 112, 112, 111, 111, 112, 122, 125, 127, 128, 129, 130, 130, 131, 131, 133, 137, 139, 144, 146, 150, 151, 151, 150, 148, 146, 143, 141, 137, 116, 101, 87, 71, 58, 42, 30, 15, 19, 30, 53, 68, 91, 105, 126, 138, 151, 166, 182, 196, 211, 224, 236, 256, 268, 276, 259, 240, 225, 203, 189, 171, 159, 147, 139, 130, 129];
corners(x_coords, y_coords, 2*pi/3);
saveas(gcf, 'example2', 'png');

generates the following figures.

Total curvature of a closed polygonal curve

Let the vertices $\ v_1, v_2, \dots, v_n\ $ lie at the points $\ z_1, z_2, \dots, z_n\ $ in the complex plane, $\ z_0=z_n\ $, $\ z_{n+1}=z_1\ $, and $\ z_j-z_{j-1}= r_je^{i\theta_j}\ $ for $\ j=1,2, \dots, n+1\ $ with $\ r_j > 0\ $ and $\ 0\le \theta_j < 2\pi\ $. If $\ \psi_j\ $ is the external angle at $\ v_j\ $ then $\ \psi_j \equiv \theta_{j+1}-\theta_j\ \left(\hspace{-0.5em}\mod 2\pi\right)$, and \begin{eqnarray} \tau&=&\sum_\limits{j=1}^n\psi_j\\ &\equiv&\sum_\limits{j=1}^n \left(\theta_{j+1}-\theta_j\right)\left(\hspace{-1em}\mod 2\pi\right)\\ &\equiv& \theta_{n+1}-\theta_1\left(\hspace{-1em}\mod 2\pi\right)\ . \end{eqnarray}

But $\ r_{n+1}e^{i\theta_{n+1}}=z_{n+1}-z_n=z_1-z_0=r_1e^{i\theta_1}\ $, so $\ \theta_{n+1}\equiv\theta_1\left(\hspace{-0.5em}\mod 2\pi\right)\ $ and $\ \tau\equiv 0\left(\hspace{-0.5em}\mod 2\pi\right)\ $.

Best Answer

Related Solutions

[Math] Discrete points curvature analysis

Total curvature of a closed polygonal curve

Related Question