# Maximum perpendicular distance between lines

13 views (last 30 days)
Michiel on 1 Feb 2023
Commented: Image Analyst on 3 Feb 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
Michiel on 1 Feb 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?

Torsten on 1 Feb 2023
Edited: Torsten on 1 Feb 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
Michiel on 2 Feb 2023
Edited: Michiel on 2 Feb 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!

KSSV on 1 Feb 2023
Edited: Torsten on 1 Feb 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 CommentsShowHide 1 older comment
Torsten on 1 Feb 2023
@KSSV forgot to set axis equal.

Image Analyst on 1 Feb 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 CommentsShowHide 1 older comment
Image Analyst on 3 Feb 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)