Perhaps, the form given by Lagrange inversion theorem cannot be much simplified here. It expresses the $n$-th coefficient of a series reversion as the sum of $n-1$ values of exponential Bell polynomials.
From the practical perspective, since $\theta_3-1$ contains nonzero coefficients only at square powers, computation of the $n$-th reversion coefficient amounts to iterating over the partitions of $n-1+k$ into $k$ nonzero squares (OEIS A243148) for $k\in\{1,2,\dots,n-1\}$.