We can write
$$
\frac{f(n)}{n^2} = \frac{\delta_n}{n} + \frac{1}{n} \sum_{k=0}^{n-1} \left(1-\sqrt{1-k^2/n^2}\right),
$$
where $0 \leq \delta_n \leq 1$. When $n$ is not the hypotenuse of a pythagorean triple we have the formula
$$
\delta_n = 1 - \frac{1}{n} - \frac{1}{n} \sum_{k=0}^{n-1} \left\{n - \sqrt{n^2-k^2}\right\}, \tag{$*$}
$$
with $\{x\}$ denoting the fractional part of $x$.
It can be deduced from this answer of mine that
$$
\frac{1}{n} \sum_{k=0}^{n-1} \left(1-\sqrt{1-k^2/n^2}\right) = \int_0^1 \left(1-\sqrt{1-x^2}\right)dx - \frac{1}{2n} + O\left(n^{-3/2}\right),
$$
so at least we know that
$$
\frac{f(n)}{n^2} = 1 - \frac{\pi}{4} + \frac{2\delta_n - 1}{2n} + O\left(n^{-3/2}\right).
$$
The behavior of $\delta_n$ is harder to get a handle on. Numerically it seems to tend to $1/2$ as $n \to \infty$, as can be seen from the following plot.

Because of this I would suspect that
$$
\frac{f(n)}{n^2} = 1 - \frac{\pi}{4} + o(n^{-1})
$$
as $n \to \infty$. Unfortunately without more information about $\delta_n$ this can't be made more precise.
The behavior $\delta_n \to 1/2$ is what we would expect if the summands $\left\{n - \sqrt{n^2-k^2}\right\}$ from $(*)$ were roughly uniformly distributed over the interval $[0,1]$. Perhaps this is the case in some specific sense, and I would be interested if someone could say something about it.