[Math] Discrete points curvature analysis

Question

I have a set of discrete points in the x-y domain and I need to find the points of abrupt change and gradual change as well.

Approach 1:

I tried finding the curvature formed by 3 points using this equation:

where a, b, c stands for distance between all three line segments.

This worked fine in some cases with soft curve, but failed to detect some points with very less acute angle as attached in the sample output:

I have selected a threshold which segregates the cases of a break point from a smooth curve. If I lower that threshold value from the current value then it may report some false cases as well. The threshold chosen is 5 and the value of curvature at the sharp turn is 2, while for the blunt turns the curvature value is 11. I have chosen the value 5 by hit and trial.

Approach 2:

I read about the double derivative of discrete points from this blog and tested it on my data-set, But it is also not detecting that sharp edge, instead it detects seemingly vertical points as break point as:

Approach 3:

I calculated the angle subtended by 3 points at the centre and then threshold it to detect the break points. It works fine in case of sharp edges but fails in case of blunt points.

Proposed solution:

From all of the above learning implementations I concluded that I can use Approach 3 to detect sharp edges and then use Approach 1 to detect blunt edges. But it would be appreciable if get this done more efficiently in a single pass. I am open to any other solution that you may think would work in this case.

Sample cases:

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]

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]

Answer 1

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.