Roundest ellipse with specified tangents Suppose you are given two points $\mathbf p_1$ and $\mathbf p_2$ in the plane with associated vectors $\mathbf v_1$ and $\mathbf v_2$. You want to find an ellipse passing through $\mathbf p_1$ tangent to $\mathbf v_1$ and through $\mathbf p_2$ tangent to $\mathbf v_2$. There are infinitely many solutions,


but a natural choice is to take the one that is the closest to a circle, i.e. the one with the minimum eccentricity. Is there an elegant way to compute such an ellipse given $(\mathbf p_1, \mathbf v_1)$ and $(\mathbf p_2, \mathbf v_2)$?
 A: I found a way to do this construction in a purely geometrical way, without coordinates. Let $A$, $B$ be the given tangency points (called $\mathbf p_1$ and $\mathbf p_2$ in the question) and $AV$, $BV$ the corresponding tangents (the case when tangents are parallel can be dealt with in a simpler way, see below).
If $M$ is the midpoint of $AB$, for a well-known property of the ellipse line $VM$ is a diameter, and we are free to choose on it a point $D$, to construct the ellipse touching the given tangents at $A$, $B$ and passing through $D$. The construction is quite simple, because the line through $D$ parallel to $AB$, meeting line $VA$ at $E$, is another tangent: if $F$ is the midpoint of $AD$ it follows that $EF$ is also a diameter, and the intersection $C$ of $VM$ and $EF$ is thus the center of the ellipse (see diagram). With the center and three points, the ellipse can be easily constructed.

But of course we must choose point $D$ such that the eccentricity of the ellipse be minimum. To that end, draw line $GC$ parallel to $AB$ and notice that lines $CD$ and $CG$ form a pair of conjugate diameters. Apollonius' equations connect semidiameters $CD$ and $CG$ with standard semiaxes $a$ and $b$ of the ellipse:
$$
CD^2+CG^2=a^2+b^2,\quad CD\cdot CG\sin\theta = ab,
$$
where $\theta=\angle VMA$. From there it is not difficult to find an expression for the eccentricity $\epsilon$ as a function of $t=CD/CG$:
$$
\epsilon^2={2\over1+1/\xi},
\quad\hbox{where}\quad
\xi=\sqrt{1-\left({2t\over1+t^2}\right)^2\sin^2\theta}.
$$
It follows that

the minimum value of the eccentricity, 
  $\epsilon=\sqrt{2\over1+1/|\cos\theta|}$, is obtained when
  $CD=CG$.

To exploit that condition, notice that $(AM/CG)^2+(CM/CD)^2=1$ and if $CD=CG$ this is equivalent to
$$
AM^2+CM^2=CD^2.
$$
Playing around with similar triangles one also finds $CD/MC=1+VM/MD$: combining those equalities we finally get
$$
MD={AM^2+AM\sqrt{AM^2+VM^2}\over VM}.
$$
From that, one can construct point $D$ and then find center $C$ as explained above.
If tangent lines are parallel, $A$ and $B$ are the endpoints of a diameter and center $C$ is their midpoint. Just construct then the conjugate diameter through $C$ parallel to the tangents and take on it point $D$ such that $CD=AC$.
A: At the moment I can offer an analytical solution, which is quite cumbersome but manageable. I hope someone will be able to give a simpler solution.
To reduce the number of parameters, we can set up a coordinate system such that $\mathbf{p}_1=(\alpha,0)$, $\mathbf{p}_2=(0,\beta)$, $\mathbf{v}_1=(0,1)$, and define then $m=v_{2y}/v_{2x}$.
Let's start with the generic equation for a conic section, where we can set to $1$ the coefficient of $x^2$ because we are dealing with an ellipse:
$$
x^2+By^2+Cxy+Dx+Ey+F=0.
$$
We have four conditions: the conic passes through $\mathbf{p}_1$, $\mathbf{p}_2$ and is there tangent to $\mathbf{v}_1$, $\mathbf{v}_2$. These conditions translate into four equations:
$$
\begin{align}
\cases{
\alpha^2+\alpha D+F=0\\
\beta^2B+\beta E+F=0\\
2m\beta B+\beta C+D+mE=0\\
\alpha C+E=0
}
\end{align}
$$
From there, one can find expressions for $B$, $C$, $D$, $E$ as a function of $F$. Plugging these into the formulas for the semiaxes $a$ and $b$ we can then compute the eccentricity $\epsilon=\sqrt{1-b^2/a^2}$. The smaller the eccentricity, the rounder the ellipse, so we must find the minimum of $\epsilon$ as a function of $F$.
I computed $d\epsilon/dF$ with Mathematica and found that it vanishes for
$$
F=-\frac{\alpha ^2 \beta  \left(\alpha ^2 \beta +\beta ^3-2 \alpha ^2 \beta  m^2+2 \alpha ^3
   m\right)}{\beta ^2 \left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta
   ^2\right)+2 \alpha  \beta  m \left(\alpha ^2+2 \beta ^2\right)}.
