Tiago,
No seu exemplo, o volume seria:
x <- y <- z <- 0:100
r <- 30
xyz <- expand.grid(x,y,z)
dados <- cbind(xyz,as.numeric(sqrt(rowSums((xyz-50.5)^2))<=r))
colnames(dados) <- c("x","y","z","var")
mean(dados$var)*100^3 # Aproximação
4*pi*30^3/3 # Exato
A idéia aqui é que a proporção de pontos que caem dentro da esfera é aproximadamente a razão entre o volume da esfera e
do cubo. Essa aproximação pode ser melhorada aumentando-se o número de pontos da malha.
Não tinha visto na sua postagem original, mas se vc já tem
uma variável indicadora com essa informação, o mesmo princípio se aplica.
Att.,
Rubem