There is a way to calculate the location of the horizontal line that goes on top of the square.
There are actually two ways to do it, but I'm not an expert on matrices so I'm not going to go into detail on how 3D rendering/calculating.
The old way to do it is based off of the idea of having a pane of glass between your eye and the object and drawing straight lines from your eye to the object to find the intersect points on the pane of glass. This pane of glass (canvas plane) has the points correctly placed.
This post contains detailed information on this perspective drawing technique and proofs for it.
http://perspectiveresources.blogspot.com/p/principles.html
No numbers are actually required, just angles and straight lines. More euclidian geometry than anything.
Note: This technique works for all shapes and rotations of shapes from a plane into perspective. To make shapes have height, that takes more work.
I made an interactive graph on Desmos that allows you to drag points $B$ and $D$ and watch how points $A$ and $C$ change.
While it might not be the quickest way to get there, you can certainly get the solution in the way that you were attempting - with your function $g$ and the distance function, as you were on track to do!
To solve the problem this way, just set up the two equations, and solve. You want the two points $X \in \{A, C\} = \{(x_A, y_A), (x_C, y_C)\}$ that satisfy two conditions. First, they lie along the line $g$, i.e.
\begin{align}
y &= \frac{a}{b}x + \frac{b^2 - a^2}{2b}
\end{align}
Second, the distance from $B$ to $D$ is twice the distance from each $X$ to $M$, i.e.
\begin{align}
\sqrt{a^2 + b^2} &= 2 \sqrt{\left(x - \frac{a}{2}\right)^2 + \left(y - \frac{b}{2}\right)^2} \\ \Rightarrow
a^2 + b^2 &= 4 \left(x - \frac{a}{2}\right)^2 + 4 \left(y - \frac{b}{2}\right)^2 \\ \Rightarrow
a^2 + b^2 &= \left(2x - a\right)^2 + \left(2y - b\right)^2
\end{align}
Substitute the value for $y$ from the first equation into the second, and solve for $x$. You need two points, but the equation is quadratic, so that's a hint that you're probably gonna be OK and end up with 2 possible solutions which are your two points.
\begin{align}
a^2 + b^2 &= \left(2x - a\right)^2 + \left(2\left( \frac{a}{b}x + \frac{b^2 - a^2}{2b} \right) - b\right)^2 \\
a^2 + b^2 &= \left(2x - a\right)^2 + \frac{a^2}{b^2} \left(2x - a\right)^2 \\
a^2 + b^2 &= \frac{a^2 + b^2}{b^2} \left(2x - a\right)^2 \\
b^2 &= \left(2x - a\right)^2 \Rightarrow
\pm b = 2x - a
\end{align}
$$
\Rightarrow x = \frac{a \pm b}{2}
$$
The two solutions ($\pm$) give you the two values for $x$:
\begin{align}
x_A &= \frac{a + b}{2}\;,\; x_C = \frac{a - b}{2}
\end{align}
Substitute each solution for $x$ into the equation for $y$, to find the two associated values of $y$.
\begin{align}
y_A &= \frac{a}{b}\left(\frac{a+b}{2}\right) + \frac{b^2-a^2}{2b} = \frac{a + b}{2} \\
y_C &= \frac{a}{b}\left(\frac{a-b}{2}\right) + \frac{b^2-a^2}{2b} = \frac{-a + b}{2}
\end{align}
Therefore points $A$ and $C$ are given by
\begin{align}
A &= (x_A, y_A) = \left(\frac{a + b}{2} , \frac{a + b}{2}\right) \\
C &= (x_C, y_C) = \left(\frac{a - b}{2} , \frac{-a + b}{2}\right)
\end{align}
Looking at the interactive graph, the way we found $A$ and $C$ was to solve for the intersections between the circle centered at $M$ - which is the set of points at the same distance from $M$ as $B$ and $D$ - and the red line $g$ which is perpendicular to $\overline{BD}$ and which passes through $M$.
Best Answer
So we have two points: $A(x_1,y_1)$ and $B(x_2,y_2)$. The first step is to find midpoint $M(x_c,y_c)$. This is straightforward: $x_c=\frac{x_1+x_2}{2}, y_c=\frac{y_1+y_2}{2}$
The next step is to move the origin to point $M$. Our points $A$ and $B$ will have new coordinates: $A'(x_1-x_c,y_1-y_c)$ and $B'(x_2-x_c,y_2-y_c)$.
Finally, we rotate $A'B'$ about the origin $90$ degrees to get the other two vertices. The rotated points will be $A''(y_c-y_1, x_1-x_c)$ and $B''(y_c-y_2, x_2-x_c)$. Moving the origin back to the original location will give us coordinates of the vertices: $C(y_c-y_1+x_c, x_1-x_c+y_c)$ and $D(y_c-y_2+x_c, x_2-x_c+y_c)$. We don't really need to calculate half-diagonal but it's easy to see that $y_c-y_1=-y_d$ and so forth.