What is the best way to check all pixels within certain radius?

Question

I'm currently developing an application that will alert users of incoming rain. To do this I want to check certain area around user location for rainfall (different pixel colours for intensity on rainfall radar image). I would like the checked area to be a circle but I don't know how to do this efficiently.

Let's say I want to check radius of 50km. My current idea is to take subset of image with size 100kmx100km (user+50km west, user+50km east, user+50km north, user+50km south) and then check for each pixel in this subset if it's closer to user than 50km.

My question here is, is there a better solution that is used for this type of problems?

Answer 1

If the occurrence of the event you are searching for (rain or anything) is relatively rare, then there's nothing wrong with scanning a square or pixels and then, only after detecting rain in that square, checking whether that rain is within the desired 50km circle. Note that the key point here is that you don't need to check each pixel of the square for being inside the circle (that would be very inefficient), you have to search for your event (rain) first and only when you found it, check whether it falls into the 50km circle. To implement this efficiently you also have to develop some smart strategy for handling multi-pixel "stains" of rain on your image.

However, since you are scanning a raster image, you can easily implement the well-known Bresenham circle algorithm to find the starting and the ending point of the circle for each scan line. That way you can easily limit your scan to the desired 50km radius.

On the second thought, you don't even need the Bresenham algorithm for that. For each row of pixels in your square, calculate the points of intersection of that row with the 50km circle (using the usual schoolbook formula with square root), and then check all pixels that fall between these intersection points. Process all rows in the same fashion and you are done.

P.S. Unfortunately, the Wikipedia page I linked does not present Bresenham algorithm at all. It has code for Michener circle algorithm instead. Michener algorithm will also work for circle rasterization purposes, but it is less precise than Bresenham algorithm. If you care for precision, find a true Bresenham on somewhere. It is actually surprisingly diffcult to find on the net: most search hits erroneously present Michener as Bresenham.

Answer 2

Maybe you can try something that will speed up your algorithm.

In brute force algorithm you will probably use equation:

(x-p)^2 + (y-q)^2 < r^2
(p,q) - center of the circle, user position
r - radius (50km)

If you want to find all pixels (x,y) that satisfy above condition and check them, your algorithm goes to O(n^2)

Instead of scanning all pixels in this circle I will check only only pixels that are on border of the circle.

In that case, you can use some more clever way to define circle.

x = p+r*cos(a)
y = q*r*sin(a) 
a - angle measured in radians [0-2pi]

Now you can sample some angles, for example twenty of them, iterate and find all pairs (x,y) that are border for radius 50km. Now check are they on the rain zone and alert user.

For more safety I recommend you to use multiple radians (smaller than 50km), because your whole rain cloud can be inside circle, and your app will not recognize him. For example use 3 incircles (r = 5km, 15km, 30km) and do same thing. Efficiency of this algorithm only depends on number of angles and number of incircles.

Pseudocode will be:

checkRainDanger()
  p,q <- position
  radius[] <- array of radii
  for c = 1 to length(radius)
    a=0
    while(a<2*pi)
      x = p + radius[c]*cos(a)
      y = q + radius[c]*sin(a)
      if rainZone(x,y)
        return true
      else          
        a+=pi/10
    end_while
  end_for
  return false //no danger

Answer 3

r2=r*r
for x in range(-r, +r):
  max_y=sqrt(r2-x*x)
  for y in range(-max_y, +max_y):
    # x,y is in range - check for rain

Answer 4

There is, you can modify the midpoint circle algorithm to give you an array of for each y, the x coordinate where the circle starts (and ends, that's the same thing because of symmetry). This array is easy to compute, pseudocode below.

Then you can just iterate over exactly the right part, without checking anything.

Pseudo code:

data = new int[radius];
int f = 1 - radius, ddF_x = 1;
int ddF_y = -2 * radius;
int x = 0, y = radius;
while (x < y)
{
    if (f >= 0)
    {
        y--;
        ddF_y += 2; f += ddF_y;
    }
    x++;
    ddF_x += 2; f += ddF_x;
    data[radius - y] = x; data[radius - x] = y;
}