Vectorizing custome RGB -> Grayscale conversion

Question

As described in the headline I want to make a very specific conversion from RGB to Grayscale.

I have a bunch of images that might look like this:

and i want to convert them to an image like this .

Now you might wonder why I am not just using opencv's inbuilt functions. The reason is that I need to map each color of the RGB image to a specific intensity value in grayscale, which is not all to difficult since I have only six colors.

Red, rgb(255,0,0)           ->  25
Brown, rgb(165,42,42)       ->  120
Light Blue, rgb(0,255,255)  ->  127
Green, rgb(127,255,0)       ->  50
Yellow, rgb(255,255,255)    ->  159
Purple, rgb(128, 0, 128)    ->  90

Now I have already created an array with some objects that contain these mappings and I am simply iterating over my images to assign the new color codes. However this is very slow and i expect to grow a magnificent beard before this is finished for all the images (also I want to know this for learning purpose). This is my super slow running code and the mapping object so far:

colorMapping = [ColorMapping(RGB=[255, 0, 0], Grayscale=25),
 ColorMapping(RGB=[165, 42, 42], Grayscale=120), ... ]

def RGBtoGray(RGBimg, colorMapping):
    RGBimg = cv2.cvtColor(RGBimg, cv2.COLOR_BGR2RGB)
    row = RGBimg.shape[0]
    col = RGBimg.shape[1]
    GRAYimg = np.zeros((row, col))
    for x in range(0,row):
        for y in range(0,col):
            pixel = RGBimg[x,y,:]
            for cm in colorMapping:
                if np.array_equal(pixel, np.array(cm.RGB)):
                    GRAYimg[x,y] = cm.Grayscale    
    return GRAYimg   

I am glad for any suggestions for using built in libraries or improving this codes computation. The color map is read from a json file, which functions as a automation step since I have to do this at least for two batches of images with different encodings.

Answer 1

Method #1

Here's one vectorized based on 1D transformation + np.searchsorted, inspired by this post -

def map_colors(img, colors, vals, invalid_val=0):
    s = 256**np.arange(3)
    img1D = img.reshape(-1,img.shape[2]).dot(s)
    colors1D = colors.reshape(-1,img.shape[2]).dot(s)
    sidx = colors1D.argsort()
    idx0 = np.searchsorted(colors1D, img1D, sorter=sidx)
    idx0[idx0==len(sidx)] = 0
    mapped_idx = sidx[idx0]
    valid = colors1D[mapped_idx] == img1D
    return np.where(valid, vals[mapped_idx], invalid_val).reshape(img.shape[:2])

Sample run -

# Mapping colors array
In [197]: colors
Out[197]: 
array([[255,   0,   0],
       [165,  42,  42],
       [  0, 255, 255],
       [127, 255,   0],
       [255, 255, 255],
       [128,   0, 128]])

# Mapping values array
In [198]: vals
Out[198]: array([ 25, 120, 127,  50, 155,  90])

# Input 3D image array
In [199]: img
Out[199]: 
array([[[255, 255, 255],
        [128,   0, 128],
        [255,   0,   0],
        [127, 255,   0]],

       [[127, 255,   0],
        [127, 255,   0],
        [165,  42,  42],
        [  0,   0,   0]]]) # <= one color absent in mappings

# Output
In [200]: map_colors(img, colors, vals, invalid_val=0)
Out[200]: 
array([[155,  90,  25,  50],
       [ 50,  50, 120,   0]])

We could pre-sort the mappings and hence, get rid of sorting needed around searchsorted and this should further boost performance -

def map_colors_with_sorting(img, colors, vals, invalid_val=0):
    s = 256**np.arange(3)
    img1D = img.reshape(-1,img.shape[2]).dot(s)
    colors1D = colors.reshape(-1,img.shape[2]).dot(s)
    sidx = colors1D.argsort()
    colors1D_sorted = colors1D[sidx]
    vals_sorted = vals[sidx]

    idx0 = np.searchsorted(colors1D_sorted, img1D)
    idx0[idx0==len(sidx)] = 0
    valid = colors1D_sorted[idx0] == img1D
    return np.where(valid, vals_sorted[idx0], invalid_val).reshape(img.shape[:2])

---

Method #2

We can use a mapping array that when indexed by 1D transformed colors would lead us directly to the final "grayscale" image, as shown below -

