Minkowski Metric Signature

as you wrote, the spacetime invariant can be expressed as: $$ds^2=g_{\mu\nu}dx^{\mu}dx^{\nu}$$ and from that we normally get: $$ds^2=-c^2dt^2+dx^2+dy^2+dz^2$$ This is not because of some arbitrary imaginary time unit, this is because the metric ($g_{\mu\nu}$) is a diagonal matrix with the coefficients of each term of the $ds^2$ equation: $$g_{\mu\nu}=\left(\begin{array}{l}-1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{array}\right)$$ and the coordinates are listed as you would assume: $$dx^{\mu}=\left(\begin{array}{l}cdt\\dx\\dy\\dz\end{array}\right)$$ Then, you should note that $$g_{\mu\nu}dx^{\mu}=dx_{\nu}=\left(\begin{array}{l}-1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{array}\right)\left(\begin{array}{l}cdt\\dx\\dy\\dz\end{array}\right)=\left(-cdt~~dx~~dy~~dz\right)$$ Also, $v_{\mu}v^{\mu}$ is the inner product, meaning: $$dx_{\nu}dx^{\nu}=\left(-cdt~~dx~~dy~~dz\right)\left(\begin{array}{l}cdt\\dx\\dy\\dz\end{array}\right)=-c^2dt^2+dx^2+dy^2+dz^2$$ This is the equation you want without any imaginary unit omission. The reason for the $-1$ in the $g_{\mu\nu}$ is that it makes the system Lorentz invariant; it maintains $ds^2$ as a spacetime invariant quantity.

Let me be historical. In Euclidean 3-D coordinates, you find the interval between positions as $$\Delta d_{Eucl}^2=(X_2-X_1)^2+(Y_2-Y_1)^2+(Z_2-Z_1)^2$$ When incorporating relativity and time, the interval becomes a spacetime quantity. Because relativity sets the maximum speed of information as $c$, we make the interval $$\Delta s^2=\Delta d_{Eucl}^2-c^2(t_2-t_1)^2$$ This represents the original interval - the distance between the two events - minus the maximum distance the information could travel in the time between the two events. That difference lets us determine if the events happened in a definite chronological order ($\Delta s^2<0$) or if they occurred in two distinctly separate positions ($\Delta s^2>0$), since in relativity we can't always be sure. It is from this that the $-1$ in the metric originates. Space and time coordinates are given opposite signs here. We keep the metric in terms of $s^2$ because we simply can't be sure if $s$ is positive or negative. There was no original imaginary time coordinate, that was simply someone's poor interpretation and it has been (thankfully) dropped for the most part.

I should probably also point out that the imaginary time coordinate can not come out of Euclidean 4-D either. If one ignores relativity, then there is no maximum velocity. If there is no maximum velocity, there is no natural way of equating spatial and temporal coordinates. Therefore, not only would it not be right to use $c$ in the $ict$ coordinate, it also would not make sense to add time to space because there would be no agreeable conversion between them. However, if you don't ignore relativity, then you must subtract the time term from the 3-D interval in order to comply with the notion of a maximum velocity. So the Euclidean signature, $(1,1,1,1)$ can not be used to describe 4-D spacetime! So you never define the time coordinate as imaginary.

If you define $x^0=ict$, then I assume one takes $x_\mu=x^\mu$ so that the metric is actually $\eta_{\mu\nu}=\text{diag}(1,1,1,1)=\delta_{\mu\nu}$, i.e. you're dealing with a Euclidean metric. Then $$ds^2=\delta_{\mu\nu}dx^\mu dx^\nu$$ gives the usual outcome : $$ds^2=-c^2dt^2+d\vec{x}^2$$ The usual conventions are as follows:

Option one: One defines $x^\mu=(ct,\vec{x})$ and $x_\mu=\eta_{\mu\nu}x^\nu=(-ct,\vec{x})$ where $\eta_{\mu\nu}=\text{diag}(-1,1,1,1)$. This results in $$ds^2=-c^2dt^2+\vec{x}^2$$This convention is usually taken in treatments that focus on (general) relativity and/or spacetime structure.

Option two: One defines $x^\mu=(ct,\vec{x})$ and $x_\mu=\eta_{\mu\nu}x^\nu=(ct,-\vec{x})$. where $\eta_{\mu\nu}=\text{diag}(1,-1,-1,-1)$. This results in $$d\tilde{s}^2=-ds^2=c^2dt^2-\vec{x}^2$$ This approach is usually taken when the focus is on particle physics results. My personal theory is that this is because it results in the equation $p_\mu p^\mu=m^2$ (as opposed to $p_\mu p^\mu=-m^2$), which is one of the main equations in particle physics and arguably looks slightly more pleasing than the alternative. As was pointed out in a comment, it also makes time-like intervals positive, which might be preferred when dealing with particles (which, of course, always travel along time-like or light-like trajectories).

Of course, the two approaches are completely equivalent to each other. The convention with the imaginary time coordinate has fallen out of use a little bit. I can see why this would happen: It's not a particularly helpful convention for intuition, nor does it seem to carry over well to general relativity.

The only real reason to introduce $ict$ coordinates is to stress the similarity (for didactic purposes I guess) between Lorentz transformation and orthogonal rotations in more used-to Euclidian space.

Note that Minkowski pseudo-Euclidian space obtains exactly "normal" Euclidian form if complex time is introduced,namely: metric signature becomes $++++$: exactly like if it was regular Euclidian 4-space. Also, more vivid: matrix of Lorentz transformation obtains exactly form of real orthogonal matrix due to $\cos(ix) = \cosh x$ (similar for $\sin x$). Thus, you rotate through complex angle, but matrix looks like regular orthogonal, for example boost in $x$-direction, $$ \left( \begin{matrix} \cos z & \sin z & 0 & 0\\ -\sin z & \cos z & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ \end{matrix} \right) $$

where $z$ is now strictly imaginary.