It is an old question but I recently solved a programming problem regarding this question, andthe article was shared above but the method was not explicitly shared here.
Following this article: Pythagorean triangles with legs less than n there is quite a nice solution
We know that each primitive Pythagorean triple can be represented by \begin{align*} a &= x^2 - y^2 \\ b &= 2xy \\ c &= x^2 + y^2 \end{align*} Following the paper lets define a few sets. \begin{align*} R(n) &= \{(x, y) \in \mathbb{Z} \times \mathbb{Z} \mid x^2 + y^2 \leq n, 0 < x < y\}\\ P(n) &= \{(x, y) \in R(n) \mid \gcd(x, y) = 1, x + y \equiv 1 \pmod{2}\} \\ Q(n) &= \{(x, y) \in R(n) \mid \gcd(x, y) = 1\} \end{align*} Note that $|P(n)|$ is exactly the amount of Pythagorean triples with hypotenuse less than or equal to $n$, which is what we are trying to find. Now according to the paper Lemma 3 states \begin{equation*} |P(n)| = \sum_{k \geq 0} (-1)^k|Q(\frac{n}{2^k})| \end{equation*} Note that we changed the region but the proof is the same.
Now all that is left is to calculate $Q(n)$. I define a few equations which are identical to the article, but easier to read. \begin{align*} L(t) &= |\{(x, y) \in R(t^2)\}| = |R(t^2)| \\ L'(t) &= |\{(x, y) \in R(t^2) \mid \gcd(x, y) = 1\}| = |Q(t^2)| \end{align*} Following the exact same proof as the article we have \begin{align*} L(t) &= \sum_{1 \leq d \leq t} L'(\frac{t}{d}) \\ L'(t) &= \sum_{1 \leq d \leq t} \mu(d) L(\frac{t}{d}) \end{align*} Therefore we have, \begin{align*} |Q(n)| &= L'(\sqrt{n}) \\ &= \sum_{1 \leq d \leq \sqrt{n}} \mu(d) L(\frac{\sqrt{n}}{d}) \\ &= \sum_{1 \leq d \leq \sqrt{n}} \mu(d) |R(\frac{n}{d^2})| \end{align*}
Now, using brute force a computer can easily find all the lattice points in $R(n)$, and with a good Möbius sieve we can get a nice easy solution.
I have a python package where I add functions like these if anyone is interested to see the code.