Discrete version of pillbox function
I think you should first construct a smooth function of $r$ (where $r=\sqrt{x^2+y^2}$). For example, your original solution was based on the following function of one variable:

$$f_1(r) = \begin{cases} 1 &amp; r &lt; R \\
0 &amp; \text{otherwise} \end{cases}$$

but you might try using a "softer" step function such as:

$$f_2(r) = \begin{cases} 1 &amp; r&lt;R-\frac{s}{2} \\
0 &amp; r\ge R+\frac{s}{2} \\
\frac{1}{2} - \frac{r-R}{s} &amp; otherwise
 \end{cases}$$

and you can create smoother variants, such as $$f_3(r) = 3f_2(r)^2-2f_2(r)^3$$ since $f_3$ is simply formed by mapping the linear ramp function of $f_2$ to smoother cubic shape.

The transition width of the ramp, $s$, should be set to $\frac{1}{2}$ or $\frac{1}{\sqrt2}$ if the linear ramp, $f_2$, is used, and you may need a larger value of $s$ (say, $s=1$) if the cubic function, $f_3$ is used.

Note that, strictly speaking, what you are trying to do is sample the original pillar-box function at discrete sample-points, so if you think about the function prior to discrete sampling (let's call this the pre-sampling function, $f(x,y)$ being defined over the Cartesian plane, for all real $x$,$y$)  the problem with your original pre-sampling function was that it contained a step-discontinuity (a $0^{th}$order discontinuity), which contains large amounts of high-frequency energy (in 2D spatial frequency), so you need to remove some of this high-frequency energy.

The radial function, $f_2(r)$, removes the $0^{th}$order discontinuity, resulting in far less high-(spatial)-frequency energy in the pre-sampling signal (but the resulting signal still has a $1^{st}$order disconinuity, since it's first derivative is not continuous). Likewise, the cubic function, $f_3(r)$, has even more reduced high-frequency energy, since it has no $0^{th}$ or $1^{st}$order disconinuities (but it's second derivative is discontinuous).

In general, if you make sure your pre-sampling function has it's first $n-1$ derivatives being continuous (and it's $n^{th}$ is disconinuous), then the amplitude of it's high frequencies will fall of like $\left(\frac{1}{f}\right)^{n+1}$

I hope this helps...