Scikit-Image, Numpy, and Selecting Colors (python)

Question

I'm trying to mask colored features from a photograph so that I can do some other processing on them. I've played with a few packages (scikit-image, mahotas, and openCV) and have settled on Scikit-image because it plays well with scikit-learn.

I'd like to accomplish two things, eventually: (1) get colored features so that I can compute things like length, and (2) remove colored features from the image while retaining the white spidering-looking veins (the lattice-looking stuff around the margin of the image will need to go, but that can be a later endeavour).

I'm new to image processing in Python, and I think I've mislead myself by using strange color spaces. Perhaps it would have been better to define an upper and lower boundary for each colour, and mask the image array? I'm definitely open to any flexible solution to the problem. Eventually I'll be applying this to hundreds of photos..

Original Photo

Code

from skimage import io, img_as_float, color, exposure

img = img_as_float(io.imread('./images/testimage2.JPG'))

# Isolate paint marks

# Put image into LAB colour space
image_lab = color.rgb2lab(img)
img = exposure.rescale_intensity(img)

# Colours of interest
color_array = np.array([
        [[[255, 255, 0.]]],    # Yellow stuff
        [[[255, 190, 200.]]],    # Pink stuff   
        [[[255, 165, 0.]]],    # Orange stuff
        [[[255, 0, 0.]]],    # Red stuff        
        ])

# Loop through the color array and pick out the colored features
for i in range(0,color_array.ndim):
    # Compute distance between the color of interest and the actual image
    # http://scikit-image.org/docs/dev/api/skimage.color.html#skimage.color.deltaE_cmc
    # "The usual values are kL=2, kC=1 for “acceptability”"
    distance_color = color.deltaE_ciede2000(color_array[i], image_lab, kL=2, kC=1, kH=0.5)    

    # Normalise distance
    distance_color = exposure.rescale_intensity(distance_color)

    # Mask image
    image_filtered = img.copy()
    image_filtered[distance_color > 0.5] = 0

    # Plot it up
    print ("Filtered to: ", color_array[i])
    f, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(20, 10))
    ax0.imshow(img)
    ax1.imshow(distance_color, cmap='gray')
    ax2.imshow(image_filtered)
    plt.show()

Results

Answer 1

You can look at your picture in a different color space, e.g. the HSV space. In this space, each pixel has 3 components: hue, saturation and value. Hue defines the color's tone, saturation the color's strength and V the brightness of the pixel. For your particular problem, you can search for pixels that have e.g. a decent saturation (i.e. are reasonably colored and not just gray). The following code snippet uses OpenCV, but I believe you will be able to translate into scikit-image:

I = cv2.imread("/tmp/XApdz.jpg")  # Load the image

hsv = cv2.cvtColor(I, cv2.COLOR_BGR2HSV);  # Convert to HSV space
h, s, v = cv2.split(hsv);  # extract the single components

plt.figure(1, figsize=(10, 10))
plt.subplot(3,2,1); plt.imshow(h, cmap='gray'); plt.title('Hue')
plt.subplot(3,2,2); plt.imshow(s, cmap='gray'); plt.title('Saturation')
plt.subplot(3,2,3); plt.imshow(v, cmap='gray'); plt.title('Value')

mask = (s > 80).astype(np.uint8);  # Thresholding in the Saturation-channel
plt.subplot(3,2,4); plt.imshow(mask); plt.title('mask')

disk = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(3,3))  # remove very small regions
opened = cv2.morphologyEx(mask, cv2.MORPH_OPEN, disk);
plt.subplot(3,2,5); plt.imshow(opened); plt.title('Opened mask')

I2 = I.copy()
I2[opened.astype(bool), :] = 0; # Set the pixels to zero, where 
plt.subplot(3,2,6); plt.imshow(I2) 

Answer 2

Thanks to Maximilian Matthé, I have an answer. Below is his OpenCV code translated to Scikit-image. Note that I changed some parameters (e.g., reduce the size of the morphological disk) to replicate his results.

from skimage import color
from skimage.morphology import disk, opening, dilation

img_hsv = color.rgb2hsv(img) # Image into HSV colorspace
h = img_hsv[:,:,0] # Hue
s = img_hsv[:,:,1] # Saturation
v = img_hsv[:,:,2] # Value aka Lightness

plt.figure(1, figsize=(15, 15))

plt.subplot(4,2,1); plt.imshow(h, cmap='gray'); plt.title('Hue')
plt.subplot(4,2,2); plt.imshow(s, cmap='gray'); plt.title('Saturation')
plt.subplot(4,2,3); plt.imshow(v, cmap='gray'); plt.title('Value')
plt.tight_layout()

mask = (s > 0.35).astype(np.uint8);  # Thresholding in the Saturation-channel
plt.subplot(4,2,4); plt.imshow(mask); plt.title('mask')

disk_elem = disk(1) # Remove small regions
opened = opening(mask, selem=disk_elem)
plt.subplot(4,2,5); plt.imshow(opened); plt.title('Opened mask')

square_elem = square(2) # rejoin colored features
dilated = dilation(opened, selem=square_elem)
plt.subplot(4,2,6); plt.imshow(dilated); plt.title('Opened mask')

img2 = img.copy()
img2[dilated.astype(bool), :] = 0; # Set the pixels to zero, where 
plt.subplot(4,2,7); plt.imshow(img2); plt.title('Final Image')