Here's a reasonable method: translate everything such that the center of the ellipse is at the origin. Consider the intersection of the ellipse with major axis $2a$ and minor axis $2b$ with the polar equation
$$r=\frac{a b}{\sqrt{b^2\cos^2\theta+a^2\sin^2\theta}}$$
and the line $\tan\,\theta=\dfrac{y_2-y_1}{x_2-x_1}$. (When solving the last equation for $\theta$, you will want to use the two-argument arctangent that is implemented in most computing environments.) Once having computed the corresponding values of $r$ at $\theta$ and $\pi+\theta$, convert to rectangular coordinates and translate back to your initialal origin.
Matrix derivation
The original equation can be described in a more concise form as
$$\vec x A \vec x^T + \vec x \vec b^T + \vec b \vec x^T + f =0$$
or
$$\vec x A \vec x^T + 2\vec x \vec b^T + f =0$$
where $\vec b=(d,e)$ and $A=\begin{pmatrix}a&b\\b&c\end{pmatrix}$.
If the vector $\vec x_0$ represents the center then in a coordinate system with the origin in this point the equation will be
$$(\vec x-\vec x_0) A (\vec x-\vec x_0)^T + g =0$$
which can be rewritten as
$$\vec x A \vec x^T - 2\vec x A\vec x_0^T + \vec x_0 A \vec x_0^T +g.$$
Which means the we need to choose $\vec x_0$ in such way that $-2A\vec x_0^T=2\vec b^T$, i.e. $$A\vec x_0^T=-\vec b^T,$$
$$\begin{pmatrix}a&b\\b&c\end{pmatrix}
\begin{pmatrix}x\\y\end{pmatrix}=
\begin{pmatrix}d\\e\end{pmatrix}.$$
This is precisely the following system of linear equations written in the matrix form:
$$
\begin{align*}
ax + by &= -d\\
bx + cy &= -e
\end{align*}
$$
(To give credit where credit is due, I was shown this derivation by a colleague of mine who taught labs for the same subject as I did, but with another group of students.)
Derivatives
We can view the curve as a level set of the function
$$F(x,y)=ax^2+2bxy+cy^2+2dx+2ey.$$
The level sets of this function create the family of curves $ax^2+2bxy+cy^2+2dx+2ey+f=0$, where the parameter $f$ is changing.
For example, if $\delta>0$ we get concentric ellipses, like in this example - the plots are taken from WolframAlpha:
And a different plot:
If $\delta>0$, we get a family of hyperbolas. WolframAlpha:
And a different plot:
In the first case, the function $F$ has minimum at the center of each of the ellipse. In the second case there is a saddle point. In both cases, the center fulfills
$$\frac{\partial F}{\partial x}= \frac{\partial F}{\partial y} =0$$
which leads to the system of equations
\begin{align*}
ax + by +d&= 0\\
bx + cy +e&= 0
\end{align*}
which is precisely the linear system described above.
So far we have only drawn some pictures (which is good for intuition, but picture is not a proof.) If we want to give a more rigorous reason why this works we can use some known facts about quadratic forms.
We know from Principal axis theorem that there is a rotation and translation such that in the new coordinates the curve has
$$\lambda_1 x^2 + \lambda_2 x^2=const$$
Clearly $F(x,y)=\lambda_1 x^2 + \lambda_2 x^2$ has partial derivatives equal to zero at the origin. (And depending on the signs of $\lambda_{1,2}$ we can say whether it has minimum, maximum or saddle point there.)
Best Answer
Assume WLOG that $A>0$ (otherwise multiply the equation by $-1$). We have $$ \left[\matrix{x\\ y}\right]^T\left[\matrix{A & B\\ B & C}\right]\left[\matrix{x\\ y}\right]+2\left[\matrix{D\\ E}\right]^T\left[\matrix{x\\ y}\right]+F=0. $$ After completing the square we get $$ \left(\left[\matrix{x\\ y}\right]+\left[\matrix{A & B\\ B & C}\right]^{-1}\left[\matrix{D\\ E}\right]\right)^T\left[\matrix{A & B\\ B & C}\right]\underbrace{\left(\left[\matrix{x\\ y}\right]+\left[\matrix{A & B\\ B & C}\right]^{-1}\left[\matrix{D\\ E}\right]\right)}_{=z}=\\ =\underbrace{\left[\matrix{D\\ E}\right]^T\left[\matrix{A & B\\ B & C}\right]^{-1}\left[\matrix{D\\ E}\right]-F}_{=\alpha^2}. $$ If we factorize the matrix of the quadratic form (which is positive definite) $$ M=\left[\matrix{A & B\\ B & C}\right]^{1/2} $$ the ellipse equation becomes $$ \left\|\frac{1}{\alpha}Mz\right\|=1 $$ which means that the image of the ellipse under the linear mapping $y=\frac{1}{\alpha}Mz$ is the unit circle with area $\pi$. Area deformation under linear/affine maps is the determinant, i.e. the area of the ellipse $S$ satisfies $$ \pi=\det(\frac{1}{\alpha}M)S=\frac{1}{\alpha^2}\det(M)S\qquad\Rightarrow\qquad S=\frac{\pi\alpha^2}{\det(M)}=\frac{\pi\alpha^2}{\sqrt{AC-B^2}}. $$ Finally, to calculate $\alpha^2$ we use the Schur determinant formula $$ \Delta=\det\left[\matrix{A & B & D\\ B & C & E\\D & E & F}\right]= \det\left[\matrix{A & B\\ B & C}\right]\cdot (-\alpha^2)=-(AC-B^2)\alpha^2 $$ which gives $$ S=\frac{-\pi\Delta}{(AC-B^2)^{3/2}}. $$