$$
Once $F$ is known, one can compute the other coefficients and find the equation of the "roundest" ellipse:
$$
\begin{align}
\cases{
B=\frac{\displaystyle\alpha ^2 \left(\beta ^2+\alpha ^2 \left(2 m^2+1\right)+2 \alpha  \beta  m\right)}{\displaystyle\beta ^2\left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha 
   \beta  m \left(\alpha ^2+2 \beta ^2\right)}\\
\\
C=\frac{\displaystyle2 \alpha ^2 m \left(-\alpha ^2+\beta ^2+2 \alpha  \beta  m\right)}{\displaystyle\beta ^2 \left(\alpha
   ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha  \beta  m
   \left(\alpha ^2+2 \beta ^2\right)}\\
\\
D=-\frac{\displaystyle2 \alpha ^2 m \left(2 \beta ^3+m \left(\alpha ^3+3 \alpha  \beta ^2\right)\right)}{\displaystyle\beta ^2
   \left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha 
   \beta  m \left(\alpha ^2+2 \beta ^2\right)}\\
\\
E=-\frac{\displaystyle2 \alpha ^3 m \left(-\alpha ^2+\beta ^2+2 \alpha  \beta  m\right)}{\displaystyle\beta ^2 \left(\alpha
   ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha  \beta  m
   \left(\alpha ^2+2 \beta ^2\right)}
}
\end{align}
$$
A check with GeoGebra confirms that this is indeed the value of $F$ giving a minimum eccentricity. 
EDIT.
Following the suggestion by Rahul, we can get a much simpler result. Choose a coordinate system such that 
$\mathbf{p}_1=(\alpha,0)$, $\mathbf{p}_2=(-\alpha,0)$, 
and define $m=v_{1y}/v_{1x}$, $n=v_{2y}/v_{2x}$.
With those choices, the equation of the ellipse can be written as:
$$
x^2+By^2-{m+n\over mn}xy+\alpha{m-n\over mn}y-\alpha^2=0.
$$
The eccentricity $\epsilon$ is then a function of $B$, and 
$d\epsilon/dB=0$ for
$$
B=1+\frac{(m+n)^2}{2 m^2 n^2}.
$$
In summary, the equation of the roundest ellipse turns out to be:
$$
x^2+\left(1+\frac{(m+n)^2}{2 m^2 n^2}\right)y^2-{m+n\over mn}xy+\alpha{m-n\over mn}y-\alpha^2=0.
$$

Added by Rahul:
Things become even simpler if we use the slope of the normals, $\mu=-1/m$ and $\nu=-1/n$:
$$
x^2+\left(1+\tfrac12(\mu+\nu)^2\right)y^2+(\mu+\nu)xy+\alpha(\mu-\nu)y-\alpha^2=0.
$$
A: In an answer to a previous not-quite-duplicate question, achille hui has given an elegant coordinate-free solution, which I summarize below.
Choose the origin at the intersection of the two desired tangent lines, so that $\mathbf p_1$ and $\mathbf p_2$ are interpreted as vectors from the intersection to the two input points. Let $\mathbf q_1,\mathbf q_2$ be the dual basis for $\mathbf p_1,\mathbf p_2$, that is, $\mathbf p_1\cdot \mathbf q_1 = \mathbf p_2\cdot \mathbf q_2 = 1$ and $\mathbf p_1\cdot \mathbf q_2 = \mathbf p_2\cdot \mathbf q_1 = 0$. The equation for the minimum-eccentricity ellipse is
$$(\mathbf{q}_1\cdot \mathbf{x} - 1)^2
+ (\mathbf{q}_2\cdot \mathbf{x} - 1)^2
+ 2\alpha(\mathbf{q}_1\cdot \mathbf{x})(\mathbf{q}_2\cdot \mathbf{x}) = 1,$$
where
$$\alpha = \frac{2\mathbf{p}_1\cdot\mathbf{p}_2}{\|\mathbf{p}_1\|^2 + \|\mathbf{p}_2\|^2}.$$
For the derivation, see the original post.
A: Least eccentric ellipses for geometric Hermite interpolation
John C. Femiani, Chia-Yuan Chuang, Anshuman Razdan
Computer Aided Geometric Design 29 (2012) 141–149
A: I needed to implement a solution to this problem and wound up frustrated with the existing StackExchange answers. I can't say whether a given answer is right or wrong; only that I tried reducing a few of the more straightforward answers to practice and couldn't get them working. 
I had more luck with the combination of (the already mentioned) "Least eccentric ellipses for geometric Hermite interpolation" by Femiani, Chuang, and Razdan (Computer Aided Geometric Design 29, 2012 pp. 141-9) and "Characteristics of conic segments in Bezier form" by Javier Sánchez-Reyes (Proceedings of the IMProVe 2011, pp. 231-4).
Below is an Octave script making use of the calculations of these papers. I calculate x1 and y1 from the first paper by different means (and they may sometimes have different signs), but the equivalence should be easy to verify. The complex-number calculations from the second paper can be easily implemented in C using complex.h. I include an alternate version of projPointOnLine() for reference. The lines below the ang calculation present a figure with the calculated ellipse. 
pkg load symbolic geometry

function rp = ppol(l1, l2, p)
  if (l2(1) == l1(1))
    rp(1) = l2(1);
    rp(2) = p(2);
  else
    m = (l2(2) - l1(2)) / (l2(1) - l1(1));
    b = l1(2) - m * l1(1);
    rp(1) = (m * p(2) + p(1) - m * b) / (m * m + 1);
    rp(2) = (m * m * p(2) + m * p(1) + b) / (m * m + 1);
  endif
endfunction

# Points and slopes
p0 = [ -4, 3 ];
h0 = [ -.8, .3 ];
p2 = [ 2, -1 ];
h2 = [ 2.2, -3 ];

# "Least eccentric ellipses for geometric Hermite interpolation"
# Femiani, Chuang, and Razdan
# Computer Aided Geometric Design 29, 2012 pp. 141-9
# (x1 and y1 calculated by different means)
l0 = createLine([ p0, h0 ]);
l2 = createLine([ p2, h2 ]);
p1 = intersectLines(l0, l2)

l3 = createLine(p0, p2);
p1p = projPointOnLine(p1, l3);
# p1p = ppol(p0, p2, p1);
ph = (p0 + p2) / 2;
vn = vectorNorm(p2-p0);
x1 = 2 * vectorNorm(p1p-ph) / vn
y1 = 2 * vectorNorm(p1p-p1) / vn
w = 1 / sqrt(x1**2 + y1**2 + 1)

# "Characteristics of conic segments in Bezier form"
# Javier Sanchez-Reyes
# Proceedings of the IMProVe 2011, pp. 231-4
# (all calculations between b0 and F2 are complex)
b0 = p0(1) + p0(2) * i;
b1 = p1(1) + p1(2) * i;
b2 = p2(1) + p2(2) * i;

alpha = 1/(1-w**2)
m = (b0 + b2)/2;
C = (1-alpha) * b1 + alpha * m
d = (1-alpha) * b1**2 + alpha * b0 * b2;
c = sqrt(C**2 - d);
F1 = C + c
F2 = C - c
a = (abs(F1 - b0) + abs(F2 - b0))/2
b = sqrt(a**2 - abs(c)**2) # Note: This is reversed in paper
ang = rad2deg(arg(c))

figure;
hold on;
axis equal;
drawEllipse(real(C), imag(C), a, b, ang, 'color', 'magenta');
drawLine(l0, 'color', 'cyan');
drawLine(l2, 'color', 'cyan');
drawLine(l3, 'color', 'blue');
#drawPoint(ph, 'color', 'blue');
drawPoint(p0, 'color', 'black');
drawPoint(p1, 'color', 'blue');
drawPoint(p2, 'color', 'black');
waitforbuttonpress();

```

