What is a fast(er) way to get the center points of objects represented in a 2D numpy array?-CodePudding

I have an image mask stored as a 2D numpy array where the values indicate the presence of objects that have been segmented in the image (0 = no object, 1..n = object 1 through n). I want to get a single coordinate for each object representing the center of the object. It doesn't have to be a perfectly accurate centroid or center of gravity. I'm just taking the mean of the x and y indices of all cells in the array that contain each object. I'm wondering if there's a faster way to do this than my current method:

for obj in np.unique(mask):
    if obj == 0:
        continue
    x, y = np.mean(np.where(mask == obj), axis=1)

Here is a reproducible example:

import numpy as np
mask = np.array([
    [0,0,0,0,0,2,0,0,0,0],
    [0,1,1,0,2,2,2,0,0,0],
    [0,0,1,0,2,2,2,0,0,0],
    [0,0,0,0,0,0,0,0,0,0],
    [0,3,3,3,0,0,4,0,0,0],
    [0,0,0,0,0,4,4,4,0,0],
    [0,0,0,0,0,0,4,0,0,0],
])

points = []
for obj in np.unique(mask):
    if obj == 0:
        continue
    points.append(np.mean(np.where(mask == obj), axis=1))
print(points)

This outputs:

[array([1.33333333, 1.66666667]),
 array([1.28571429, 5.        ]),
 array([4., 2.]),
 array([5., 6.])]

CodePudding user response：

I came up with another way to do it that seems to be about 3x faster:

import numpy as np
mask = np.array([
    [0,0,0,0,0,2,0,0,0,0],
    [0,1,1,0,2,2,2,0,0,0],
    [0,0,1,0,2,2,2,0,0,0],
    [0,0,0,0,0,0,0,0,0,0],
    [0,3,3,3,0,0,4,0,0,0],
    [0,0,0,0,0,4,4,4,0,0],
    [0,0,0,0,0,0,4,0,0,0],
])

flat = mask.flatten()
split = np.unique(np.sort(flat), return_index=True)[1]
points = []
for inds in np.split(flat.argsort(), split)[2:]:
    points.append(np.array(np.unravel_index(inds, mask.shape)).mean(axis=1))
print(points)

I wonder if the for loop can be replaced with a numpy operation which would likely be even faster.

CodePudding user response：

You can copy this answer (give them an upvote too if this answer works for you) and use sparse matrices instead of np arrays. However, this only proves to be quicker for large arrays, with increasing speed boosts the larger your array is:

import numpy as np, time
from scipy.sparse import csr_matrix

def compute_M(data):
    cols = np.arange(data.size)
    return csr_matrix((cols, (np.ravel(data), cols)),
                      shape=(data.max()   1, data.size))

def get_indices_sparse(data,M):
    #M = compute_M(data)
    return [np.mean(np.unravel_index(row.data, data.shape),1) for R,row in enumerate(M) if R>0]

def gen_random_mask(C, n, m):
    mask = np.zeros([n,m],int)
    for i in range(C):
        x = np.random.randint(n)
        y = np.random.randint(m)
        mask[x:x np.random.randint(n-x),y:y np.random.randint(m-y)] = i
    return mask

N = 100
C = 4
for S in [10,100,1000,10000]:
    mask = gen_random_mask(C, S, S)
    print('Time for size {:d}x{:d}:'.format(S,S))
    s = time.time()
    for _ in range(N):
        points = []
        for obj in np.unique(mask):
            if obj == 0:
                continue
            points.append(np.mean(np.where(mask == obj), axis=1))
    points_np = np.array(points)
    print('NP: {:f}'.format((time.time() - s)/N))
    mask_s = compute_M(mask)
    s = time.time()
    for _ in range(100):
        points = get_indices_sparse(mask,mask_s)
    print('Sparse: {:f}'.format((time.time() - s)/N))
    np.testing.assert_equal(points,points_np)

Which results in the timings of:

Time for size 10x10:
NP: 0.000066
Sparse: 0.000226
Time for size 100x100:
NP: 0.000207
Sparse: 0.000253
Time for size 1000x1000:
NP: 0.018662
Sparse: 0.004472
Time for size 10000x10000:
NP: 2.545973
Sparse: 0.501061

CodePudding user response：

The problem likely comes from np.where(mask == obj) which iterates on the whole mask array over and over. This is a problem when there are a lot of objects. You can solve this problem efficiently using a group-by strategy. However, Numpy do not yet provide such an operation. You can implement that using a sort followed by a split. But a sort is generally not efficient. An alternative method is to ask Numpy to return the index in the unique call so that you can then accumulate the value regarding the object (like a reduce-by-key where the reduction operator is an addition and the key are object integers). The mean can be obtained using a simple division in the end.

objects, inverts, counts = np.unique(mask, return_counts=True, return_inverse=True)

# Reduction by object
x = np.full(len(objects), 0.0)
y = np.full(len(objects), 0.0)
xPos = np.repeat(np.arange(mask.shape[0]), mask.shape[1])
yPos = np.tile(np.arange(mask.shape[1]), reps=mask.shape[0])
np.add.at(x, inverts, xPos)
np.add.at(y, inverts, yPos)

# Compute the final mean from the sum
x /= counts
y /= counts

# Discard the first item (when obj == 0)
x = x[1:]
y = y[1:]

If you need something faster, you could use Numba and perform the reduction manually (and possibly in parallel).

EDIT: if you really need a list in output, you can use points = list(np.stack([x, y]).T) but this is rather slow to use lists instead of Numpy arrays (and not memory efficient either).