Borrar filtros
Borrar filtros

Maximum perpendicular distance between lines

20 visualizaciones (últimos 30 días)
Michiel
Michiel el 1 de Feb. de 2023
Comentada: Image Analyst el 3 de Feb. de 2023
I have data from a cycling lactate test and plotted a third polynomial line (black). Now I'd like to calculate the maximum perpendicular distance between the red line and the black line (which is method to determine lactate threshold).
Can anyone tell me what the best way is to calculate this? Thanks!
The variables are defined as follows:
x = [120 150 180 210 240 270 300 330 360]; % power
y= [ 1.4000 1.3000 1.3000 1.4000 1.7000 2.4000 3.9000 6.2000 13.0000]; % lactate
fit = polyfit(x,y,3);
x_fit = min(x):0.01:max(x); %x coordinates black line
y_fit = polyval(fit, x_fit); %y coordinates black line
dmax_coeff = polyfit([x(1) x(end)],[y(1) y(end)],1);
dmax_x = min(x):0.01:max(x); % x coordinates red line
dmax_y = polyval(dmax_coeff, dmax_x); % y coordinates red line
  3 comentarios
John D'Errico
John D'Errico el 1 de Feb. de 2023
I would offer an answer, but it is not at all clear what the question is.
What is a perpendicular distance?
Does it correspond to the length of a line that is perpendicular to BOTH curves? In that case, there will be exactly two such lines when one curve is a straight line, and the second is a cubic polynomial.
Is it the maximum distance to the cubic in a direction perpendicular to the straight line?
Is it the maximum VERTICAL distance between the curves? So perpendicular to the x axis? (That is what KSSV did.)
Michiel
Michiel el 1 de Feb. de 2023
Thanks, let me clarify. The line should be perpendicular to the red line. So basically draw a line that is perpendicular to the red line (for each datapoint of the red line) and measure the distance between the intersection of the red line with the black line. Of all these distances I need the maximum value. Is that clear?

Iniciar sesión para comentar.

Respuesta aceptada

Torsten
Torsten el 1 de Feb. de 2023
Editada: Torsten el 1 de Feb. de 2023
x = [120 150 180 210 240 270 300 330 360]; % power
y= [ 1.4000 1.3000 1.3000 1.4000 1.7000 2.4000 3.9000 6.2000 13.0000]; % lactate
fit = polyfit(x,y,3);
xlin = 120:360;
ylin = polyval(fit,xlin);
hold on
plot(x,y,'o')
plot(xlin,ylin)
plot([x(1),x(end)],[y(1),y(end)])
lb = 0;
ub = 1;
lambda0 = 0.5;
lambda = fmincon(@(X)obj(X,x,y,fit),lambda0,[],[],[],[],lb,ub)
Local minimum possible. Constraints satisfied. fmincon stopped because the size of the current step is less than the value of the step size tolerance and constraints are satisfied to within the value of the constraint tolerance.
lambda = 0.6534
Pstart = [x(1)+lambda*(x(end)-x(1)),y(1)+lambda*(y(end)-y(1))]
Pstart = 1×2
276.8198 8.9796
Pend = compute_intersection(Pstart,x,y,fit)
Pend = 1×2
277.1405 2.3452
distance = norm(Pend-Pstart)
distance = 6.6422
plot([Pstart(1);Pend(1)],[Pstart(2),Pend(2)])
hold off
axis equal
function v = obj(X,x,y,fit)
Pstart = [x(1)+X*(x(end)-x(1)),y(1)+X*(y(end)-y(1))];
Pend = compute_intersection(Pstart,x,y,fit);
v = -norm(-Pstart+Pend)^2;
end
function Pend = compute_intersection(Pstart,x,y,fit)
g = @(mu) Pstart + mu*[-(y(end)-y(1)),x(end)-x(1)];
h = @(u) [u,polyval(fit,u)];
fun = @(mu,u) g(mu)-h(u);
sol = fsolve(@(x)fun(x(1),x(2)),[0.5,0.5*(x(1)+x(end))],optimset('Display','off'));
Pend = h(sol(2));
end
  1 comentario
