Felix
Felix

Reputation: 739

python: calculate center of mass

I have a data set with 4 columns: x,y,z, and value, let's say:

x  y  z  value
0  0  0  0
0  1  0  0
0  2  0  0
1  0  0  0
1  1  0  1
1  2  0  1
2  0  0  0
2  1  0  0
2  2  0  0

I would like to calculate the center of mass CM = (x_m,y_m,z_m) of all values. In the present example, I would like to see (1,1.5,0) as output.

I thought this must be a trivial problem, but I can't find a solution to it in the internet. scipy.ndimage.measurements.center_of_mass seems to be the right thing, but unfortunately, the function always returns two values (instead of 3). In addition, I can't find any documentation on how to set up an ndimage from an array: Would I use a numpy array N of shape (9,4)? Would then N[:,0] be the x-coordinate?

Any help is highly appreciated.

Upvotes: 11

Views: 50138

Answers (4)

gb96
gb96

Reputation: 1694

Why did ndimage.measurements.center_of_mass not give the expected result?

The key is in how the input data masses was represented by an array of 4-tuples (x, y, z, value)

# x   y   z   value
[[0,  0,  0,  0],
 [0,  1,  0,  0],
 [0,  2,  0,  0],
 [1,  0,  0,  0],
 [1,  1,  0,  1],
 [1,  2,  0,  1],
 [2,  0,  0,  0],
 [2,  1,  0,  0],
 [2,  2,  0,  0]]

The array masses here represents the 3-D position and weights of each mass. Note however that this python array structure is only a 2-D array. It's shape is (9, 4).

The input you need to pass to ndimage to get the expected result is a 3-D array containing zeros everywhere and the weight of each mass at the appropriate coordinates within the array, like this:

from scipy import ndimage
import numpy

masses = numpy.zeros((3, 3, 1))
#      x  y  z    value
masses[1, 1, 0] = 1
masses[1, 2, 0] = 1

CM = ndimage.measurements.center_of_mass(masses)
#  x    y    z
# (1.0, 1.5, 0.0)

Which is exactly the expected output.

Note the limitation of this solution (and the ndimage library) is it requires non-negative integer coordinates. Also will not be efficient for large and/or sparse volumes because each "pixel" of the ndimage needs to be instantiated in memory.

Upvotes: 1

Bachbold
Bachbold

Reputation: 434

Another option is to use the scipy center of mass:

from scipy import ndimage
import numpy

masses = numpy.array([[0,  0,  0,  0],
[0,  1,  0,  0],
[0,  2,  0,  0],
[1,  0,  0,  0],
[1,  1,  0,  1],
[1,  2,  0,  1],
[2,  0,  0,  0],
[2,  1,  0,  0],
[2,  2,  0,  0]])

ndimage.measurements.center_of_mass(masses)

Upvotes: 5

Bi Rico
Bi Rico

Reputation: 25813

How about:

#                   x      y     z  value
table = np.array([[ 5. ,  1.3,  8.3,  9. ],
                  [ 6. ,  6.7,  1.6,  5.9],
                  [ 9.1,  0.2,  6.2,  3.7],
                  [ 2.2,  2. ,  6.7,  4.6],
                  [ 3.4,  5.6,  8.4,  7.3],
                  [ 4.8,  5.9,  5.7,  5.8],
                  [ 3.7,  1.1,  8.2,  2.2],
                  [ 0.3,  0.7,  7.3,  4.6],
                  [ 8.1,  1.9,  7. ,  5.3],
                  [ 9.1,  8.2,  3.3,  5.3]])

def com(xyz, mass):
    mass = mass.reshape((-1, 1))
    return (xyz * mass).mean(0)

print(com(table[:, :3], table[:, 3]))

Upvotes: 2

Aleksander Lidtke
Aleksander Lidtke

Reputation: 2926

The simplest way I can think of is this: just find an average of the coordinates of mass components weighted by each component's contribution.

import numpy
masses = numpy.array([[0,  0,  0,  0],
[0,  1,  0,  0],
[0,  2,  0,  0],
[1,  0,  0,  0],
[1,  1,  0,  1],
[1,  2,  0,  1],
[2,  0,  0,  0],
[2,  1,  0,  0],
[2,  2,  0,  0]])

nonZeroMasses = masses[numpy.nonzero(masses[:,3])] # Not really necessary, can just use masses because 0 mass used as weight will work just fine.

CM = numpy.average(nonZeroMasses[:,:3], axis=0, weights=nonZeroMasses[:,3])

Upvotes: 16

Related Questions