From here, the inverse of a $2\times2$ quaternion matrix is given by
$$ \begin{bmatrix} a & b \\ c & d \end{bmatrix}^{-1} =\frac{1}{\mathrm{Nrd}} \begin{bmatrix} |d|^2\overline{a}-\overline{c}d\overline{b} & |b|^2\overline{c}-\overline{a}b\overline{d} \\ |c|^2\overline{b}-\overline{d}c\overline{a} & |a|^2\overline{d}-\overline{b}a\overline{c} \end{bmatrix} \tag{1}$$
where the so-called "reduced norm" of the matrix is given by
$$ \mathrm{Nrd}\Big(\begin{bmatrix} a & b \\ c & d \end{bmatrix}\Big)=|a|^2|d|^2+|b|^2|d|^2-2\mathrm{Re}(\overline{b}a\overline{c}d). \tag{2}$$
This answer surprised me: I, proceeding like the OP, found two formulas involving noncommutative rational functions (so to speak) depending on whether the diagonal or antidiagonal entries were both nonzero. I was not expecting this.
But the user that posted this is long gone and gave no explanation. Via searching, I found the reduced norm $\mathrm{Nrd}(a)$ of an element $a\in A$ of a central simple algebra $A$ of dimension $n^2$ over a field $K$ is characterized by $\det(L_a)=(\mathrm{Nrd}\,a)^n$, where $L_a(x)=ax$ is considered a $K$-linear endomorphism of $A$ as a $K$-vector space. Or, equivalently, there is some field extension $L/K$ which is a "splitting field" for $A$ in the sense that $A\otimes L\cong M_n(L)$ and then $\mathrm{Nrd}$ is just $\det$ taken in $M_n(L)$ (but the result necessarily lands in $K$). The "reduced trace" $\mathrm{Trd}\,a$ can be defined similarly. User "johgi" used the notation $\mathrm{tr}$ for the reduced trace, and I figured out this is twice the real part.
Q1. Where does the formula $(2)$ for the reduced norm $\mathrm{Nrd}\,\big[\begin{smallmatrix}a&b\\c&d\end{smallmatrix}\big]$ come from?
Of course, a tedious calculation can verify the formula, but I expect cannot demystify its origin. I notice for reals it factors as $(ad-bc)^2$, so $(1)$ looks like a "rehomogenized" version of the adjugate formula for a matrix inverse.
Q2. Where does the formula $(1)$ for the $2\times2$ matrix inverse come from?
I can check that it is indeed an inverse, using $2\mathrm{Re}(x)=x+\overline{x}$ and $\mathrm{Re}(xy)=\mathrm{Re}(yx)=\mathrm{Re}(\bar{x}\bar{y})$ and $\mathrm{Re}(xyz)=\mathrm{Re}(zxy)=\mathrm{Re}(yzx)$. Multiplying the matrix with its inverse, it seems like the reduced norm has "pieces" which come together, so I wouldn't be surprised if the answers to Q1 and Q2 were intertwined.
I notice if we define $L_a(x)=ax$ and $R_a(x)=xa$ along with $X(p)=L_dR_a-L_cR_b$ and its adjoint $X^\ast(p)=L_{\overline{d}}R_{\overline{a}}-L_{\overline{c}}R_{\overline{b\,}}$ (viewing $\mathbb{H}$ as an inner product space), we have
$$ \begin{bmatrix} a & b \\ c & d \end{bmatrix}^{-1} = \frac{1}{\mathrm{Nrd}} \begin{bmatrix} \phantom{-}X^\ast(d) & -\overline{X(\overline{b})} \\ -X^\ast(c) & \phantom{-}\overline{X(\overline{c})} \end{bmatrix}, $$
which looks a lot like an adjugate. Dunno what to make of that.
Q3. Does this generalize to other dimensions and algebras?
$ \def\mc#1{\left[\begin{array}{r|rrr}#1\end{array}\right]} \def\m#1{\left[\begin{array}{r}#1\end{array}\right]} \def\qif{\quad\iff\quad} \def\M{{\cal M}} $Transform each component quaternion into its $4\times 4$ matrix representation, e.g. $$\eqalign{ a = \mc{a_1\\a_2\\a_3\\a_4} \qif &A = \mc{ a_1 & -a_2 & -a_3 & -a_4 \\\hline a_2 & a_1 & -a_4 & a_3 \\ a_3 & a_4 & a_1 & -a_2 \\ a_4 & -a_3 & a_2 & a_1 }, \quad A^{-1} = \frac{A^T}{|a|^2} \\ }$$ Construct the equivalent block matrix and calculate its inverse $$\eqalign{ &M = &\m{a&b\\c&d} &\qif \M = \mc{A&B\\\hline C&D} \\ &&\qquad M^{-1} &\qif \M^{-1} \\ }$$ via the Schur complements of $\M$. This will recover the formula of interest.
For quaternionic matrices larger than $2\times 2$, the block matrix approach remains valid for numerical purposes, but the corresponding symbolic formulas quickly become intractable.