center of mass in each cell of N-dimensional rectilinear grid

Question

I have an n-Dimensional rectilinear grid, such that the edges in each dimension i are given by x_i = {x[i, 0], x[i,1], ..., x[i, Ni-1], x[i, Ni]}, with N_i edges in that dimension. I then have some density y, with shape (N0, N1, ... Ni, ... Nn-1) defined at each grid vertex. We can assume that the density varies smoothly, and the density between vertices (i.e. within a cell) can be calculated by smoothly (linearly) interpolating between vertices. How do I find the center of mass in each dimension, for each cell/bin?

Note that the answer is not just the y-weighted average of the edges of each bin. For example, consider a 2D grid with x-coordinates [0.0, 1.0], and a density [[0.0, 0.0], [1.0, 2.0]]. The y-weighted x position of vertices is 1.0, but clearly the center of mass needs to be somewhere mid-way between the edges, not up against the edge:

0.0--------2.0  -y2
 |          |
 |        * |
 |          |
 |          |
 |          |
0.0--------1.0  -y1

 |        |
 x1=0.0   x2=1.0

where the * approximates the center of mass.

Answer 1

Let us assume that density is distributed over cell using bilinear interpolation with normalized cell coordinates 0..1 and vertices values f[ij]. Every real cell coordinates might be normalized to the range above using linear transformation (subtract left bottom coordinate, divide by size).

To get mass center coordinates, we must apply the next formula

Practically: for 2D case we calculate denominator ("mass") as definite double integral of density(x,y) over x=0..1 and y=0..1 ranges and two nominators (vector components) as integrals of density(x,y)*x and density(x,y)*y. (3D (and nD) cases in vector form look similar, but formulae become more complex.)

Evaluation of these integrals using Maple gives the next result:

which might be easily implemented in Python as code below:

def mass_center_bilinear(f00, f10, f01, f11):
    mass = 1/4*f00+1/4*f10+1/4*f01+1/4*f11
    if mass == 0:
        return None
    x_int = 1/12*f00+1/6*f10+1/12*f01+1/6*f11
    y_int = 1/12*f00+1/12*f10+1/6*f01+1/6*f11
    return (x_int/mass, y_int/mass)

print(mass_center_bilinear(1,1,1,1))
print(mass_center_bilinear(0,0,0,1))
print(mass_center_bilinear(0,0,1,1))
print(mass_center_bilinear(0,1,2,3))
print(mass_center_bilinear(0,1,0,2))

(0.5, 0.5)
(0.6666666666666666, 0.6666666666666666)
(0.5, 0.6666666666666666)
(0.5555555555555555, 0.611111111111111)
(0.6666666666666666, 0.5555555555555555)  #mass center for your example