Find specified number of points on curve with identical chord length (not arc length)

Hello!

I have the user draw a line and obtain the points on that line:

% Load example image
I = imread('peppers.png');
imshow(I)
% Draw line
roi = drawassisted('Closed', false);
% Obtain positions of points
points_input = roi.Position; 


% Number of target points
num_points_target = 10;

% Find these
points_output

I want to find a specified number of points (num_points_target) along the described curve, so that the chord length between the outputs points (points_output) would be identical. I am specifically looking for chord length, not arc length. I have seen this but it uses arc length ( https://uk.mathworks.com/matlabcentral/fileexchange/34874-interparc/files/interparc.m ).

I assume some interpolation - linear would likely be okay - would be in order as, depending on the drawing speed, the points_input might not be spaced very evenly.

Code would obviously be great but if someone felt like just explaining the mathematical approach to this, that would already be very helpful. At the moment, I think I would be able to cobble something dreadful with many loops together myself but it seems to me that there must be an elegant solution to this problem surely.

[Image Analyst]

There may be more than one points along the drawn curve that are the same distance from the starting point. Are you planning on handling that situation? Just take the first one in the list of (x,y) pixel locations that they drew?

[Matt J]

The ambiguity mentioned by @Image Analyst's can also lead to multiple valid solutions. Consider the path ABCD below. If the specified number of sub-chords is two, then all sequences AQD are possible solutions, where Q is any point on line segment BC.

[Julia]

That's a very valid point. I would tend towards @Image Analyst's suggestion and use the first point in the list (so the first drawn point) if there are multiple solutions.

[Julia]

Can I say thank you for all your thoughts and comments so far - I really appreciate it! :)

[Matt J]

I would tend towards @Image Analyst's suggestion and use the first point in the list (so the first drawn point) if there are multiple solutions.

No matter what rule you choose, the solution will be a discontinuous function of the input points, meaning that the results can vary unstably across different computers and CPUs. For example, using the "first point" rule, the first point in my original configuration above that gives you a solution is B, and the path you would choose is ABD. However, if I were to rotate BC by even the smallest epsilon>0 degrees, either clockwise or counterclockwise (see below), then now the only solution is the path through the center, AOD. Not only that, but the chord length drops discontinuously from AB to the circle radius AO.

Answer 1

[Matt J]

Here's an approach using fminbnd to minimize the attached cost function circleMarch().

t=linspace(0,200,1e4)'; %%%
XY=[cosd(t), sind(t)];  %Example input data
Nsteps=5;               %%%

fun=@(d)circleMarch(XY,d,Nsteps);

tic;
    interDists=vecnorm(diff(XY),2,2);
    Lmax=sum(interDists);
    Lmin=min(interDists);
    d=fminbnd(fun, Lmin,Lmax, optimset('MaxIter',inf,'TolX',1e-6));
toc

Elapsed time is 0.125987 seconds.

[res,~,xy]=fun(d); xy(end,:)=XY(end,:)

res = 8.6503e-07

xy = 6×2

1.0000 0 0.7660 0.6428 0.1736 0.9848 -0.5000 0.8660 -0.9397 0.3420 -0.9397 -0.3420

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

chordLengths=vecnorm(diff(xy),2,2),

chordLengths = 5×1

0.6840 0.6840 0.6840 0.6840 0.6840

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

plot(XY(:,1),XY(:,2),'.', xy(:,1),xy(:,2),'r'); axis equal; axis padded
legend 'Original Curve' 'Chords'

Warning: Hardware-accelerated graphics is unavailable. Displaying fewer markers to preserve interactivity.

[Julia]

This looks very exciting even though I must admit I have no clue what is going on.

I tried it with the input data below and got the following error message:

Error using circleMarch (line 62)
Final segment has zero length.

Error in @(d)circleMarch(XY,d,Nsteps)

Error in fminbnd (line 238)
fx = funfcn(x,varargin{:});
     ^^^^^^^^^^^^^^^^^^^^^

XY = [94	310
94	308
94	300
86	254
86	212
86	182
92	164
108	152
130	150
144	162
158	188
172	214
180	224
186	232
210	234
216	234
224	230
230	210
230	202
242	200
252	200
272	216
282	232
286	234
304	236
348	214
348	208
350	186
340	164
326	150
318	138
318	112
332	98
356	92
396	90
408	90
410	90
410	90];

[Matt J]

It is because your final two points are the same. You need to weed out duplicates (but keep the points in order, see unique(___stable))

XY = [94	310
94	308
94	300
86	254
86	212
86	182
92	164
108	152
130	150
144	162
158	188
172	214
180	224
186	232
210	234
216	234
224	230
230	210
230	202
242	200
252	200
272	216
282	232
286	234
304	236
348	214
348	208
350	186
340	164
326	150
318	138
318	112
332	98
356	92
396	90
408	90
410	90];


