Sum of the inverse squares of the hypotenuse of Pythagorean triangles What is the sum of the series

$$ S = \frac{1}{5^2} + \frac{1}{13^2} + \frac{1}{17^2} +
\frac{1}{25^2} + \frac{1}{29^2} + \frac{1}{37^2} + \cdots $$

where the sum is taken over all hypotenuse of primitive Pythagorean triangles.
By numerical computation, I found the sum to be $0.056840308812554488$ correct to $18$ decimal places. I would like to know if this sum has a closed form. 

Using the general formula for primitive Pythagorean triangles, $$ S = \sum_{r>s\ge 1, \\
\gcd(r,s)= 1}\frac{1}{(r^2 + s^2)^2} $$

Trivially, for all primitive and non primitive Pythagorean triangles, the sum will be $\zeta(2) = \pi^2/6$ times the corresponding sum for primitive Pythagorean triangles which turn out to be about $0.09349856033594433852$.
Motivation: We have equated the sum of the square of the sides of a right triangle to the square of the hypotenuse, so I was curious to know what the sum of the reciprocal of the square of the hypotenuse would be. Also since $\zeta(2)$ converges, and the density of hypotenuse is smaller than that of natural numbers, this sum must trivially converge.
Related question: What is the sum of the reciprocal of the hypotenuse of Pythagorean triangles?
 A: Let us use the notations from
http://mathworld.wolfram.com/DoubleSeries.html
We fix a positive definite binary quadratic form $q$ given by $q(m,n)=am^2+bmn+cn^2$, $a,b,c$ integers.
We use sumations over the index set $$J=\Bbb Z\times\Bbb Z-\{(0,0)\}\ .$$
We define
$$
\begin{aligned}
S(q;s) = S(a,b,c;s)          &=\sum_{(m,n)\in J} q(m,n)^{-s}=\sum_{(m,n)\in J} (am^2+bmn+cn^2)^{-s}\ ,\\
S_1(q;s) = S_1(a,b,c;s)      &=\sum_{(m,n)\in J} \color{blue}{(-1)^m}\; q(m,n)^{-s}\ ,\\
S_2(q;s) = S_2(a,b,c;s)      &=\sum_{(m,n)\in J} \color{blue}{(-1)^n}\; q(m,n)^{-s}\ ,\\
S_{12}(q;s) = S_{12}(a,b,c;s) &=\sum_{(m,n)\in J}\color{blue}{(-1)^{m+n}}\; q(m,n)^{-s}\ .
\end{aligned}
$$
The last three sums are "twisted versions" of the first sum,
the "twist" occurs by using a character for the first parameter, for the second one, for both.
In our case, $q(m,n)=m^2 +n^2$, and $(a,b,c)=(1,0,1)$, we have a symmetric case (w.r.t. the exchange $a\leftrightarrow c$).
We will drop $q$ below from notations in $S_?(q,s)$, since we use only the above quadratic form $q$. I decided during the  edit operation that should bring us quick to numbers we can compute that it is better for the check to introduce the versions $S^+$ for all sums, where the plus index indicates a further restriction to $(m,n)\in J$ with 
$$(+)\qquad m,n>0\ .$$ 
From loc. cit. we extract the following relations:
$$
\begin{aligned}
S(s) &= \sum_{(m,n)\in J}(m^2+n^2)^{-s} 
\\
&= 4\beta(s)\;\zeta(s)\ ,\\
%S_1(s) =S_2(s) &= \sum_{(m,n)\in J}(-1)^m(m^2+n^2)^{-s} =\dots
%= 2^{-s}b_2(2s) = -2^{-s}\cdot 4\beta(2s)\; \eta(2s) 
%\ ,
%\\
S_{12}(s) &= \sum_{(m,n)\in J}(-1)^{m+n}(m^2+n^2)^{-s} 
\\
&= 
-4 \beta(s) \;\eta(s)=-4\beta(s) \;(1-2^{1-s})\; \zeta(s)\ .
\\[2mm]
&\qquad\text{ Then the plus versions are:}
\\[3mm]
S^+(s) &= \beta(s)\;\zeta(s) - \zeta(2s)\ ,
\\
-S_{12}^+(s) &= \beta(s)\;\eta(s) - \eta(2s)
\\
&=
\beta(s)\;(1-2^{1-s})\zeta(s) - (1-2^{1-2s})\zeta(2s)\ ,
\\
&\qquad \text{ which gives}
\\
S^+(s)-S_{12}^+(s)
&=2\beta(s)\;(1-2^{-s})\zeta(s) - 2(1-2^{-2s})\zeta(2s)\ .
\end{aligned}
$$
Now let us search for a linear combination of the above sums
that correspond to summing $q(m,n)^{-s}$ over the set of $K$ of all $(m,n)$ with positive (components with) different parity.
This is
$$\frac 12(\ S^+(s)-S^+_{12}(s)\ )\ .$$
So far can we write:
$$
\begin{aligned}
&\beta(s)\;(1-2^{-s})\zeta(s) - (1-2^{-2s})\zeta(2s)
\\ 
&\qquad=\frac 12(\ S^+(s)-S^+_{12}(s)\ ) \\
&\qquad=\sum_{\substack{(m,n)\in K\\m,n> 0}} q(m,n)^{-s}\\
&\qquad=2\sum_{\substack{(m,n)\in K\\m>n> 0}} q(m,n)^{-s}\\
&\qquad=2\sum_{\substack{(m,n)\in K\\m>n> 0\\ d=(m,n)\text{ odd}}} q(m,n)^{-s}\qquad\text{ and with }M=m/d,\ N=n/d\\
&\qquad=2\sum_{d>0\text{ odd}}d^{-2s}
\sum_{\substack{(M,N)\in K\\M>N> 0\\ (M,N)=1}} q(M,N)^{-s}\\
&\qquad=
2(1-2^{-2s})\; \zeta(2s)\;
\color{blue}{
\sum_{\substack{(M,N)\in K\\M>N> 0\\ (M,N)=1}} q(M,N)^{-s}
}
\ .
\end{aligned}
$$
The isolated sum in the last expresson is the sum we need, let us take it for $s=2$.
The value we obtain is:
$$
\color{brown}{
\frac{\beta(2)\;\zeta(2)}{2(1+2^{-2})\zeta(4)}
-\frac 12\ = \frac{6C}{\pi^2} - \frac 12.}
$$
$$
\color{brown}{
\frac{\beta(2)\zeta(2)}{2(1+2^{-2})\zeta(4)} - \frac{1}{2}
= \frac{6C}{\pi^2} - \frac 1 2.}
$$
where $C$ is the Catalan constant.
Numerically:
sage: E = catalan * zeta(2) / 2 / (1+2^-2) / zeta(4) - 1/2 
sage: E.n()
0.0568403090661582

