Maximum perpendicular distance between lines
14 次查看(过去 30 天)
显示 更早的评论
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 个评论
John D'Errico
2023-2-1
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.)
采纳的回答
Torsten
2023-2-1
编辑:Torsten
2023-2-1
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)
Pstart = [x(1)+lambda*(x(end)-x(1)),y(1)+lambda*(y(end)-y(1))]
Pend = compute_intersection(Pstart,x,y,fit)
distance = norm(Pend-Pstart)
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
更多回答(2 个)
KSSV
2023-2-1
编辑:Torsten
2023-2-1
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)
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
Image Analyst
2023-2-1
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 个评论
Image Analyst
2023-2-3
@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)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!