Michiel
Michiel el 2 de Feb. de 2023
Editada: Michiel el 2 de Feb. de 2023
Thanks, but I'm not sure what it does; the intersecting line which is plotted is not at all perpendicular to the (in this case) yellow line?
Edit: ah I think it had to do with axis not set to equal. When set to equal it looks perpendicular. Thanks!

Iniciar sesión para comentar.

Más respuestas (2)

KSSV
KSSV el 1 de Feb. de 2023
Editada: Torsten el 1 de Feb. de 2023
You have to use the foot of the perpendicular formula.
clc; clear all ;
x = [120 150 180 210 240 270 300 330 360]; % power
y= [ 1.4000 1.3000 1.3000 1.4000 1.7000 2.4000 3.9000 6.2000 13.0000]; % lactate
fit = polyfit(x,y,3);
x_fit = min(x):0.01:max(x); %x coordinates black line
y_fit = polyval(fit, x_fit); %y coordinates black line
dmax_coeff = polyfit([x(1) x(end)],[y(1) y(end)],1);
dmax_x = min(x):0.01:max(x); % x coordinates red line
dmax_y = polyval(dmax_coeff, dmax_x); % y coordinates red line
% GEt line equation (REd line)
p = polyfit(dmax_x,dmax_y,1) ;
m = p(1) ; c = p(2) ;
a = m ; b = -1 ; % converting y = mx+c into ax+by+c = 0
% Get foot of the perendiculars form each blacl line to the red line
% Formula: (h-p)/a = (k-q)/b =-(ap+bq+c)/(a2+b2)
p = x_fit ; q = y_fit ;
h = -(a*p+b*q+c)/(a^2+b^2)*a+p ;
k = -(a*p+b*q+c)/(a^2+b^2)*b+q ;
% Get perpendicular distance
d = sqrt((p-h).^2+(q-k).^2) ;
% Pick max distance
[val,idx] = max(d)
val = 6.6422
idx = 15715
plot(dmax_x,dmax_y,'r',x_fit,y_fit,'k')
hold on
plot(p(idx),q(idx),'*r',h(idx),k(idx),'*r')
plot([p(idx) h(idx)],[q(idx) k(idx)],'r')
axis equal
  2 comentarios
Michiel
Michiel el 1 de Feb. de 2023
Thanks, but that's not what I meant. See my explanation with John D'Errico.
Torsten
Torsten el 1 de Feb. de 2023
@KSSV forgot to set axis equal.

Iniciar sesión para comentar.


Image Analyst
Image Analyst el 1 de Feb. de 2023
What you describe has a name. It's called the "triangle threshold method". It's fairly well known in the Image Processing community and is good for finding the "corners" in skewed curves like you show in your example. It's in ImageJ but I don't think it's in MATLAB yet so I wrote my own. It's especially good for finding thresholds for images where the histogram shape is a skewed Gaussian.
I'm attaching the function. It will work for your data.
  2 comentarios
Michiel
Michiel el 2 de Feb. de 2023
Thanks, I had a look, but I'm not sure how I should use this function on my data?
Image Analyst
Image Analyst el 3 de Feb. de 2023
@Michiel I made a few changes to make the graphics better for non-histogram/line plot data like yours. I'm attaching the new version. Now I also show the distance to the line for every data point with the one in red being the longest distance.
Test code to assign your data and call the function:
% Create data
x = [120 150 180 210 240 270 300 330 360]; % power
y= [ 1.4000 1.3000 1.3000 1.4000 1.7000 2.4000 3.9000 6.2000 13.0000]; % lactate
plot(x, y, 'b.-', 'LineWidth', 2, 'MarkerSize', 30)
grid on;
% Get the index of the data point with the longest distance to the hypoteneuse.
index = triangle_threshold(y, 'L', true)
% Put lines on plot showing the result.
xline(x(index), 'Color','r', 'LineWidth',2)
yline(y(index), 'Color','r', 'LineWidth',2)

Iniciar sesión para comentar.

Categorías

Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by