This question asks how to construct a rasterized two-dimensional kernel for a circular uniform distribution.
The underlying mathematical problem is this: given a planar circle $R$ and a rectangle $X,$ find the area of $X\cap R$. The kernel's values will be these areas divided by the area of $R$ (so that it integrates to unity).
This is an elementary integration problem but it's messy. Borrowing statistical methods for describing probability distributions in two dimensions, we can construct the solution out of simpler ones. These are:
Shift and scale. Choose the circle's center as the origin and choose units of measurement in which it has a unit radius.
Create a "distribution function." This function $F$ of two arguments $x$ and $y$ returns the area of the circle located to the left of $x$ and beneath $y$. Using it we may express the area of the circle covered by the rectangle with lower left corner $(x_0,y_0)$ and upper right corner $(x_1,y_1)$ as
$$f(x_0,y_0; x_1,y_1) = F(x_1,y_1) - F(x_0,y_1) - F(x_1,y_0) + F(x_0,y_0).$$
Exploit symmetries. By negating $x$ and or $y$ and swapping $x$ with $y$ as necessary we can relate $F(x,y)$ to the value of $F(u,v)$ where $u \ge v\ge 0$.
The computation requires the "marginal" distributions $F_X(x)$, the area of the circle to the left of $x$, and $F_Y(y)$, the area beneath $y$. By symmetry $F_X=F_Y$.
By integration (or elementary trigonometry) we find that for $-1 \le x \le 1$,
$$F_X(x) = \pi/2 + \arcsin(x) + x\sqrt{1-x^2}.$$
(The pieces are twice the areas of a quarter circle, a sector, and a triangle, respectively.) This function rises from $0$ to $\pi$.
I leave the details of (1) - (3) to the interested reader (they can be found in the implementation of pucirc below). They can all be worked out by carving $R$ along the vertical lines at $0$ and $x$ and the horizontal lines at $0$ and $y$. There are really only two distinct cases to consider: $(x,y)$ lies inside the circle or outside it (that is, $x^2+y^2$ is less than $1$ or greater than $1$). All the resulting pieces can be expressed in terms of $F_X$, simple constants, and simple functions (like $xy$, which is the area of the rectangle from the origin to the point $xy$).
To illustrate, here are the overlap areas in a unit grid for a central circle of radius $2.8$, as requested in the question.
The circle's boundary is overlaid in orange. The areas of each cell are posted (rounded to $0.001$). They were found by applying $f$ repeatedly to each of the cells: see the "Application" section in the code for details.
As a quick check, the last line of this code sums all the cell areas and compares them to the theoretical area of the circle. The resulting discrepancy should be zero (and it is in this example). Several tests (which graph various functions) have been left in but commented out.
#
# Area of a circle of radius `r`, centered at `origin`,
# to the left of `x` and below `y`.
#
pucirc <- function(x, y, radius=1, origin=c(0,0)) {
#
# Helper function.
#
xsqrt <- function(x) {
y <- pmax(-1, pmin(1, x))
y * sqrt(1 - y^2) + asin(y)
}
#
# Area of the unit circle to the left of `x`.
#
pucirc_margin <- function(x) {
ifelse(x >= 1, pi, ifelse(x <= -1, 0, pi/2 + xsqrt(x)))
}
# curve(pucirc_margin(x), from=-1.5, to=1.5)
#
# Area of the unit circle below `y` and to the left of `x` for
# 0 <= x,y
#
pucirc_ <- function(x, y) {
u <- pmax(x, y)
v <- pmin(x, y)
ifelse(u >= 1, pucirc_margin(v),
ifelse(u^2+v^2 > 1, pucirc_margin(u) + pucirc_margin(v) - pi,
pucirc_margin(u)/2 + pucirc_margin(v)/2 - pi/4 + u*v))
}
# curve(pucirc_margin(x), from=0, to=1.5, ylim=c(0, pi), lwd=2)
# curve(pucirc_(x, 0), add=TRUE, col="Gray", lty=3, lwd=2)
# for (y in c(1/3, 0.9, 1)) {
# curve(pucirc_(x, y), from=0, to=1.5, add=TRUE, col=hsv(y * 2/3), lwd=2)
# }
#
# Area of the unit circle below `y` and to the left of `x` for
# any `x` and `y`.
#
pucirc0 <- function(x, y) {
ifelse(x < 0 & y < 0, pi - pucirc_margin(-x) - pucirc_margin(-y) + pucirc_(-x, -y),
ifelse(x < 0 & y >= 0, pucirc_margin(y) - pucirc_(-x, y),
ifelse(x >= 0 & y < 0, pucirc_margin(x) - pucirc_(x, -y),
pucirc_(x, y))))
}
# curve(pucirc_margin(x), from=-1.5, to=1.5, ylim=c(0, pi), lwd=2)
# curve(pucirc0(x, 0), add=TRUE, col="Gray", lty=3, lwd=2)
# for (y in c(-1, -1/2, 1/3, 0.9, 1)) {
# curve(pucirc0(x, y), add=TRUE, col=hsv((y+1) * 1/3), lwd=2)
# }
pucirc0((x-origin[1])/radius, (y-origin[2])/radius) * radius^2
}
#
# Area of a circle overlapping a rectangle with lower left corner `ll`
# and upper right corner `ur`.
#
pucirc_rect <- function(ll, ur, radius=1, origin=c(0,0)) {
pucirc(ur[1], ur[2], radius, origin) -
pucirc(ur[1], ll[2], radius, origin) -
pucirc(ll[1], ur[2], radius, origin) +
pucirc(ll[1], ll[2], radius, origin)
}
#------------------------------------------------------------------------------#
#
# Application: construct a rasterized circular uniform kernel.
#
library(data.table) # Fast calculation
library(ggplot2) # Graphical display
radius <- 2.8
rmax <- floor(2*radius) / 2 # Smallest rectangle coordinate
rmin <- floor(-2*radius-1) / 2 # Larges rectangle coordinate
#
# Create a table of all lower-left cell coordinates.
#
X <- as.data.table(expand.grid(i=seq(rmin, rmax, by=1), j=seq(rmin, rmax, by=1)))
#
# Compute the overlap areas for the cells.
#
f <- Vectorize(function(i,j) pucirc_rect(c(i,j), c(i,j)+1, radius=radius))
invisible(X[, Overlap := f(i,j)])
#
# Create data to plot the circle itself.
#
theta <- seq(0, 2*pi, length.out=361)
P <- data.table(x=cos(theta)*radius, y=sin(theta)*radius)
#
# Plot everything.
#
g <- ggplot(X, aes(i,j)) +
geom_raster(aes(fill=Overlap, alpha=Overlap), hjust=1, vjust=1) +
geom_text(aes(label=round(Overlap, 3)), nudge_x=1/2, nudge_y=1/2, size=3) +
geom_path(data=P, aes(x, y), size=1.5, alpha=2/3, color="Orange") +
guides(alpha="none", fill="none") +
coord_fixed(ratio=1)
print(g)
#
# Check.
#
sum(X$Overlap) / (pi * radius^2) - 1 # Should be zero.
