Thanks, helps me see what you’re aiming for. But otherwise that was a great idea to paste a minimal working example. I think this algorithm should scale pretty well to 3D. The example shows how you would do it for a single grid size, then just wrap into a function taking n as argument:

```
# start with example code from above
x <- rnorm(500, 50, 10)
y <- rnorm(500, 50, 10)
plot(x, y, xlim=c(0, 100), ylim=c(0, 100))
abline(h=seq(0, 100, 10), col="grey")
abline(v=seq(0, 100, 10), col="grey")
# now define the number of bins we want
n <-10
output <- matrix(NA, 10, 10)
# now define a function that takes original data and rescales it onto axis 0...n
rescale <- function(x, x0, xm, n) {
(x - x0)/(xm - x0)*n
}
xi <- rescale(x, 0, 100, 10)
yi <- rescale(y, 0, 100, 10)
plot(xi, yi, xlim=c(0, 10), ylim=c(0, 10), col="grey", pch=16)
points(xi, yi)
abline(h=seq(0, 10, 1), col="grey")
abline(v=seq(0, 10, 1), col="grey")
# now we can do something similar, but this time taking the index part as the exact bin
xii <- ceiling(rescale(x, 0, 100, 10))
yii <- ceiling(rescale(y, 0, 100, 10))
points(xii-0.5, yii-0.5, col="red", pch=16, cex=0.5)
# finally, use table to count the number in occupied cells
(ret <- table(xii, yii))
# Total number occupied cells
sum(ret > 0)
# And just to be sure let's plot our results
x2 <- as.numeric(rownames(ret))[row(ret)]
y2 <- as.numeric(colnames(ret))[col(ret)]
text(x2-0.4, y2-0.4, as.vector(ret), col="red")
```

A benefit of this approach is that the last step only counts filled cells, as you can see in the image.

So you’re eventual solution might be:

```
rescale <- function(x, x0, xm, n) {
(x - x0)/(xm - x0)*n
}
count_filled_cells <- function(x, y, xmin, xmax, ymin, ymax, ngridx, ngridy) {
ret <- table(ceiling(rescale(x, xmin, xmax, ngridx)), ceiling(rescale(y, ymin, ymax, ngridy)))
sum(ret > 0)
}
count_filled_cells(x, y, 0, 100, 0, 100, 10, 10)
# now apply across range of grid sizes
sizes <- c(2, 5, 10, 20)
answer <- sapply(sizes, function(n) count_filled_cells(x, y, 0, 100, 0, 100, n, n))
```