Array slice maximum that depends on the index of the previous axis

Question

So I have a large 2D array, coming from a tiff image, in which I want to calculate the center of mass. To do that, I am using the indices of the image (as coordinates) and the average function:

from PIL import Image
from numpy import *

Im    = Image.open("32bit_grayscale.tif")
imArr = array(Im, dtype='float32')
indx  = indices(imArr.shape)

cenMassX = average(indx[0,:,:],weights=imArr[:,:])
cenMassY = average(indx[1,:,:],weights=imArr[:,:])

In some other similar images, there are two separate centers of mass that I want to calculate. Both regions of calculation are separated by a straight oblicuous line, of which I have its equation.

I would like to use the average method again, as it is very efficient, but I would need to set the maximum of the slice of the second axis of the indx array to a function of the current first axis value. If the line was something like y=slope*x+interY, I would need something like this:

cenMassX_A = average(indx[0,:,:slope*row+interY],weights=imArr[:,:slope*row+interY])
cenMassY_A = average(indx[1,:,:slope*row+interY],weights=imArr[:,:slope*row+interY])
cenMassX_B = average(indx[0,:,slope*row+interY:],weights=imArr[:,slope*row+interY:])
cenMassY_B = average(indx[1,:,slope*row+interY:],weights=imArr[:,slope*row+interY:])

Where row represents the current value of the first axis index (the "x" axis). Disregard the fact that I could go off the array limit depending on the equation.

I can do this using for loops, but it is very inefficient (by a factor 20) and not very "pythonic":

cenMassX_A = 0
cenMassY_A = 0
cumSum     = 0
for row in range(0,imArr.shape[0]):
    for col in range(0,int(round(slope*row+interY))):
        cenMassX_A += row*imArr[row,col]
        cenMassY_A += col*imArr[row,col]
        cumSum     += imArr[row,col]
cenMassX_A /= cumSum
cenMassY_A /= cumSum

cenMassX_B = 0
cenMassY_B = 0
cumSum     = 0
for row in range(0,imArr.shape[0]):
    for col in range(int(round(slope*row+interY)),imArr.shape[1]):
        cenMassX_B += row*imArr[row,col]
        cenMassY_B += col*imArr[row,col]
        cumSum     += imArr[row,col]
cenMassX_B /= cumSum
cenMassY_B /= cumSum

So, is there a way to do this, or am I stuck with the for loops? I have been reading about masks and rolling windows but still cannot figure out a solution. Thanks in advance!

Answer 1

What happens if you set imArr[i,j]=0 for all points on one side or the other of your line? This is the simplest masking approach.

I = indx[0,...]*slope + indx[1,...]>=M
imArr1 = imArr.copy()
imArr1[I]=0
print np.average(indx[0,...],weights=imArr1)
print np.average(indx[1,...],weights=imArr1)

imArr1 = imArr.copy()
imArr1[~I]=0

print np.average(indx[0,...],weights=imArr1)
print np.average(indx[1,...],weights=imArr1)

This works fine if I take a simple 'image', and concatenate it to itself (horizontally or diagonally).

Answer 2

Why don't you use the code you wrote:

cenMassX_A = average(indx[0,:,:slope*row+interY],weights=inArr[:,:slope*row+interY])
cenMassY_A = average(indx[1,:,:slope*row+interY],weights=inArr[:,:slope*row+interY])

and apply it row-wise?

Keep the for-loops over the rows, but replace your col-loops with this code (changing the row-index from all rows to just the one in the for loop). And of course you still need to add the results from the rows together.