Using numpy to modify RGB values on entire image, instead of looping pixels

Question

I have a function that loops the RGB pixels of an opencv image and applies some changes to each pixel. But this approach is slow and I think it will be much faster if I can use numpy to do the changes all at once.

How can I make these filters in Numpy instead of in a loop?

def apply_image_filter(image):
    image = cv2.normalize(image, image, 0, 255, cv2.NORM_MINMAX)

    h = image.shape[0]
    w = image.shape[1]

    # loop over the image, pixel by pixel
    for y in range(0, h):
        for x in range(0, w):
            if np.all(image[y, x] == 0): # Ignore black pixels
                continue

            mean = np.mean(image[y, x])
            maxdex = np.argmax(image[y, x])
            maxval = image[y, x][maxdex]

            # Equal RGB go to white or black
            if np.all(image[y, x] == mean):
                image[y, x] = [255, 255, 255] if mean > 100 else [0, 0, 0]
                continue

            # High color equals go to white
            if np.all((mean * 0.85) < image[y, x]) and np.all((mean * 1.15) > image[y, x]):
                image[y, x] = [255, 255, 255] if mean >= 90 else [0, 0, 0]
                continue

            if maxval > mean * 1.5:
                image[y, x] = [0, 0, 0]
                image[y, x][maxdex] = 255
                continue

            median = np.median(image[y, x])

            if maxval > median * 1.5:
                image[y, x] = [0, 0, 0]
                image[y, x][maxdex] = 255
                continue

    return image

More specifically how do I move the conditional pixel updates outside of the the loop.

For example this line:

if np.all((mean * 0.85) < image[y, x]) and np.all((mean * 1.15) > image[y, x]):
    image[y, x] = [255, 255, 255] if mean >= 90 else [0, 0, 0]

It checks that the R, G, and B values in the pixel are all within 15% of the mean value of the pixel and if so then sets the pixel to either black or white depending on the mean value.

How can I move this conditional update so that it does all pixels at once in one Numpy command, instead of on a per pixel basis?

Updated from comments:

I have moved the mean, argmax and median out of the loop so they are first applied to the entire array. Now the execution speed is about 60% faster than it was. However I'm still doing conditional statements for each pixel.

def apply_image_filter(image):
    image = cv2.normalize(image, image, 0, 255, cv2.NORM_MINMAX)

    h = image.shape[0]
    w = image.shape[1]

    # These Numpy operations are now outside the loop
    means = np.mean(image, axis=2)
    maxdexes = np.argmax(image, axis=2)
    medians = np.median(image, axis=2)

    # loop over the image, pixel by pixel
    for y in range(0, h):
        for x in range(0, w):
            if means[y, x] == 0: # Ignore black pixels
                continue
            
            # assign to local variables for readability
            mean = means[y, x]
            maxdex = maxdexes[y, x]
            maxval = image[y, x][maxdex]

            # Equal RGB go to white or black
            if np.all(image[y, x] == mean):
                image[y, x] = [255, 255, 255] if mean > 100 else [0, 0, 0]
                continue

            # High color equals go to white
            if np.all((mean * 0.85) < image[y, x]) and np.all((mean * 1.15) > image[y, x]):
                image[y, x] = [255, 255, 255] if mean >= 90 else [0, 0, 0]
                continue

            if maxval > mean * 1.5:
                image[y, x] = [0, 0, 0]
                image[y, x][maxdex] = 255
                continue

            if maxval > medians[y, x] * 1.5:
                image[y, x] = [0, 0, 0]
                image[y, x][maxdex] = 255
                continue

    return image

Update #2 (my solution)

After reading the comments and answer I have spent time trying to figure out how to "vectorize" my code, and I believe that I am managed to move every one of my conditional pixel statements into a numpy vector operation, here is my updated code, that does the same thing as the origional code, but runs about 100 times faster.

