How to find the minimum distance between two fitted curves

I have the following problem:

I use the fit commands

curve1 = fit(x1,y1,'smoothspline')
curve2= fit(x1,y2,'smoothspline')

Now I obtain a plot like the attached image. How do I find the minumum difference between the two curves? I can't use

[M,I] = min(abs(curve1-curve2));  

because matlab says: Undefined operator '-' for input arguments of type 'cfit'

If you see the attached image I want the min distance between the green and the red line.

[Image Analyst]

Can you attach a .mat file with all your x and y arrays?

[Raffaele Russo]

this is my script with all plots, here my aim is to find the minimum vertical distance between f5 and f6

Answer 1

[John D'Errico]

Not impossible to do. I can think of two ways to do this. You give no data, so let me make some up.

x1 = sort(rand(100,1))-.5;
y1 = sin(6*x1)/2 + randn(size(x1))/1000;
x2 = sort(rand(100,1))-.5;
y2 = cos(x2) + randn(size(x2))/1000;
curve1 = fit(x1,y1,'smoothingspline');
curve2 = fit(x2,y2,'smoothingspline');

plot(curve1)
hold on
plot(curve2)

Note that fit requires 'smoothingspline'. 'smoothspline' will just give you an error.

As I said, the closest point between the two curves can be found in several ways. Simple is to just sample a zillion points on the two curves, then find the one which minimizes the Euclidean distance. Thus...

nint = 1000;
xint = linspace(-.5,.5,nint)';
y1int = curve1(xint);
y2int = curve2(xint);

% interpoint distance matrix:
% Note: this uses a feature of R2016b or later.
% Older releases will find bsxfun useful.
d = sqrt((xint - xint.').^2 + (y1int - y2int.').^2);
d = d + diag(repmat(inf,nint,1));
[mind,ind] = min(d(:));
[i1,i2] = ind2sub([nint,nint],ind);

plot([xint(i1),xint(i2)],[y1int(i1),y2int(i2)],'-bo')
axis equal

Note my use of axis equal in there, to force the axes to have the same units. That makes it clear this really is the pair of points with minimium distance between the curves.

Another point worth noting, is IF you really wanted to just find the point with minimum vertical distance between the curves, it would have been far easier to compute. But your statement said the minimum distance between the curves.

Can you do this more exactly? Well, yes. But is the exact solution it really, truly important? That is especially true, when those curves are just smoothing splines.

How would I find the exact solution?

For that, I would need to extract each polynomial segment from the splines. Then compute the minimum distance between each pair of those segments. This reduces to essentially a polynomial rootfinding problem, although the polynomials involved will be of sufficiently high order that a numerical solver will be needed. I seriously doubt you need that much accuracy though.

[Raffaele Russo]

this is the code, i need the minimum vertical distance between f5 and f6

close all; clc; clear all;

global pch A_noz_ratio t_stay

%% Propulsion system spect

pch = 3;

A_noz_ratio = 200;

t_stay = 0.067;

Thrust = 0.15;

filename = 'Isp - Tch.csv'; %dati estrapolati dal grafico di Zakirov mediante sito

M = csvread(filename);

Tchamber = M(:,1);

Isp = M(:,2);

figure(1)

f1= fit(Tchamber,Isp,'smoothingspline'); %fit dei dati estrapolati

plot(f1,Tchamber,Isp);

xlabel('Tch'); ylabel('Isp');

h = 196;

Tch_LF = csvread('t - mp.csv'); %estrapolazione dati dal grafico di Zakirov

mflowrate = Tch_LF(:,1);

Tch = Tch_LF(:,2);

f2= fit(mflowrate,Tch,'smoothingspline');

figure(2)

plot(f2,mflowrate,Tch);

xlabel('mflowrate'); ylabel('Tch');

LF = mflowrate*(1e-3/((0.25*pi)*0.013^2)); %calcolo l'LF corrispondente ai valori estrapolati

LF_sperimentali = [0.61;0.79;1.43;2.04]; %sapendo che Z ha usato un catalizzatore a Phi=13mm

Tch_sperimentali = [1050;1195;1057;433];

Tch = [Tch;Tch_sperimentali]; %unisco i dati sperimentali ottenuti da prove in laboratorio con un

LF = [LF;LF_sperimentali]; %sistema propulsivo con le seguenti caratteristiche(inserire dati Peppe)

figure(3);

f3 = fit(LF,Tch,'smoothingspline'); %fitto Zakirov + dati sperimentali

plot(f3,LF,Tch);

xlabel('LF'); ylabel('Tch');

figure(4);

Isp_vec = f1(Tch); %ricavo l'impulso specifico come funzione della temperatura alle Tch da me conosciute

f4 = fit(LF,Isp_vec*1e-4,'smoothingspline'); %lo fitto ed effettuo una riduzione di scala di 1:10000

plot(f4,LF,Isp_vec*1e-4);

xlabel('LF'); ylabel('Isp/10000');

H = csvread('Hr.csv'); %estrazione di dati su entalpia di formazione a Tstandard e T qualsiasi

T = H(:,1);

Hr_std = H(:,2);

DHr = H(:,3);

DHr_Tdec= interp1(T,DHr,Tch+273,'linear','extrap'); %estrapolo il delta entralpia avendo interpolato linearmente i dati in tabella

Hr_Tdec = (82.1+DHr_Tdec)/0.044013; %calcolo Hr alle Tch in camera

% a questo punto posso calcolare il diametro minimo con la formula ...

... di Zakirov essendo noti Hr e Tch

dchamber_min = t_stay*(4*h*(8.31441/0.044013).*(Tch+273).*((Tch+273)-298))./(3e5.*Hr_Tdec*1e3);

figure(5)

f5 = fit(LF,dchamber_min,'smoothingspline');

plot(f5,LF,dchamber_min);

xlabel('LF'); ylabel('Dmin');

mp = Thrust./(Isp_vec*9.81); %calcolo il mass flow rate ai vari Isp

Area = mp./LF; %calcolo le aree corrispondenti agli LF considerati

D = ((4/pi).*Area).^0.5; %posso calcolare il diametro banalmente

figure(6)

f6 = fit(LF,D,'smoothingspline');

plot(f6,LF,D);

xlabel('LF'); ylabel('Dmin');

figure(7)

plot(f4,'m'); hold on

plot(LF,Isp_vec*1e-4,'m.');

hold on

plot(f5,'r'); hold on

plot(LF,dchamber_min,'.r'); hold on

plot(f6,'g');

plot(LF,D,'.g');

hold on

axis([0 max(LF) 0 0.025]);

legend('Isp','Isp data','Dmin','dmin data','D ','D data');

xlabel('LF');

[Raffaele Russo]

unfortunately i have a problem with your solution, when i apply it to my script i can find only one point of minimum distance but because my curves have 2 intersection points the points are two and i need to know both (for me would be good to have the ascissa of minimum distance as output of the function).

[John D'Errico]

I hate moving targets. First minimum distance. Then we are told it is vertical distance you wanted. But in reality, maybe you want the intersection points, and you want to find them all.

So, which is it? Do you want to find the minimum vertical distance between f5 and g6? Or the intersection points of f5 and f6? I'll assume that now you mean to look for intersection points.

Use Doug Schwarz's utility: intersections. It is found on the file exchange for download.

Since I still have seen no data from you, you will need to do it yourself.

xint = linspace(min(LF),max(LF),1000);
y5int = f5(xint);
y6int = f6(xint);
[X0,Y0] = intersections(xint,y5int,xint,y6int);

https://www.mathworks.com/matlabcentral/fileexchange/11837-fast-and-robust-curve-intersections