Here is an expanded version of the generating function argument I sketched in a comment.
For $i=1,2,3$, define the generating functions $F_i(x,y) := \sum_{n=0}^\infty \sum_{q=0}^\infty R(n,3q+i) x^n y^q$, which are well defined for $x,y$ small. If one starts with the recursive identities
\begin{align*} R(n,3q) &= \sum_{0 \leq r \leq q} R(n-1,3r) + \sum_{0 \leq r \leq q} R(n-1,3r+1) \\
&\quad + \sum_{0 \leq r \leq q-1} R(n-1,3r+2)\\
R(n,3q+1) &= \sum_{0 \leq r \leq q+1} R(n-1,3r) + \sum_{0 \leq r \leq q} R(n-1,3r+1) \\
&\quad + \sum_{0 \leq r \leq q} R(n-1,3r+2)\\
R(n,3q+2) &= \sum_{0 \leq r \leq q+1} R(n-1,3r) + \sum_{0 \leq r \leq q+1} R(n-1,3r+1) \\
&\quad + \sum_{0 \leq r \leq q+1} R(n-1,3r+2)
\end{align*}
for $n \geq 1$, multiplies by $x^n y^q$, and then sums using the geometric series formula and the initial condition $R(0,3q+i)=1$, one obtains after some calculation the equations
\begin{align*}
F_0(x,y) &= \frac{1}{1-y} + \frac{x}{1-y}(F_0(x,y) + F_1(x,y) + y F_2(x,y))\\
F_1(x,y) &= \frac{1}{1-y} + \frac{x}{1-y}\left(\frac{F_0(x,y)}{y} + F_1(x,y) + F_2(x,y)\right) - \frac{x\alpha(x)}{y}\\
F_2(x,y) &= \frac{1}{1-y} + \frac{x}{1-y}\left(\frac{F_0(x,y)}{y} + \frac{F_1(x,y)}{y} + \frac{F_2(x,y)}{y}\right) - \frac{x\beta(x)}{y}
\end{align*}
for almost all small $x,y$, where
$$ \alpha(x) := \sum_{n=0}^\infty R(n,0) x^n$$
and
$$ \beta(x) := \sum_{n=0}^\infty (R(n,0)+R(n,1)+R(n,2)) x^n.$$
Note that it is $\alpha$ that we want to understand. So the strategy will be to eliminate the other unknowns $F_0,F_1,F_2,\beta$ to isolate a formula for $\alpha$.
We have a linear system of three equations in three unknowns $F_0,F_1,F_2$. Solving this system using a standard symbolic algebra package, one can eliminate these unknowns, obtaining for instance
$$ F_1(x,y) = \frac{(x^3-yx^2 - y^2 x + yx)\alpha(x) +yx^2 \beta(x) - y^2+x^2}{P(x,y)}$$
for almost all small $x,y$, where $P$ is the (irreducible) cubic
$$ P(x,y) := y^3 - (1-2x) y^2 + xy - x^3;$$
there are similar formulae for $F_0$ and $F_2$ that we shall discard (they give equivalent constraints to the one (1) we will end up using).
Since $F_1$ is analytic at the origin, we conclude the constraint
$$ (x^3-yx^2 - y^2 x + yx)\alpha(x) +yx^2 \beta(x) - y^2+x^2 = 0 \quad (1)$$
whenever $x,y$ are small and $P(x,y)=0$. So now the main challenge is to use this relation (1) to eliminate $\beta$.
When $x=0$, the equation $P(x,y)$ has a double zero at $y=0$. Thus for small $x$, there are two small solutions $y_1,y_2$ to $P(x,y)=0$ and one large solution $y_3$ (which is near $y=1$, since $P(0,1)=0$). Since (1) holds for $y=y_1$ and $y=y_2$, we may eliminate $\beta(x)$ to conclude after some algebra that
$$ \alpha(x) = -\frac{x^2 + y_1 y_2}{x (x^2 + y_1 y_2 - x)}.$$
However, as $y_1,y_2,y_3$ are the roots of $P(x,y)=0$ we have $y_1 y_2 y_3 = x^3$, so we can simplify to
$$ \alpha(x) = \frac{y_3+x}{y_3 - xy_3 - x^2}.$$
From the implicit function theorem $y_3$ is an analytic function of $x$ for $x$ small, so this in fact describes $\alpha$ completely as an element of ${\bf Q}(x,y_3) \equiv {\bf Q}(x,y)/(P(x,y))$, which is a cubic extension of ${\bf Q}(x)$ and should therefore obey a cubic equation with coefficients in ${\bf Q}(x)$. Indeed, using a symbolic algebra package, one can verify the identity
$$ x^3 \alpha(x)^3 + (4x^2-3x+1) \alpha(x)^2 + (5x-3) \alpha(x) + 2 = \frac{x P(x,y_3)}{(y_3 - xy_3 - x^2)^3};$$
but $P(x,y_3)$ vanishes, hence
$$ x^3 \alpha(x)^3 + (4x^2-3x+1) \alpha(x)^2 + (5x-3) \alpha(x) + 2 = 0.$$
Writing $A(x) := x\alpha(x)$, we then get the generating function identity
$$ x^2 A(x)^3 + (4x^2-3x+1) A(x)^2 + (5x^2-3x) A(x) + 2x^2 = 0$$
for the sequence $a(n)$ in A301897. This uniquely specifies $A$ if we enforce the asymptotics $A(x) = x + O(x^2)$, so we obtain $R(n,0)=a(n+1)$ as claimed. With a similar effort one could obtain explicit formulae for $\beta, F_0, F_1, F_2$ which would lead eventually to some combinatorial formula for $R(n,q)$; one could also analyze the singularities of these generating functions to obtain asymptotics for these sequences using standard analytic combinatorics methods (e.g., the residue theorem).
EDIT: I have some further commentary regarding the means I arrived at this answer here.
Best Answer
TL;DR: I have a proof of your conjectured asymptotic formula, modulo the correctness of a certain alternative description of your $c_n$ sequence.
I tried to complete Iosif Pinelis's elegant analysis by finding a way to derive the value of the constant $2/\pi$ (denoted $a$ in Iosif's answer) from the quadratic recurrence relation. The problem with that relation is that the behavior it implies for $c_n$ for large $n$ seems to depend in a very sensitive way on the initial values of the sequence, so I concluded that this approach has little chance of working.
Fortunately, I've now discovered another, linear recurrence relation for the same sequence $c_n$ that has better behavior and gives your claimed asymptotics without much effort. The relation is:
$$ c_0=1, \qquad c_n = g_n + \sum_{k=1}^n h_k c_{n-k} \quad (n\ge 1), \qquad (*) $$ where I define \begin{align} g_n &= \frac{((2n)!)^2}{2^{3n}(n!)^3}, \\ h_n &= \frac{((2n)!)^2}{2^{3n}(n!)^3}\cdot \frac{2n-1}{2n+1}. \end{align} (Edit: this corrects a small typo from the earlier version. As you pointed out in a comment, $g_n$ can also be expressed as $\frac{2^n \Gamma(n+1/2)^2}{\pi \Gamma(n)}$.)
I haven't verified rigorously that this relation is equivalent to your quadratic relation, but numerically it gives the correct sequence 1, 2, 6, 24, 126, 864, 7596, ..., and I believe this should be straightforward to prove. The reasoning that led me to it involves your description of the sequence as coming from an asymptotic expansion for a ratio of two Bessel functions. I started with the relation $$ K_0(-1/4z) = K_1(-1/4z) \times \sum_{n} c_n z^n, $$ and, expanding both Bessel functions in a power series (actually not quite a traditional power series because of some nasty-looking transcendental terms, but those can be factored out), massaged this into a linear system of equations satisfied by the $c_n$'s, which gave me the recurrence after a bit of additional guesswork.
To rigorously prove the relation, one can work with this Bessel function picture and do the analysis more carefully, or one can try to prove directly that the linear recurrence is equivalent to the quadratic recurrence without any reference to Bessel functions. I suspect this is doable through an inductive argument, probably involving formulating and proving some auxiliary hypergeometric summation identities.
Finally, if we assume that $(*)$ is correct, we can prove your claim that $$ c_n \sim \frac{2}{\pi} 2^n (n-1)!. $$ Observe that $$ \frac{c_n}{2^n (n-1)!} = \frac{g_n}{2^n (n-1)!} + \frac{h_n}{2^n (n-1)!} + \frac{1}{2^n (n-1)!} \sum_{k=1}^{n-1} h_k c_{n-k} $$ Using Stirling's formula, you can check that each of the first two terms in this expression converges to $1/\pi$. The third term (the normalized sum) can be easily shown to be $O(1/n)$ (the first summand is $O(1/n)$, and the remaining summands are $O(1/n^2)$ and there are $n-2$ of them). Here I am using the fact that the sequence $c_n/2^n (n-1)!$ is bounded, as shown in Iosif Pinelis's answer (and as can probably also be shown from the linear recurrence without much effort).