The problem of finding the geometric (minimum) distance between two arbitrary ellipses is surprisingly difficult. The general problem of finding all stationary points (minimum/maximum/saddle, no less than 12 possible points) has indeed been solved, but that algorithm is horrificly complex and would require thousands if not millions of operations once implemented.
The function distanceEllipseEllipse() is based on a somewhat more practical algorithm which limits itself to minima in the distance function. It is based on repeatedly finding the distance between a point and an ellipse, which can be done analytically (see my other post, distanceEllipsePoints.m). The algorithm by itself is not very robust (it frequently finds a local minimum, which is *not* the true distance).
This function executes the algorithm 4 times, for 4 different initial values, which greatly improves its robustness. Some numerical experimentation (comparing with a bruteforce search) has shown that the true minimum distance is returned in more than 95% of the cases. The algorithm is implemented such that MATLAB's JITaccelerator can accelerate it to the fullest extent, which makes it pretty fast and wellsuited to handle large datasets requiring this calculation.
This is an implementation of the algorithm described in
IkSung Kim: "An algorithm for finding the distance between two
% ellipses". Commun. Korean Math. Soc. 21 (2006), No.3, pp.559567.
A copypastable example (also in header of Mfile):
% Ellipse1 Ellipse2(=circle)
a = [2.0 1.0];
b = [0.5 1.0];
c = {[0,0,0], [2,2,0]}; % location of centers
u = {[1,0,0], [1,0,0]}; % both oriented in XYplane
v = {[0,1,0], [0,1,0]}; % to visualize them more easily
% plot the ellipses
f = 0:0.01:2*pi;
E1 = [a(1)*cos(f) + c{1}(1); b(1)*sin(f) + c{1}(2)];
E2 = [a(2)*cos(f) + c{2}(1); b(2)*sin(f) + c{2}(2)];
figure, hold on
plot(E1(1,:),E1(2,:),'r', E2(1,:),E2(2,:),'b')
axis equal
% run routine
[min_dist, fp_min, fs_min] = ...
distanceEllipseEllipse(a,b,c,u,v)
% plot the minimum distance returned
x = [a(1)*cos(fp_min) + c{1}(1), a(2)*cos(fs_min) + c{2}(1)];
y = [b(1)*sin(fp_min) + c{1}(2), b(2)*sin(fs_min) + c{2}(2)];
line(x,y,'color', 'k')
This should generate the screen shot given here.
If you find this work useful, please consider a small donation:
https://www.paypal.me/RodyO/3.5
1.2.0.0  Description update 

1.2.0.0  [linked to Github] 

1.1.0.0  Updated contact info 
Create scripts with code, output, and formatted text in a single executable document.
Rody Oldenhuis (view profile)
@xiaoqing Yes, indeed; I'll update the example to allow for easier tinkering with the input data :)
xiaoqing (view profile)
okay, so for the example in your .m file to correctly visualize the result, u,v need to be normalized before calculating the coordinates of those two points.
xiaoqing (view profile)
a = [2.0 1.2];
b = [0.5 1.0];
c = {[0,0,0],[1,3,0]}; % location of centers
u = {[1,1,0], [1,0,0]}; % both oriented in XYplane
v = {[1,1,0], [0,1,0]}; % to visualize them more easily
does not seem to calculate a correct minimum distance.
Cihan Ulas (view profile)
Hi,
I generated two plane with as follows;
par1.a=2;par1.b=2;par1.c=1;par1.d=2;
par2.a=10;par2.b=1;par2.c=2;par2.d=100;
P1=generate_plane(par1,10:10,10:10,'r.');
P2=generate_plane(par2,50:40,50:40,'b.');

function P=generate_plane(par,rangex,rangey,color)
[x y]=meshgrid(rangex,rangey);
z=(par.a*x(:)+par.b*y(:)+par.d)/par.c;
z=3*rand+z;
P=[x(:) y(:) z];

Then find the principal axes and orientations as follows,
[U1 E1]=svd(cov(P1));
[U2 E2]=svd(cov(P2));
a=[E1(1,1) E2(1,1)]; %principal axis
b=[E1(2,2) E2(2,2)] ; %minor axis
c={mean(P1), mean(P2)};
u= {U1(:,1)', U2(:,1)'};
v= {U1(:,2)', U2(:,2)'};
However the algorithm can not find the correct distance. Also if two ellipses are intersect it fails. The result should be 0.
For 3D, an explicit example is necessary.