Characterizing BV functions via the Stieltjes integral
This is a straightforward application of the Uniform Boundedness Theorem.
For any partition $P = (x_0,x_1, \ldots,x_n)$ of $[a,b]$ define the linear operator $T_P:C([a,b]) \to \mathbb{R}$ by
$$T_Pf = \sum_{j=1}^nf(x_j)(g(x_j) - g(x_{j-1}).$$
Here $C([a,b])$ is taken to be the Banach space equipped with the supremum norm $\|f\|_\infty.$
Since every $f$ is integrable with respect to $g$, we have
$$\lim_{\|P\| \to 0} T_Pf = \int_a^b f \, dg.$$
and, hence, for every $f$ there exists $M_f > 0$ such that $\sup_P|T_P f| \leqslant M_f$.
By the Uniform Boundedness Theorem we have uniform boundedness in operator norm and
$$\sup_{f \in C([a,b]), \|f\|_\infty = 1} \sup_P|T_P f| \leqslant M < \infty.$$
For any partition $P$ we can construct a continuous (piecewise linear function) $f$ such that $f(x_j) = 1 $ if $g(x_j) - g(x_{j-1}) \geqslant 0$ and $f(x_j) = -1$ if $g(x_j) - g(x_{j-1}) < 0.$ Thus
$$\sum_{j=1}^n |g(x_j) - g(x_{j-1})| \leqslant M < \infty$$
and $g$ must be of bounded variation.
Here is a proof following what you had in mind. It's a bit lengthy, there is probably a much shorter proof using some powerful theorem of functional analysis, like Riesz representation theorem for signed measures for example.
Take some function $g:[0;1]\to \mathbb R$ to be not of bounded variation. This means that tere exists some $x\in I=[0;1]$ such that for every $\varepsilon>0$ the restriction of $g$ to $[x-\varepsilon;x+\varepsilon]$ is not of bounded variation. To prove this we assume that the oposite is true, for every $x$ there is some $\varepsilon>0$ such that the restriction of $g$ to $]x-\varepsilon;x+\varepsilon[$ is of bounded variation. This gives us a covering of $[0;1]$ by open sets, by compacity we can extract a finite sub-couvering of $I$ by open intervals $I_1,\ldots , I_n$. The total variation of $g$ over $I$ is less than the sum of the variations of $g$ over the intervals $I_i$, which is finite. This contradict the fact that $g$ is not of bounded variation and proves the existence of such an $x$. Moreover we must have that every restriction of $g$ to $[x;x+\varepsilon]$ is not of bounded variation or that every restriction of $g$ to $[x-\varepsilon;x]$ is not of bounded variation (we can eventually have both). Otherwise, the restriction of $g$ to some $[x-\varepsilon;x+\varepsilon] $ would be of bounded variations.
Without loss of generality (and to simplify notations) we can assume that $x=0$ and that every restriction of $g$ to some $[0,\varepsilon]$ is not of bounded variation. Now take a finite number of points $0<t^1_1<\ldots <t^1_{n(1)}\leq 1$ such that $t_1^1<\frac 1 {2^1}$ and $$\sum_{i=1}^{n(1)-1} |g(t^1_{i+1})-g(t^1_i)|\geq 1.$$ Notice that we have $t^1_1>0$, such a collection of points must exists because $g$ is not of bounded variation over $[0;1]$. After that we take a finite number of points $0<t^2_1<\ldots< t^2_{n(2)}<t^1_1$ such that $t^2_1<\frac{1}{2^2}$ and $$\sum_{i=1}^{n(2)-1} |g(t^2_{i+1})-g(t^2_i)|\geq 1.$$ Again this set of point exists because $g$ is not of bounded variation over $[0;t_1^1]$. Repeting this process ad vitam eternam we get an infinite sequence of numbers $$1\geq \tau_1>\tau_2>\ldots >\tau_n>\ldots=1\geq t^1_{n(1)}>\ldots>t^1_1>t^2_{n(2)}>\ldots>t^2_1>\ldots t^n_m>\ldots$$ (We just ordered the $t^i_j$ and renamed them as $\tau_k$) This sequence $(\tau_i)_{i\in \mathbb N}$ is strictly decreasing, has limit $0$ (because $t^i_i<2^{-i}$) and we have $$\sum_{i=1}^\infty |g(\tau_i)-g(\tau_{i+1})|=\infty. $$
We are now ready to construct the $f$ such that $\int f \mathrm d g $ does not exists. First we use a lemma telling us that if $(a_n)_n$ is a sequence of positive real numbers such that $\sum a_n=+\infty$ then there exists some sequence $(\varepsilon_n)_n$ of positive real numbers with limit $0$ such that $\sum a_n \varepsilon_n= +\infty$. Applying this lemma to the sequence $(|g(\tau_{n+1})-g(\tau_{n})|)_n$ we get a sequence $(\varepsilon_n)_n$ with limit $0$ such that $$\sum_{n=1}^\infty \varepsilon_n\mathrm{sgn}\big(g(\tau_{n+1})-g(\tau_{n})\big)\big(g(\tau_{n+1})-g(\tau_{n})\big)=+\infty$$ where $\mathrm{sgn}(x)$ is $1$ if $x>0$, $-1$ if $x<0$ and $0$ otherwise. We define $f$ by $f(0)=0$, $f(\tau_n)=\varepsilon_n\mathrm{sgn}\big(g(\tau_{n+1})-g(\tau_{n})\big)$ for every $n$ and we connect these points by affine functions. By construction we have $$ \sum_{n=1}^\infty f(\tau_n)\big(g(\tau_{n+1})-g(\tau_{n})\big)=+\infty $$ so $\int f\mathrm d g$ doesn't exists. We just have to verify that $f$ is indeed continuous : $f$ is continuous on every $[\varepsilon,1]$ because it is piecewise linear on these intervals. $f$ is also continuous in $0$ because $\lim \varepsilon_n =0$ and $f(0)=0$. This in turn imply that $f$ is continuous over $[0;1]$ and finishes the proof.