finding the circle with maximal rainfall

Question

Let's say that I have the rainfall in millimeters for each point on a 2-dimensional map.

I want to find the (x,y) coordinate of a circle of radius R that has received the maximum amount of rain.

Something like (for a square actually, probably in polar coordinate it will be easy to express a circle)

$$mm(x,y,R)=\int\limits_{x-R}^{x+R}\int\limits_{y-R}^{y+R}f(s,t)dsdt$$

and I want to find the point (x,y) that gives the maximum value of mm, given R fixed.

I'd like to frame the problem from a mathematical point of view.

From a computational point of view, I'll be dealing with a discrete function f as I'll have values for rainfall on a discrete grid.

I've seen some algorithms trying to calculate the mm by actually calculating the integral over many values of (x,y), but that looks to me like a brute force approach.

Is there any other mathematical way of solving this? I wouldn't mind if I find a square or a circle, the concept will be the same

Answer 1

I don't expect a closed-form solution exists that will work for you, since it sounds like you are working with experimental data, meaning the function $f(x,y)$ will be discrete and probably discontinuous. If you had a continuous function $f(x,y)$, then you could use partial derivatives to find optimum values.

But for the computational problem, I don't think it's as bad as you think. Once you compute the sum of the rainfall at all points in a square of width $2W+1$ at a point $(x_0,y_0)$, you can easily compute the sum of rainfall at a point $(x_0,y_0+1)$, because you just need the difference between the two. That is, the new square includes all the same points as the old square, except the new square also has $(x,y_0+1+W)$ for $x$ between $x_0-W$ and $x_0+W$, and the old square has points $(x,y_0-W)$ for $x_0-W\leq x\leq x_0+W$ which are not in the new square. Thus, all you have to do is take the old sum, add the values of the $2W+1$ new points, then subtract the values of the $2W+1$ points which were removed.

So, once you compute $mm(x_0,y_0,2W+1)$, which costs $O(W^2)$ additions, you only need to compute another $2(2W+1)$ additions to find $mm(x_0,y_0+1,2W+1)$. If your entire map has $L$ points, this is an $O(W^2+WL^2)$ algorithm to find the optimum.

If $O(WL^2)$ is still too much computation, you could store the values of $$ B(x_0,y_0)=\sum_{x=x_0-W}^{x_0+W}f(x,y_0)$$ for all the values of $y$, and then you can compute $B(x_0+1,y_0)$ with only 2 additions, as

$$ B(x_0+1,y_0)=B(x_0,y_0) + f(x_0+W+1,y_0) - f(x_0-W,y_0)$$

Overall this would require $O(L)$ memory but only $O(L^2)$ computation.

For a circle it will be the same principle, but the boundary of added/removed points will be slightly more difficult to compute.

Answer 2

Convolution with respect to the indicator function of the circle, using FFT

background

• Put your data on a grid of size $n\times n$, make sure that your data is in the center region of the grid, and that you have "zero padding" around the data (that is, towards the edge of the grid, there is a region with zero data). call this grid $D$ (for data).
• On a new grid (also $n\times n$), draw a circle of radius $R$ around the middle point of the grid. Put ones inside the circle, and zeros outside. Call this grid $I$ (for indicator).
• Calculate the 2d Fourier transform of each of the grids, call then $\hat D$, $\hat I$.
• Calcluate the product of the Fourier transforms (pixel by pixel): $\hat C = \hat D \cdot \hat I$.
• The inverse Fourier transform of $\hat C$ (call it $C$) is (up to discretization errors and normalization) is the convolution of $I$ with $D$, and is the object that you have called $mm$.

Note: Calculating an Fourier transform of the grid takes the order of $(n\times n) \log(n)$ operations. This is probably as efficient as you can get.