Nsteps=15;

fun=@(d)circleMarch(XY,d,Nsteps);

tic;
    interDists=vecnorm(diff(XY),2,2);
    Lmax=sum(interDists);
    Lmin=min(interDists);
    d=fminbnd(fun, Lmin,Lmax, optimset('MaxIter',inf,'TolX',1e-6));
toc

Elapsed time is 0.106696 seconds.

[res,~,xy]=fun(d); xy(end,:)=XY(end,:)

xy = 16×2

94.0000 310.0000 88.0935 266.0376 86.0000 221.7296 87.5338 177.3987 122.9142 150.6442 154.6001 181.6859 176.8483 220.0603 219.4976 232.2512 249.9515 200.0000 281.4928 231.1885 325.4540 225.2730 349.8005 188.1943 326.4688 150.4688 323.8088 106.1912 365.6685 91.5166

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

chordLengths=vecnorm(diff(xy),2,2),

chordLengths = 15×1

44.3574 44.3574 44.3574 44.3574 44.3574 44.3574 44.3574 44.3574 44.3574 44.3574 44.3574 44.3574 44.3574 44.3574 44.3574

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

plot(XY(:,1),XY(:,2),'.', xy(:,1),xy(:,2),'r-x'); axis equal; axis padded
legend 'Original Curve' 'Chords'

[Julia]

Dear Matt,

Sorry, this was a stupid mistake to make - obviously the points list must not have replicate points... Thanks for pointing it out!

Your code works a treat and really does exactly what I hoped it would do. I understand the raised concerns about there not being a clean mathematical solution but for what I want it is good enough and I can guard against it failing simply by prompting new input (if necessary - I expect this will be rare).

In our case, this will be used to distribute a set number of evenly spaced measurement points on the user-drawn line. The line itself does not matter, it just guides the point placement. The important thing is that the points are positioned on the (interpolated) line, so experimentally in the area we want to measure. The readout will be a point grid later and thus it's important that we have equal chord spacing (as no one will later on see the line the user had drawn originally, just the measurement point locations along some anatomical structure).

Answer 2

[Image Analyst]

Use pdist2 to find the distance of every point to every other point. Then run along every row, starting with the diagonal element in each row and going to the right side of the array and use min to find the other element that is closest to your desired chord length. Pseudocode (untested):

xy = [x(:), y(:)];
distances = pdist2(xy, xy);
for row = 1 : height(xy)
    thisRow = distances(row, :);
    thisRow(1:row) = inf; % So we don't find any points earlier in the tracing.
    [min, column] = min(abs(thisRow - desiredChordLength));
    % Save the column
    nextPoint(row) = column;
end

[Julia]

Dear ImageAnalyst,

Thank you so much for your time and thoughts! I can see you have thought about this for a while :) What you suggest is also the first idea that sprang to my mind.

One practical problem might be that one would need to make sure that the data are sufficiently interpolated before starting this to avoid encountering very large steps which then require recomputing the entire matrix. This, on the other hand, means that one might create a rather large matrix in the first place. Possible but maybe not optimal.

More importantly, however, this works fine if you know what chord length step you are looking for (so chord length is the defining factor). But I am looking for a set number of steps and the problem, as I see it, is exactly that we don't know what chord length step is needed if we want to have, for example, exactly 10 steps from start to finish.

Answer 3

[John D'Errico]

Admittedly, interparc does not solve the problem you want to solve. At least, not the probem you CLAIM you want to solve. As the author, I'll admit that. The problem is, the question you have posed is not well-posed, as there are infinitely many solutions. Consider this curve:

x = [0 1 1.5 2 3];
y = [0 0 10 0 0];

plot(x,y,'o-')
grid on

The actual arc length along that curve is fairly long, at 22.025 units.

dt = sum(sqrt(diff(x).^2 + diff(y).^2))

dt = 22.0250

However, if I choose a curve with unit chordal spacing, I might get this:

x1 = [0 1 2 3]; y1 = [0 0 0 0];
hold on
plot(x1,y1,'rs-')

But I could have also chosen a curve that follows that hat shape, also with unit chordal spacing. And worse, I could have done so in various ways. This is why I claim your question is ill-posed, because it would have multiple solutions, all equally valid, yet all very different.

[Julia]

Dear John,

That is very valid and thank you for pointing it out. Please see my explanation at end of Matt's code. For practical purposes, I think we just need a valid solution - mathematically not very satisfactory, agreed, but in practical terms we are very unlikely to encounter the mentioned problematic cases frequently.