It seems that this has been answered here by S. Mustonen: PointsInGrid.pdf
There $f(n,m)$ is denoted by $L_n(m)$ and a formula would be $$L_n(m)=\frac{1}{2}[f(m,n+1)-2f(m,n)+f(m,n-1)]$$ for $$f(m,k)=\sum_{\substack{-m<kx<m\\-m<ky<m\\(x,y)=1}}(m-|kx|)(m-|ky|)$$