I came across this question while making some calculations.
QUESTION. Can you find some transformation to "decouple" the double integral as follows? $$\int_0^{\frac{\pi}2}\int_0^{\frac{\pi}2}\frac{d\alpha\,d\beta}{\sqrt{1-\sin^2\alpha\sin^2\beta}} =\frac14\int_0^{\frac{\pi}2}\frac{d\theta}{\sqrt{\cos\theta\,\sin\theta}}\int_0^{\frac{\pi}2}\frac{d\omega}{\sqrt{\cos\omega\,\sin\omega}}.$$
(Thanks go to Etanche and Jandri) \begin{align}J&=\int_0^{\frac{\pi}{2}}\int_0^{\frac{\pi}{2}} \frac{1}{\sqrt{1-\sin^2(\theta)\sin^2 \varphi}}d\varphi d\theta\\ &\overset{z\left(\varphi\right)=\arcsin\left(\sin(\theta)\sin \varphi\right)}=\int_0^{\frac{\pi}{2}} \left(\int_0^ \theta\frac{1}{\sqrt{\sin(\theta-z)\sin(\theta+ z)}}dz\right)d\theta\tag1\\ &=\frac{1}{2}\int_0^{\frac{\pi}{2}} \left(\int_{u}^{\pi-u}\frac{1}{\sqrt{\sin u\sin v}}dv\right)du \tag2\\ &=\frac{1}{2}\int_0^{\frac{\pi}{2}} \left(\int_{u}^{\frac{\pi}{2}}\frac{1}{\sqrt{\sin u\sin v}}dv\right)du+\underbrace{\frac{1}{2}\int_0^{\frac{\pi}{2}} \left(\int_{\frac{\pi}{2}}^{\pi-u}\frac{1}{\sqrt{\sin u\sin v}}dv\right)du}_{w=\pi-v}\\ &=\int_0^{\frac{\pi}{2}} \left(\int_{u}^{\frac{\pi}{2}}\frac{1}{\sqrt{\sin u\sin v}}dv\right)du\\ &\overset{u\longleftrightarrow v}=\int_0^{\frac{\pi}{2}} \left(\int_{0}^{u}\frac{1}{\sqrt{\sin u\sin v}}dv\right)du\\ &=\boxed{\frac{1}{2}\int_0^{\frac{\pi}{2}} \int_0^{\frac{\pi}{2}}\frac{1}{\sqrt{\sin u\sin v}}dudv} \end{align}
and to obtain the form in the OP, finally substitute $u=2\theta $, $v=2\omega $.
$(1)$: $\displaystyle dz=\dfrac{\sqrt{\sin^2\theta-\sin^2 z}}{\sqrt{1-\sin^2 z}}d\varphi$, $\sin^2 a-\sin^2 b=\sin(a-b)\sin(a+b)$
$(2)$: Change of variable $u=\theta-z,v=\theta+z$
It might help to put the question in a bit broader perspective. Elliptic integrals provide a variety of such "decoupling" expressions, as listed in Definite Integrals of the Complete Elliptic Integral.
For example, the following identity holds for any $z>$1, $$ \int_0^{\frac{\pi}2}\int_0^{\frac{\pi}2}\frac{d\alpha\,d\beta}{\sqrt{1-\sin^2\alpha\sin^2\beta}}\frac{\sin\alpha}{\sqrt{z^2-\sin^2\alpha}} =\left(\int_0^{\frac{\pi}2}\frac{d\theta}{\sqrt{z+\cos 2\theta}}\right)^2,\qquad(1)$$ and there is a similar expression for $0<z<1$.
There are also closed form expressions for double integrals that include the one in the OP (for $s=1$), $$\int_0^{\frac{\pi}{2}}\int_0^{\frac{\pi}{2}}\frac{\sin^{s-1}2\alpha\;d\alpha\,d\beta}{\sqrt{1-\sin^2\alpha\sin^2\beta}}=\left(\frac{\pi\Gamma(s/2)}{2\Gamma(3/4)\Gamma(s/2+1/4)}\right)^2,\;\;{\rm Re}\,s>0.\qquad (2)$$
For $s=1$ the right-hand-side equals $K(1/2)^2$, and it might be interesting to see if there is a more general expression in terms of the square of an elliptic integral (perhaps by evaluating $K$ at one of the elliptic singular values, when the complete elliptic integrals may be computed in analytic form in terms of gamma functions).
Following Iosif Pinelis's suggestion, the last integral can be worked out by substituting the series expansion of the elliptic integral (and transforming to $u=\sin\alpha$) $$\int_0^{\frac{\pi}{2}}\int_0^{\frac{\pi}{2}}\frac{\sin^{s-1}2\alpha\;d\alpha\,d\beta}{\sqrt{1-\sin^2\alpha\sin^2\beta}}=2^{s-1}\frac{\pi}{2}\sum_{n=0}^\infty\left(\frac{(2n)!}{2^{2n}(n!)^2}\right)^2\int_{0}^1 u^{2 n} u^{s-1} \left(1-u^2\right)^{\frac{s}{2}-1}\,du$$ $$\qquad=2^{s-3}\Gamma \left(\frac{s}{2}\right)\sum_{n=0}^\infty\frac{\Gamma \left(n+\frac{1}{2}\right)^2 \Gamma \left(n+\frac{s}{2}\right)}{\Gamma (n+1)^2 \Gamma (n+s)}.\qquad (3)$$ Mathematica evaluates the sum over $n$ of the gamma functions to give the expression in formula 2. I don't know how Mathematica does it, and have asked here --- with several instructive answers!
$\newcommand\al\alpha\newcommand\be\beta\newcommand\th\theta\newcommand\om\omega$Making substitutions $\sin^2\al=x$, $\sin^2\be=y$, $\sin2\th=s$, and $\sin2\om=t$, we rewrite the equality in question in the following "algebraic" way: $$\int_0^1\int_0^1\frac{dx\,dy}{\sqrt{xy(1-x)(1-y)(1-xy)}} =\int_0^1\int_0^1\frac{2ds\,dt}{\sqrt{st(1-s)(1-t)(1+s)(1+t)}}.$$ Perhaps this could help, especially because the integrands on the left and on the right now look "more similar" to each other than in the original setting.
A further comment: The main difficulty here appears to be with the integral, say $L$, on the left. Expanding $\dfrac1{\sqrt{1-xy}}$ in powers of $xy$:
$$\dfrac1{\sqrt{1-xy}}=\sum_{m=0}^\infty a_mx^my^m$$
with $a_m:=(-1)^m\binom{-1/2}m$ and using the expression of the beta function in terms of the gamma function, we have
$$L=\sum_{m=0}^\infty a_m\Big(\int_0^1 dx\, x^{m-1/2}(1-x)^{-1/2}\Big)^2
=\frac{\pi^2}4\,N,$$
where
$$N:=\sum_{m=0}^\infty\frac1{4^{3m}}\,\binom{2m}m^3.$$
According to Mathematica,
$$N=\frac\pi{\Gamma(3/4)^4}=1.393203929\ldots\quad\text{and hence}\quad L=\frac{\pi^3}{4\Gamma(3/4)^4};$$
however, I do not know how Mathematica does this.
One may also note that $N$ is the expected number of times a symmetric simple random walk on $\mathbb Z^3$ visits its starting point. That is, $$N=\sum_{m=0}^\infty P(X_{2m}=Y_{2m}=Z_{2m}=0) =E\sum_{m=0}^\infty 1(X_{2m}=Y_{2m}=Z_{2m}=0),$$ where $X_n:=R_1+\cdots+R_n$, the $R_i$'s are iid Rademacher random variables (with $P(R_i=\pm1)=1/2$), and $(Y_n)$ and $(Z_n)$ are iid copies of $(X_n)$. It seems amusing that, whereas the dimension $3$ of $\mathbb Z^3$ is odd, $N$ is the square of a certain nice expression.
Closely related to $N$ is Pólya's random walk constant, equal $1-1/N$, which is the probability that the random walk will ever return to its starting point.
The equality $$N=\frac\pi{\Gamma(3/4)^4}=1.393203929\ldots$$ seems to be proved by Montroll, formulas (6.11), (6.12), (6.1).
