This is an image processing question.

In 2-D Euclidean space a pillbox function, $f(x)$ is given by:

$$f(x,y) = 1\text{ if } x^2 + y^2 < R, 0 \text{ otherwise } $$

$R$ is come constant.

If I try to approximate the top hat when I'm looking at an image matrix of points, a naive algorithm would make some of the image pixels 1 inside a given radius and the others 0 outside that radius.

However, if I'm looking at a 4 x 4 grid of points, the shape I'd approximate would look like a square if I tried to naively create a discrete top-hat.

How would I create a discrete top-hat function that accurately approximates it's continuous counterpart, even when my image matrix is small?

EDIT: I believe the comment about the edge image pixels having a value equal to the area of the pixel contained within the pillbox is correct. However, implementing this algorithm seems difficult. A reference would be appreciated.

EDIT 2: Matlab has an algorithm for computing a discrete pillbox function in it's fspecial function (disk option). This is the algorithm I'm after, though I still need a reference or an explanation of the geometry, as the function is quite complicated.

  • $\begingroup$ I'm debating assigning a 1 to entirely inside pixels, 0 to entirely outside pixels and an intermediate value to partially in, partially out pixels. $\endgroup$ – ncRubert Sep 18 '14 at 21:51
  • 1
    $\begingroup$ I would set the coefficient equal to the area of the pixel which is enclosed by the circle, so it would be 1 for the inner pixels and some value < 1 for pixels that the circle intersects. The sum of all the coefficients would then be equal to the area of the circle. $\endgroup$ – Paul R Sep 19 '14 at 11:43

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 & r < R \\ 0 & \text{otherwise} \end{cases}$$

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

$$f_2(r) = \begin{cases} 1 & r<R-\frac{s}{2} \\ 0 & r\ge R+\frac{s}{2} \\ \frac{1}{2} - \frac{r-R}{s} & 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...

| improve this answer | |

Building on your comment about assigning 1 to pixels inside the circle and 0 to pixels entirely outside, with values in between for ones that are on the edge, I would try this:

Create you circle as 10N*10N instead of N*N. Maybe even use a higher value in place of 10. Draw your circle with your naive algorithm - 1 inside, 0 outside. Now scale this image down from 10N*10N to N*N. This will automatically give you fractional values for the pixels that are on the border at the lower resolution. You could also at this point split it into 0 and 1 again with a threshold function. Doing that with 4*4 gives you something like this (colors inverted:)

enter image description here

| improve this answer | |

You will need to use continuous values if you want a better approximation on small images. As Paul R mentioned a natural choice here would be to choose the ratio of each pixel's area inside the circle to the area of the pixel.

You could perform an integral to give you the area of a pixel inside the circle, however, a Riemann sum is much less tedious. Partition each pixel into a grid of N x N equal-sized squares. If the center of a square is inside the circle then its area is added to the sum.

The discrete pixel coordinates represent the center of a $1 \times 1$ square in the Cartesian coordinate system. The origin of the coordinate system is chosen to be the center of the circle.

Coordinate system: pillbox-3x3-grid
Labeled diagram: pillbox-3x3-labelled

Everything is labeled using $(x,y)$ Cartesian coordinates to make the math consistent, however, since a circle is symmetric you don't need to worry about it. The relationship is of course that $(x,y)$ coordinates correspond to $[-y,x]$ pixel coordinates.

From the diagram you can set up the appropriate integrations to find out the areas within each pixel's box. Then you just need to normalize.

If you need to create pillboxes in code then it would be simpler to use a Riemann sum albeit at the cost of computational time. From a signal processing standpoint this Riemann sum would be equivalent to:

$Upsample \rightarrow Calculation \rightarrow Downsample$

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.