Maximum perpendicular distance between lines

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

[Torsten]

To which curve should the distance line be perpendicular ? To the red line or to the fitted data curve ? Or to the horizontal line (vertical distance) ?

[John D'Errico]

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]

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?

Answer 1

[Torsten]

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]

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!

Answer 2

[KSSV]

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

[Michiel]

Thanks, but that's not what I meant. See my explanation with John D'Errico.

[Torsten]

@KSSV forgot to set axis equal.

Answer 3

[Image Analyst]

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.

[Michiel]

Thanks, I had a look, but I'm not sure how I should use this function on my data?

[Image Analyst]

@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)