def map_colors_with_mappingar_solution(img):
    # Edit the custom colors and values here
    colors = np.array([
        [  0,   0, 255],
        [ 42,  42, 165],
        [255, 255,   0],
        [  0, 255, 127],
        [255, 255, 255],
        [128,   0, 128]], dtype=np.uint8) # BGR format
    vals = np.array([25, 120, 127, 50, 155, 90], dtype=np.uint8)      
    return map_colors_with_mappingar(img, colors, vals, 0)            

def map_colors_with_mappingar(img, colors, vals, invalid_val=0):
    s = 256**np.arange(3)
    img1D = img.reshape(-1,img.shape[2]).dot(s)
    colors1D = colors.reshape(-1,img.shape[2]).dot(s)
    
    N = colors1D.max()+1
    mapar = np.empty(N, dtype=np.uint8)
    mapar[colors1D] = vals
    
    mask = np.zeros(N, dtype=bool)
    mask[colors1D] = True
    
    valid = img1D < N
    valid &= mask[img1D]
    
    out = np.full(len(img1D), invalid_val, dtype=np.uint8)
    out[valid] = mapar[img1D[valid]]
    return out.reshape(img.shape[:2])

This should scale well as you increase the number of custom colors.

Let's time it for the given sample image -

# Read in sample image
In [360]: im = cv2.imread('blobs.png')

# @Mark Setchell's solution
In [362]: %timeit remap2(im)
7.45 ms ± 105 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# Method2 from this post
In [363]: %timeit map_colors_with_mappingar_solution(im)
6.76 ms ± 46.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Further perf. boost

Going one step further, we could do the 1D reduction in a more performant way and hence achieve further perf. boost, like so -

# https://stackoverflow.com/a/57236217/ @tstanisl
def scalarize(x):
    # compute x[...,2]*65536+x[...,1]*256+x[...,0] in efficient way
    y = x[...,2].astype('u4')
    y <<= 8
    y +=x[...,1]
    y <<= 8
    y += x[...,0]
    return y

def map_colors_with_mappingar(img, colors, vals, invalid_val=0):    
    img1D = scalarize(img)
    colors1D = scalarize(colors)
    
    N = colors1D.max()+1
    mapar = np.empty(N, dtype=np.uint8)
    mapar[colors1D] = vals
    
    mask = np.zeros(N, dtype=bool)
    mask[colors1D] = True
    
    valid = img1D < N
    valid &= mask[img1D]
    
    out = np.full(img1D.shape, invalid_val, dtype=np.uint8)
    out[valid] = mapar[img1D[valid]]
    return out

# On given sample image
In [10]: %timeit map_colors_with_mappingar_solution(im)
5.45 ms ± 143 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Answer 2

This is probably as easy as anything else. It does make 6 passes over your image, so some clever Numpy folk may know a better way, but it will be a lot faster than your loops.

#!/usr/bin/env python3

import numpy as np
import cv2

# Load image
im = cv2.imread('blobs.png')

# Make output image
res = np.zeros_like(im[:,:,0])

res[np.all(im == (0, 0, 255),   axis=-1)] = 25
res[np.all(im == (42,42,165),   axis=-1)] = 120
res[np.all(im == (255,255,0),   axis=-1)] = 127
res[np.all(im == (0,255,127),   axis=-1)] = 50
res[np.all(im == (255,255,255), axis=-1)] = 159
res[np.all(im == (128,0,128),   axis=-1)] = 90

# Write image of just palette indices
cv2.imwrite('indices.png', res)

---

You can make it run in 5ms versus 30ms by converting each RGB triplet into a single 24-bit number as inspired by this answer.

#!/usr/bin/env python3

import numpy as np
import cv2
def remap2(im):
    # Make output image
    res = np.zeros_like(im[:,:,0])

    # Make a single 24-bit number for each pixel
    r = np.dot(im.astype(np.uint32),[1,256,65536]) 

    c0 =   0 +   0*256 + 255*65536
    c1 =  42 +  42*256 + 165*65536
    c2 = 255 + 255*256 +   0*65536
    c3 =   0 + 255*256 + 127*65536
    c4 = 255 + 255*256 + 255*65536
    c5 = 128 +   0*256 + 128*65536

    res[r == c0] = 25
    res[r == c1] = 120
    res[r == c2] = 127
    res[r == c3] = 50
    res[r == c4] = 159
    res[r == c5] = 90
    return res

# Load image
im = cv2.imread('blobs.png')
res = remap2(im)
cv2.imwrite('result.png',res)