[FIXED] Whats the fastest method for iou in numpy?

Issue

I'm currently trying to do an iou analysis for a 3d image in different time points. These images contain around 1500 objects (cells) with an unique id in an around 2000x2000x2000 image.

I found logical_and and logical_or functions of numpy which take one variable at a time from each image so I made a very basic double for loop, to feed every combination of every value into the iou analysis. The code looks like this:

for i in [x for x in image_1_ids if x != 0]:
    for j in [y for y in image_2_ids if y != 0]:
    
        intersection = np.logical_and(image_1 == i, image_2 == j)
        union = np.logical_or(image_1 == i, image_2 == j)

        iou = np.sum(intersection) / np.sum(union)

        df.at[i, j] = iou

This code takes forever to run due to the many variable feed one at a time. Which makes it a combination of basically 1500x1500. Is there a more efficient way in doing this in numpy?

Solution

A faster algorithm is possible by:

• Exploiting (more) the fact that object ids can be used as indices
• Calculating the union_ij basically as count_i + count_j - intersection_ij

The algorithm basically works like this:

n_1 = np.max(image_1) + 1
n_2 = np.max(image_2) + 1

counts_1 = np.zeros(n_1, int)
counts_2 = np.zeros(n_2, int)
intersection = np.zeros((n_1, n_2), int)
for i, j in zip(image_1.flat, image_2.flat):
    counts_1[i] += 1
    counts_2[j] += 1
    intersection[i, j] += 1

union = counts_1[:, np.newaxis] + counts_2 - intersection

iou = intersection / union

When we understand how that works it's time to look at some improvements:

• Get the counts without a for-loop
• Avoid divide-by-zero (in case object ids are not consecutive)
• Correct the result where the images contain zero

In code:

counts_1 = np.bincount(image_1.ravel('K'))
counts_2 = np.bincount(image_2.ravel('K'))

shape = (len(counts_1), len(counts_2))

linear_ids = np.ravel_multi_index((image_1, image_2), shape)
intersection = np.bincount(linear_ids.ravel(), minlength=np.product(shape))
intersection = intersection.reshape(shape)

union = counts_1[:, np.newaxis] + counts_2 - intersection

iou = np.zeros(shape, float)
np.divide(intersection, union, out=iou, where=(union != 0))
iou[0, :] = iou[:, 0] = 0

This doesn't take into account image_1_ids and image_2_ids yet. It's a bit unclear what needs to happen here, but with the result from above you can probably do something like this:

iou = iou[image_1_ids[:, np.newaxis], image_2_ids]

Answered By - user7138814