Maximum perpendicular distance between lines

조회 수: 13 (최근 30일)
Michiel
Michiel 2023년 2월 1일
댓글: Image Analyst 2023년 2월 3일
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
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.)
Michiel
Michiel 2023년 2월 1일
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?

댓글을 달려면 로그인하십시오.

채택된 답변

Torsten
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)
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
Michiel
Michiel 2023년 2월 2일
편집: Michiel 2023년 2월 2일
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!

댓글을 달려면 로그인하십시오.

추가 답변 (2개)

KSSV
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)
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
Michiel
Michiel 2023년 2월 1일
Thanks, but that's not what I meant. See my explanation with John D'Errico.
Torsten
Torsten 2023년 2월 1일
@KSSV forgot to set axis equal.

댓글을 달려면 로그인하십시오.


Image Analyst
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
Michiel
Michiel 2023년 2월 2일
Thanks, I had a look, but I'm not sure how I should use this function on my data?
Image Analyst
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)

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Solver Outputs and Iterative Display에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by