I feel like it loses some readability, but for 100x speed up i'll just have to add good comments.

def apply_image_filter(image):
    image = cv2.normalize(image, image, 0, 255, cv2.NORM_MINMAX)

    means = np.mean(image, axis=2)
    maxdexes = np.argmax(image, axis=2)

    maxes = np.amax(image, axis=2)
    mins = np.amin(image, axis=2)

    image[np.logical_and(means == maxes, means > 100)] = 255
    image[np.logical_and(means == maxes, means <= 100)] = 0

    up_means = means * 1.15
    down_means = means * 0.85

    image[np.logical_and(means >= 90, np.logical_and(up_means > maxes, down_means < mins))] = 255
    image[np.logical_and(means < 90, np.logical_and(up_means > maxes, down_means < mins))] = 0

    up_means = means * 1.5
    new_means = np.mean(image, axis=2)

    higher_than_mean_logic = np.logical_and(maxes > up_means, new_means < 255)
    image[higher_than_mean_logic] = 0
    image[higher_than_mean_logic, maxdexes[higher_than_mean_logic]] = 255

    new_maxes = np.amax(image, axis=2)
    medians = np.median(image, axis=2)
    medians_logic = np.logical_and(new_maxes > medians * 1.5, new_maxes < 255)
    image[medians_logic] = 0
    image[medians_logic, maxdexes[medians_logic]] = 255

    return image

Answer 1

Python loops like the one you do are slow because of the CPython interpreter but also because scalar access to Numpy are very slow since Numpy is not optimized for that (and even if it would be, CPython prevent the code to be fast due to the interpreter itself and the way it is designed). The general solution is to vectorize your code using high-level Numpy calls instead of scalar access in loops. However, while this is possible in your case, it would results in many temporary arrays to be created which is not efficient (although it should be at least one order of magnitude faster than your current code). A faster and simpler way is simply to use Numba for such task when there is no simple Numpy functions that can be used.

Here is the resulting (untested) Numba code:

@nb.njit('void(uint8[:,:,::1], float64[:,::1], int64[:,::1], float64[:,::1])')
def computeLoop(image, means, maxdexes, medians):
    assert image.shape[0] == 3

    # loop over the image, pixel by pixel
    for y in range(0, h):
        for x in range(0, w):
            if means[y, x] == 0: # Ignore black pixels
                continue
            
            # assign to local variables for readability
            mean = means[y, x]
            maxdex = maxdexes[y, x]
            maxval = image[y, x, maxdex]

            # Equal RGB go to white or black
            if image[y, x, 0] == mean and image[y, x, 1] == mean and image[y, x, 2] == mean:
                value = 255 if mean > 100 else 0
                image[y, x, :] = value
                continue

            mini, maxi = mean * 0.85, mean * 1.15

            # High color equals go to white
            if mini < image[y, x, 0] < maxi and mini < image[y, x, 1] < maxi and mini < image[y, x, 2] < maxi:
                value = 255 if mean >= 90 else 0
                image[y, x, :] = value
                continue

            if maxval > mean * 1.5 or maxval > medians[y, x] * 1.5:
                image[y, x, :] = 0
                image[y, x, maxdex] = 255
                continue

def apply_image_filter(image):
    image = cv2.normalize(image, image, 0, 255, cv2.NORM_MINMAX)

    h = image.shape[0]
    w = image.shape[1]

    # These Numpy operations are now outside the loop
    means = np.mean(image, axis=2)
    maxdexes = np.argmax(image, axis=2)
    medians = np.median(image, axis=2)

    computeLoop(image, means, maxdexes, medians)

    return image

This code should certainly be several order of magnitude faster than the initial one. Moreover, note that the code can be parallelized by adding the njit flag parallel=True and using np.prange instead of range in the encompassing y-indexed loop. Note that it is your responsibility to check whether the code can be parallelized or not (otherwise the behavior is undefined: results could be wrong and the program can crash).