Compute the max sum circular area

Question

I have an n by n matrix of integers and I want to find the circular area, with origin at the top left corner, with maximum sum. Consider the following grid with a circle imposed on it.

This is made with:

import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import numpy as np
plt.yticks(np.arange(0, 10.01, 1))
plt.xticks(np.arange(0, 10.01, 1))
plt.xlim(0,9)
plt.ylim(0,9)
plt.gca().invert_yaxis()
# Set aspect ratio to be equal
plt.gca().set_aspect('equal', adjustable='box')
plt.grid()   
np.random.seed(40)
square = np.empty((10, 10), dtype=np.int_)
for x in np.arange(0, 10, 1):
    for y in np.arange(0, 10, 1):
        plt.scatter(x, y, color='blue', s=2, zorder=2, clip_on=False) 
for x in np.arange(0, 10, 1):
    for y in np.arange(0, 10, 1):
        value = np.random.randint(-3, 4)
        square[int(x), int(y)] = value
        plt.text(x-0.2, y-0.2, str(value), ha='center', va='center', fontsize=8, color='black')

r1 = 3
circle1 = Circle((0, 0), r1,  color="blue", alpha=0.5, ec='k', lw=1)
plt.gca().add_patch(circle1)

In this case the matrix is:

[[ 3  0  2 -3 -3 -1 -2  1 -1  0]
 [-1  0  0  0 -2 -3 -2  2 -2 -3]
 [ 1  3  3  1  1 -3 -1 -1  3  0]
 [ 0  0 -2  0  2  1  2  2 -1 -1]
 [-1  0  3  1  1  3 -2  0  0 -1]
 [-1 -1  1  2 -3 -2  1 -2  0  0]
 [-3  2  2  3 -2  0 -1 -1  3 -2]
 [-2  0  2  1  2  2  1 -1 -3 -3]
 [-2 -2  1 -3 -2 -1  3  2  3 -3]
 [ 2  3  1 -1  0  1 -1  3 -2 -1]]

When the circle has radius 3, there are 11 points in the grid within the circle. As the radius increases, more and more points fall into the circle.

I am looking for a fast way to find a radius which maximizes the sum of the integers of grid points within it. The radius will not be unique so any one that maximizes the sum is ok. I will ultimately want to do this with much larger matrices.

This question is related but I am not sure how to extend it to my new question.

Answer 1

Same O(n^2 log n) method as others: accumulate weight sums by distance, then compute cumulative sums. But in pure python which should be trivial to translate to C++. Also PyPy may run surprisingly fast.

from collections import Counter

g = [[-3, 2, 2],
     [ 2, 0, 3],
     [-1, -2, 0]]
n = 3

sum_dist = Counter()
for i in range(n):
    for j in range(n):
        dist = i**2 + j**2
        sum_dist[dist] += g[i][j]

sorted_dists = sorted(sum_dist.keys())
for i in range(1, len(sorted_dists)):
    sum_dist[sorted_dists[i]] += sum_dist[sorted_dists[i-1]]

print(sum_dist)
print(max(sum_dist, key=sum_dist.get))