For any non-negative integer $k$, $x^2+y^2=5^k$ has $4(k+1)$ integer solutions.
Suppose that a polygon has integer coordinates for all of its vertices. Let $i$ be the number of integer points interior to the polygon, and let $b$ be the number of integer points on its boundary (including both vertices and points along the sides). Then the area $A$ of this polygon is:
$$A = i + \frac{b}{2} - 1$$
For a radius $r$, the Gauss Circle Problem determines the number of lattice points inside a circle $x^2 + y^2 \leq r^2$ as
$$N(r) = 1 + 4\sum_{i=0}^\infty\left(\left\lfloor \frac{r^2}{4i+1}\right\rfloor - \left\lfloor \frac{r^2}{4i+3}\right\rfloor\right)$$
Show $\pi \approx \frac{4\sum_{i=0}^\infty\left(\left\lfloor \frac{5^k}{4i+1}\right\rfloor - \left\lfloor \frac{5^k}{4i+3}\right\rfloor\right) - 2(k+1) }{5^k}$
Proof
When the lattice polygon closely approximates the circle, their areas $A$ are nearly equal:
$$\pi r^2 \approx i + \frac{b}{2} - 1$$
So
$$\pi \approx \frac{i + \frac{b}{2} - 1}{r^2}$$
Let $b \approx 4(k+1)$, $r=\sqrt{5^k}$, $i \approx N(\sqrt{5^k}) - 4(k+1)$
\begin{align} \pi &\approx \frac{i + \frac{b}{2} - 1}{r^2} \\ &= \frac{N(\sqrt{5^k}) - 4(k+1) + \frac{4(k+1)}{2} - 1}{5^k} \\ &= \frac{N(\sqrt{5^k}) - 2k - 3}{5^k} \\ &= \frac{4\sum_{i=0}^\infty\left(\left\lfloor \frac{5^k}{4i+1}\right\rfloor - \left\lfloor \frac{5^k}{4i+3}\right\rfloor\right) - 2(k+1) }{5^k} \end{align}
Question
Is that a valid estimate for $\pi$?
For $k=8$ and $N=250$, Wolfram gives $\pi \approx 3.1395072$.
Notes
When joining the $b=4(k+1)$ boundary lattice points to form the polygon, there are possibly other integer points in-between
$$\text{gcd}(|x_2 - x_1|, |y_2 - y_1|) - 1$$
so it's a rough approximation of $b$.
For example, a circle $x^2+y^2=5^2$ with boundary points $(10,5),(5,10)$ has four other points in-between $(9, 6)$, $(8, 7)$, $(7, 8)$, and $(6, 9)$. Does there exist a formula to compute the exact number of boundary points?
Update
To hopefully improve the numerical issues for $k>9$ as noted in the answer, I've submitted this question
Currently, the $\pi$ estimate uses $b \approx 4(k+1)$ which isn't very accurate. For example, $k=14$ has $b\approx60$ points, however, the true number of boundary points is $b=117540$.
For an N-polygon, the precise total number of boundary points is
$$b = \underbrace{4(k+1)}_{\text{Vertices}} + \underbrace{\sum_{i=0}^{N-1}\left(\gcd(|x_{(i+1)} - x_{i}|, |y_{(i+1)} - y_{i}|) - 1\right)}_{\text{Total Points In-between Vertices}}$$
So, perhaps a better $\pi$ estimate is
$$\pi \approx \frac{i + \frac{b}{2} - 1}{r^2}$$
where
$b = 4(k+1)+\sum_{i=0}^{N-1}\left(\gcd(|x_{(i+1)} - x_{i}|, |y_{(i+1)} - y_{i}|) - 1\right)$
$r=\sqrt{5^k}$
$i = N(\sqrt{5^k}) - b$
If we look at
$$a_k=\sum_{i=0}^\infty\left(\left\lfloor \frac{5^k}{4i+1}\right\rfloor - \left\lfloor \frac{5^k}{4i+3}\right\rfloor\right)$$ they generate the sequence $$\{1,5,20,100,490,2452,12269,61360,306776,1533987,\cdots\}$$ Now, for your formula $$b_k=\frac{4 a_k-2(k+1)}{5^k}$$ generate the sequence $$\left\{2,\frac{16}{5},\frac{74}{25},\frac{392}{125},\frac{78}{25},\frac{9 796}{3125},\frac{49062}{15625},\frac{245424}{78125},\frac{1227086}{390 625},\frac{6135928}{1953125},\cdots\right\}$$
Converted to exact decimal numbers $$\left( \begin{array}{cc} k & b_k \\ 0 & 2.0000000000 \\ 1 & 3.2000000000 \\ 2 & 2.9600000000 \\ 3 & 3.1360000000 \\ 4 & 3.1200000000 \\ 5 & 3.1347200000 \\ 6 & 3.1399680000 \\ 7 & 3.1414272000 \\ 8 & 3.1413401600 \\ 9 & 3.1415951360 \\ \end{array} \right)$$
For $k>9$, a lot of serious numerical issues.