The function and the region are symmetric about the line $y=x$, so you can re-write your integral as
$$8\int_0^1 \int_0^x f(x,y) \; dy \; dx.$$
So you're integrating over the triangle with vertices $(0,0),$ $(1,0),$ and $(1,1).$
The vertical line $x=1$ has polar equation $r\cos \theta =1$ or $r=\sec \theta.$ So in polar, the integral is
$$8\int_0^{\pi/4} \int_0^{\sec \theta} \sqrt{1+r^2}\; r \; dr \; d\theta.$$
Then things are straightforward.
Define:
\begin{eqnarray}
I(a,b,c,d):=\int_0^{\tan ^{-1}\left(\sqrt{\frac{b}{a}}\right)} \frac{\sin (\theta ) \tan ^{-1}\left(\sin (\theta ) \sqrt{\frac{d}{b+c \sin ^2(\theta )}}\right)}{\sqrt{a b d \left(b+c \sin ^2(\theta )\right)}} \, d\theta \quad (i)
\end{eqnarray}
Then also define:
\begin{eqnarray}
{\mathfrak F}^{(A,B)}_{a,b} &:=& \int\limits_A^B \frac{\log(z+a)}{z+b} dz\\
&=& F[B,a,b] - F[A,a,b] + 1_{t^* \in (0,1)} \left(
-F[A+(t^*+\epsilon)(B-A),a,b] + F[A+(t^*-\epsilon)(B-A),a,b] \right) \quad (ii)
\end{eqnarray}
where
\begin{eqnarray}
t^*:=-\frac{Im[(A+b)(b^*-a^*)]}{Im[(B-A)(b^*-a^*)]}
\end{eqnarray}
and
\begin{equation}
F[z,a,b] := \log(z+a) \log\left( \frac{z+b}{b-a}\right) + Li_2\left( \frac{z+a}{a-b}\right)
\end{equation}
for $a$,$b$,$A$,$B$ being complex.
Then we have:
\begin{eqnarray}
I(a,b,c,d)&=& \frac{1}{\sqrt{a}} \int\limits_0^{\sqrt{\frac{d}{a+b+c}}} \frac{u \tan ^{-1}(u)}{\left(d-c u^2\right) \sqrt{d-u^2 (b+c)}} du\\
&=& \frac{\sqrt{d}}{\sqrt{a} (b+c)}\int\limits_0^{\sin ^{-1}\left(\sqrt{\frac{b+c}{a+b+c}}\right)} \frac{\sin (\phi ) \tan ^{-1}\left(\sin (\phi ) \sqrt{\frac{d}{b+c}}\right)}{d-\frac{c d \sin ^2(\phi )}{b+c}} d\phi\\
&=&-\frac{2 i}{\sqrt{a} \sqrt{d}} \int\limits_0^{\frac{\sqrt{\frac{b+c}{a+b+c}}}{\sqrt{\frac{a}{a+b+c}}+1}} \frac{t}{b \left(t^2+1\right)^2+c \left(t^2-1\right)^2} \log \left(\frac{2 i t \sqrt{\frac{d}{b+c}}+t^2+1}{-2 i t \sqrt{\frac{d}{b+c}}+t^2+1}\right) dt\\
&=&\frac{1}{4} \frac{1}{\sqrt{a b c d}} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^4
(-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\left\lfloor \frac{\xi -1}{2}\right\rfloor }
\int\limits_0^{\frac{\sqrt{\frac{b+c}{a+b+c}}}{\sqrt{\frac{a}{a+b+c}}+1}}
\frac{\log \left(i (-1)^{\left\lfloor \frac{\xi -1}{2}\right\rfloor } \sqrt{\frac{d}{b+c}}+i (-1)^{\xi -1} \sqrt{\frac{b+c+d}{b+c}}+t\right)}{t-i (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\eta +1} e^{i
(-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor } \tan ^{-1}\left(\frac{\sqrt{c}}{\sqrt{b}}\right)}} dt \\
&=&\frac{1}{4} \frac{1}{\sqrt{a b c d}} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^4
(-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\left\lfloor \frac{\xi -1}{2}\right\rfloor }
%
{\mathfrak F}^{(0,\frac{\sqrt{\frac{b+c}{a+b+c}}}{\sqrt{\frac{a}{a+b+c}}+1})}_{i (-1)^{\left\lfloor \frac{\xi -1}{2}\right\rfloor } \sqrt{\frac{d}{b+c}}+i (-1)^{\xi -1} \sqrt{\frac{b+c+d}{b+c}},-i (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\eta +1} e^{i
(-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor } \tan ^{-1}\left(\frac{\sqrt{c}}{\sqrt{b}}\right)}}
\end{eqnarray}
In the top line we substituted for $u=\sin(\theta) \sqrt{d/(b+c \sin(\theta)^2)}$. In the second line we substituted $u = \sqrt{d/(c+b)} \sin(\phi)$. In the third line we substituted $t=\tan(\phi/2)$. In the forth line we used partial fraction decomposition and properties of the logarithm. Finally in the fifth line we used the anti-derivative defined in $(ii)$.
Clear[F]; Clear[FF];
F[z_, a_, b_] :=
Log[a + z] Log[(b + z)/(-a + b)] + PolyLog[2, (a + z)/(a - b)];
FF[A_, B_, a_, b_] :=
Module[{result, ts, zs, zsp, zsm, eps = 10^(-15)},
(*This is Integrate[Log[z+a]/(z+b),{z,A,B}] where all a,b,A,
and B are complex. *)
result = F[B, a, b] - F[A, a, b];
ts = - (Im[(A + b) (Conjugate[b] - Conjugate[a])]/
Im[(B - A) (Conjugate[b] - Conjugate[a])]);
If[0 <= ts <= 1,
zsp = A + (ts + eps) (B - A);
zsm = A + (ts - eps) (B - A);
result += -F[zsp, a, b] + F[zsm, a, b];
];
result
];
{a, b, c, d} = RandomReal[{0, 3}, 4, WorkingPrecision -> 50];
NIntegrate[
Exp[-a/2 x^2 - b/2 y^2 - c/2 z^2 - d/2 w^2], {x, 0, Infinity}, {y, 0,
x}, {z, 0, y}, {w, 0, z}]
NIntegrate[
Sin[th]/Sqrt[a b d (b + c Sin[th]^2)] ArcTan[
Sin[th] Sqrt[d/(b + c Sin[th]^2)]], {th, 0, ArcTan[Sqrt[b/a]]}]
1/Sqrt[a ] NIntegrate[
u ArcTan[u] 1/((d - c u^2) Sqrt[d - (c + b) u^2]), {u, 0, Sqrt[
d/ (a + b + c)]}]
Sqrt[d]/(Sqrt[a ] (b + c))
NIntegrate[
Sin[phi] ArcTan[
Sqrt[d/(c + b)] Sin[phi]] 1/(d - c (d/(c + b) Sin[phi]^2)) , {phi,
0, ArcSin[ Sqrt[( c + b)/ (a + b + c)]]}]
- I 2/(Sqrt[a ] Sqrt[d])
NIntegrate[
t /(c (-1 + t^2)^2 + b (1 + t^2)^2) Log[(
1 + t^2 + 2 I Sqrt[d/(b + c)] t)/(
1 + t^2 - 2 I Sqrt[d/(b + c)] t)], {t, 0, Sqrt[(b + c)/(
a + b + c)]/(1 + Sqrt[a/(a + b + c)])}]
- I 2/(Sqrt[a ] Sqrt[d])
NIntegrate[
t /(c (-1 + t^2)^2 +
b (1 + t^2)^2) Log[((1/
2 (2 I Sqrt[d/(b + c)] - Sqrt[-4 - (4 d)/(b + c)]) +
t) (1/2 (2 I Sqrt[d/(b + c)] + Sqrt[-4 - (4 d)/(b + c)]) +
t))/((1/2 (-2 I Sqrt[d/(b + c)] - Sqrt[-4 - (4 d)/(b + c)]) +
t) (1/2 (-2 I Sqrt[d/(b + c)] + Sqrt[-4 - (4 d)/(b + c)]) +
t))], {t, 0, Sqrt[(b + c)/(a + b + c)]/(
1 + Sqrt[a/(a + b + c)])}]
1/Sqrt[a b c d] 1/4 NIntegrate[
Sum[(-1)^(Floor[(eta - 1)/2]) (-1)^
Floor[(xi - 1)/2] Log[
t + (-1)^Floor[(xi - 1)/2] I Sqrt[d/(b + c)] + (-1)^(xi - 1)
I Sqrt[( b + c + d)/(b + c)]]/(
t - (-1)^(1 + eta +
Floor[(eta - 1)/2]) I Exp[(-1)^(Floor[(eta - 1)/2]) I ArcTan[
Sqrt[c]/Sqrt[b]]]), {xi, 1, 4}, {eta, 1, 4}], {t, 0, Sqrt[(
b + c)/(a + b + c)]/(1 + Sqrt[a/(a + b + c)])}]
1/Sqrt[a b c d] 1/4 Sum[(-1)^(Floor[(eta - 1)/2]) (-1)^
Floor[(xi - 1)/2] FF[0, Sqrt[(b + c)/(a + b + c)]/(
1 + Sqrt[a/(
a + b + c)]), (-1)^Floor[(xi - 1)/2] I Sqrt[d/(b + c)] + (-1)^(
xi - 1) I Sqrt[( b + c + d)/(
b + c)], -(-1)^(1 + eta +
Floor[(eta - 1)/2]) I Exp[(-1)^(Floor[(eta - 1)/2]) I ArcTan[
Sqrt[c]/Sqrt[b]]]], {xi, 1, 4}, {eta, 1, 4}]
Update: As a sanity check look at the case $a=b=c=d=1$.
Then define:
\begin{eqnarray}
M1&:=&\left(
\begin{array}{cccc}
-1+\sqrt{3} & \sqrt{2} & \sqrt{2} & -1+\sqrt{3} \\
1 & \sqrt{2-\sqrt{3}} & \sqrt{2-\sqrt{3}} & 1 \\
\sqrt{2-\sqrt{3}} & 1 & 1 & \sqrt{2-\sqrt{3}} \\
\sqrt{2} & -1+\sqrt{3} & -1+\sqrt{3} & \sqrt{2} \\
\end{array}
\right)\\
M2&:=&\left(
\begin{array}{cccc}
\frac{1}{\sqrt{2}} & \frac{1}{2} \left(1+\sqrt{3}\right) & \frac{1}{2} \left(1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} \\
\frac{1}{\sqrt{2}} & \frac{1}{2} \left(-1+\sqrt{3}\right) & \frac{1}{2} \left(-1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} \\
\frac{1}{2} \left(-1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \frac{1}{2} \left(-1+\sqrt{3}\right) \\
\frac{1}{2} \left(1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \frac{1}{2} \left(1+\sqrt{3}\right) \\
\end{array}
\right)\\
A1&:=&\left(
\begin{array}{cccc}
-\frac{\pi }{6} & \frac{\pi }{12} & -\frac{\pi }{4} & 0 \\
\frac{5 \pi }{6} & \frac{\pi }{12} & \frac{5 \pi }{12} & -\frac{\pi }{3} \\
-\frac{5 \pi }{12} & \frac{\pi }{3} & -\frac{5 \pi }{6} & -\frac{\pi }{12} \\
\frac{\pi }{4} & 0 & \frac{\pi }{6} & -\frac{\pi }{12} \\
\end{array}
\right)\\
A2&:=&\left(
\begin{array}{cccc}
-\frac{\pi }{12} & \frac{\pi }{6} & -\frac{\pi }{6} & \frac{\pi }{12} \\
\frac{7 \pi }{12} & -\frac{\pi }{6} & \frac{\pi }{6} & -\frac{7 \pi }{12} \\
-\frac{\pi }{6} & \frac{7 \pi }{12} & -\frac{7 \pi }{12} & \frac{\pi }{6} \\
\frac{\pi }{6} & -\frac{\pi }{12} & \frac{\pi }{12} & -\frac{\pi }{6} \\
\end{array}
\right)
\end{eqnarray}
and we have
\begin{eqnarray}
I(1,1,1,1)=\frac{1}{4} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^4
(-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\left\lfloor \frac{\xi -1}{2}\right\rfloor }
\left(
Li_2(M1_{\xi,\eta}\exp(\imath A1_{\xi,\eta}))-
Li_2(M2_{\xi,\eta}\exp(\imath A2_{\xi,\eta}))
\right)
\end{eqnarray}
We have checked numerically that this quantity above coincides with $\pi^2/96$ to one hundred digits. It would be interesting to prove this analytically.
Best Answer
One suggestive principle is this. To evaluate the integral $\int_a^b f(x)\;dx$ by complex contour integral methods, it is (usually) required that $a$ and $b$ are "special" points for the function $f$: for example, poles or branch points.
For example, $$ \int_{-1}^1\frac{dx}{\sqrt{1-x^2}} \tag1$$ is a good candidate, since $\pm1$ are branch points for the integrand. But $$ \int_{1/2}^1 \frac{dx}{\sqrt{1-x^2}} \tag2$$ is not a good candidate, since the point $x=1/2$ is nothing special for that integrand. Probably evaluating $(2)$ is no easier than evaluating the indefinite integral $$ \int_{a}^1 \frac{dx}{\sqrt{1-x^2}} \tag2 $$ for all $a$ with $-1 < a < 1$.