colorcovid.py

# -*- coding: utf-8 -*-
"""ColorCovid.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1fAQep4O657KhPFrd3wuZAl9krpHxGIam

# Necessary imports
"""

# Commented out IPython magic to ensure Python compatibility.
import numpy as np
# %matplotlib inline
import matplotlib.pyplot as plt
import cv2
import easygui as gui
import sys
import tkinter as tk
import skimage
from skimage.morphology import watershed
from skimage.feature import peak_local_max
from skimage.color import rgb2hsv
import scipy
from scipy import ndimage
from PIL import Image
import imutils
import math
from win32api import GetSystemMetrics
import csv
import time

"""# Process images"""

img = Image.open('tests5.tif')
img = np.array(img.convert("RGB"))
img = img[...,:3]
hsv_img = rgb2hsv(img)*255
hue_img = hsv_img[:, :, 0]
sat_img = hsv_img[:, :, 1]
value_img = hsv_img[:, :, 2]

fig, ax = plt.subplots(ncols=4, nrows=2, figsize=(15,8))

ax[0,0].imshow(img)
ax[0,0].set_title("RGB image")
ax[0,0].axis('off')

ax[0,1].imshow(hue_img, cmap='bwr')
ax[0,1].set_title("Hue channel")
ax[0,1].axis('off')
ax[1,1].hist(hue_img.ravel(), bins=32, range=[0, 256])
ax[1,1].set_xlim(0, 256);

ax[0,2].imshow(sat_img, cmap='bwr')
ax[0,2].set_title("Sat image")
ax[0,2].axis('off')
ax[1,2].hist(sat_img.ravel(), bins=32, range=[0, 256])
ax[1,2].set_xlim(0, 256);

ax[0,3].imshow(value_img, cmap='bwr')
ax[0,3].set_title("Value channel")
ax[0,3].axis('off')
ax[1,3].hist(value_img.ravel(), bins=32, range=[0, 256])
ax[1,3].set_xlim(0, 256);


plt.tight_layout()

hsv_img = rgb2hsv(img)*255
hue_img = hsv_img[:, :, 0]
sat_img = hsv_img[:, :, 1]
val_img = hsv_img[:, :, 2]

thr_hue = 50
hue_img_idx1 = hue_img < thr_hue
hue_img_idx2 = hue_img >= thr_hue
hue_img[hue_img_idx1] = 0
hue_img[hue_img_idx2] = 255

thr_sat = 250
sat_img_idx1 = sat_img < thr_sat
sat_img_idx2 = sat_img >= thr_sat
sat_img[sat_img_idx1] = 0
sat_img[sat_img_idx2] = 255

thr_val = 100
val_img_idx1 = val_img < thr_val
val_img_idx2 = val_img >= thr_val
val_img[val_img_idx1] = 0
val_img[val_img_idx2] = 255

fig, ax = plt.subplots(ncols=4, nrows=2, figsize=(15,8))

ax[0,0].imshow(img)
ax[0,0].set_title("RGB image")
ax[0,0].axis('off')

ax[0,1].imshow(hue_img, cmap='bwr')
ax[0,1].set_title("Hue channel")
ax[0,1].axis('off')
ax[1,1].hist(hue_img.ravel(), bins=32, range=[0, 256])
ax[1,1].set_xlim(0, 256);

ax[0,2].imshow(sat_img, cmap='bwr')
ax[0,2].set_title("Sat image")
ax[0,2].axis('off')
ax[1,2].hist(sat_img.ravel(), bins=32, range=[0, 256])
ax[1,2].set_xlim(0, 256);

ax[0,3].imshow(val_img, cmap='bwr')
ax[0,3].set_title("Value channel")
ax[0,3].axis('off')
ax[1,3].hist(val_img.ravel(), bins=32, range=[0, 256])
ax[1,3].set_xlim(0, 256);



if circles is not None:
    # convert the (x, y) coordinates and radius of the circles to integers
    circles = np.round(circles[0, :]).astype("int")
    # loop over the (x, y) coordinates and radius of the circles
    for (x, y, r) in circles:
        # draw the circle in the output image, then draw a rectangle
        # corresponding to the center of the circle
        cv2.circle(output, (x, y), r, (0, 255, 0), 4)
        cv2.rectangle(output, (x - 5, y - 5), (x + 5, y + 5), (0, 128, 255), -1)
    # show the output image
    
else:
    cv2.imshow("output", np.hstack([img, output]))
    cv2.waitKey(0)

circles = np.round(circles[0, :]).astype("int")
circles.shape

"""# Watershed"""

img = cv2.imread('tests5.tif')
output = img.copy()
hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
hue = hsv[:,:,0]
sat = hsv[:,:,1]
val = hsv[:,:,2]

cv2.imshow("output",output)
cv2.waitKey(0)

ret, _ = cv2.threshold(sat,0,255,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)
ret, thresh = cv2.threshold(sat,ret,255,cv2.THRESH_BINARY_INV)

cv2.imshow("output",thresh)
cv2.waitKey(0)

# noise removal
kernel = np.ones((3,3),np.uint8)
opening = cv2.morphologyEx(thresh,cv2.MORPH_OPEN,kernel, iterations = 2)

cv2.imshow("output",opening)
cv2.waitKey(0)

# sure background area
sure_bg = cv2.dilate(opening,kernel,iterations=3)

cv2.imshow("output",sure_bg)
cv2.waitKey(0)

# Finding sure foreground area
dist_transform = cv2.distanceTransform(opening,cv2.DIST_L2,5)
ret, sure_fg = cv2.threshold(dist_transform,0.5*dist_transform.max(),255,0)

cv2.imshow("output",sure_fg)
cv2.waitKey(0)

# Finding unknown region
sure_fg = np.uint8(sure_fg)
unknown = cv2.subtract(sure_bg,sure_fg)

cv2.imshow("output",unknown)
cv2.waitKey(0)

# Marker labelling
ret, markers = cv2.connectedComponents(sure_fg)

# Add one to all labels so that sure background is not 0, but 1
markers = markers+1

# Now, mark the region of unknown with zero
markers[unknown==255] = 0

markers = cv2.watershed(img,markers)
output[markers == -1] = [255,0,0]

cv2.imshow("output",output)
cv2.waitKey(0)
cv2.destroyAllWindows()

start_time = int(round(time.time() * 1000))


# ++++++++++++++++++++++ Image processing ++++++++++++++++++++++

# Open image and create copy for output
img = cv2.imread('tests5_1.tif')
img_size = img.shape
img_rgb = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
output_markers = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
# Convert color space to HSV
hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
# Get individual HSV channels (hue, saturation, value)
hue = hsv[:,:,0]
sat = hsv[:,:,1]
val = hsv[:,:,2]
# Get individual RGB channels (red, blue, green)
red = img_rgb[:,:,0]
green = img_rgb[:,:,1]
blue = img_rgb[:,:,2]

# Create background mask (use the Value component of HSV since background is dark)
# Get optimum threshold by Otsu's binarization
opt_thresh, _ = cv2.threshold(val,0,255,cv2.THRESH_OTSU)
# Get the mask with half the optimum threshold
_, bck_mask = cv2.threshold(val,opt_thresh/2,255,cv2.THRESH_BINARY_INV)
# Filter noise
kernel = np.ones((5,5),np.uint8)
bck_mask = cv2.morphologyEx(bck_mask, cv2.MORPH_OPEN, kernel)          # Delete individual noise clusters
bck_mask = cv2.morphologyEx(bck_mask, cv2.MORPH_CLOSE, kernel)         # Close wholes in mask
bck_mask = cv2.dilate(bck_mask,kernel,iterations = 5)                  # Expand inwards to clean edge noise
# cv2.imshow("output",bck_mask)
# cv2.waitKey(0)

# Get high saturation mask
_, thresh_sat = cv2.threshold(sat,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
# Filter noise
thresh_sat = cv2.morphologyEx(thresh_sat, cv2.MORPH_OPEN, kernel)          # Delete individual noise clusters
thresh_sat = cv2.morphologyEx(thresh_sat, cv2.MORPH_CLOSE, kernel)         # Close wholes in mask
# cv2.imshow("output",thresh_sat)
# cv2.waitKey(0)

# Subtract background mask from Saturation mask
samples_mask = cv2.subtract(thresh_sat,bck_mask)
# cv2.imshow("output",samples_mask)
# cv2.waitKey(0)

# compute the exact Euclidean distance from every binary
# pixel to the nearest zero pixel, then find peaks in this
# distance map
eucl_dist = ndimage.distance_transform_edt(samples_mask)
localMax = peak_local_max(eucl_dist, indices=False, min_distance=10)

# perform a connected component analysis on the local peaks,
# using 8-connectivity, then aplpy the Watershed algorithm
markers = ndimage.label(localMax, structure=np.ones((3, 3)))[0]
labels = watershed(-eucl_dist, markers, mask=samples_mask)
num_samples = len(np.unique(labels)) - 1
samples_data = []
print("[INFO] {} samples detected".format(num_samples))

# loop over the unique labels returned by the Watershed algorithm
average_radius = 0
for idx,label in enumerate(np.unique(labels)):
    # if the label is zero, we are examining the 'background' so simply ignore it
    if label == 0:
        continue
    # otherwise, allocate memory for the label region and draw it on the mask
    mask = np.zeros(samples_mask.shape, dtype="uint8")
    mask[labels == label] = 1
    # detect contours in the mask and grab the largest one
    cnts = cv2.findContours(mask.copy(), cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)
    cnts = imutils.grab_contours(cnts)
    c = max(cnts, key=cv2.contourArea)
    # Get color statistics
    nonzero_idx = mask.nonzero()
    hue_avr = int(np.mean(mask[nonzero_idx] * hue[nonzero_idx]))
    sat_avr = int(np.mean(mask[nonzero_idx] * sat[nonzero_idx]))
    val_avr = int(np.mean(mask[nonzero_idx] * val[nonzero_idx]))
    red_avr = int(np.mean(mask[nonzero_idx] * red[nonzero_idx]))
    green_avr = int(np.mean(mask[nonzero_idx] * green[nonzero_idx]))
    blue_avr = int(np.mean(mask[nonzero_idx] * blue[nonzero_idx]))
    # Get result (positive:True or negative:False) using some heuristic (can change)
    result_state = hue_avr > 15
    
    # draw a circle enclosing the object
    ((x, y), r) = cv2.minEnclosingCircle(c)
    r = int(r)
    average_radius += r
    
    # Build sample list:
    # - sample index
    # - tuple with position in original image (x,y) and radius (r)
    # - tuple with HSV values (mean over all pixels)
    # - tuple with RGB values (mean over all pixels)
    samples_data.append((idx,(x,y,r),result_state,(hue_avr,sat_avr,val_avr),(red_avr,blue_avr,green_avr)))
    
# Calculate final radius average and multiply by scalar
average_radius /= num_samples
corrected_radius = int(average_radius * 0.9)

# Create image for results
y_offset = 50
x_offset = 50
y_spacing = 90
x_spacing = 200
num_cols = 10    # number of columns
num_cols = min(num_cols, int((GetSystemMetrics(0) - 2*x_offset)/x_spacing)) # prevent expaing outside monitor screen
num_rows = math.ceil(num_samples/num_cols)
# Create a new output image (with RGB channels)
output_results = np.ones((int(num_rows * y_spacing + y_offset), int(num_cols * x_spacing + x_offset), 3), np.uint8)*255
output_results_size = output_results.shape

# Sort samples by some feature
feature_map = {'index':0, 'position':1, 'result':2, 'HSV':3, 'RGB':4}
#samples_data.sort(key=lambda tup: not tup[feature_map['result']])   # oder by results (POSITIVE first)
#samples_data.sort(key=lambda tup: tup[feature_map['result']])   # oder by results (NEGATIVE first)
samples_data.sort(key=lambda tup: tup[feature_map['HSV']][0])

# Add markers to the "output_markers" image
r = corrected_radius
dec_pl = 2 # decimal places for RGB and HSV
img_s = 50 # samples image size
for for_idx, (idx,(x,y,ind_r), result_state, (hue_avr,sat_avr,val_avr),(red_avr,blue_avr,green_avr)) in enumerate(samples_data):
    x = int(x)
    y = int(y)
    cv2.circle(output_markers, (x, y), r, (0,255,0), 1)
    cv2.putText(output_markers, "#{}".format(idx), (x-12, y+3),cv2.FONT_HERSHEY_SIMPLEX, 0.4, (0, 0, 0), 1)

    # +++++++++++++++ Add samples to "output_results" +++++++++++++++ 
    
    # Define position of each sample in the final "output_results"
    x_idx = (for_idx % num_cols)
    y_idx = (math.floor(for_idx / num_cols))
    x_px = x_idx * x_spacing + x_offset
    y_px = y_idx * y_spacing + y_offset
    
    # Copy from original image to new "output_results"
    # Cut square with skirt pixels
    sr = int(average_radius*1.5)   # skirt pixel radius (larger than circle)
    img_temp = img_rgb[y-sr:y+sr+1,x-sr:x+sr+1,:]
    sf = img_s/(2*sr)   # scale factor
    # Resize sample patch
    rs_img = cv2.resize(img_temp, None, fx=sf, fy=sf, interpolation=cv2.INTER_CUBIC)
    # Transform variables according to scale factor
    sr = int(sr*sf)
    nr = int(r*1.8*sf)  # new circle radius
    xr, yr = np.arange(0,rs_img.shape[0]), np.arange(0,rs_img.shape[1])
    mask_cut = (xr[np.newaxis,:]-sr)**2 + (yr[:,np.newaxis]-sr)**2 < (nr)**2
    mask_cut = np.stack([mask_cut]*3,axis=2)
    img_temp = rs_img[mask_cut]
    xr, yr = np.arange(0,output_results_size[1]), np.arange(0,output_results_size[0])
    mask_paste = (xr[np.newaxis,:]-x_px)**2 + (yr[:,np.newaxis]-y_px)**2 < (nr)**2
    mask_paste = np.stack([mask_paste]*3,axis=2)
    output_results[mask_paste] = img_temp
    
    # Print statistics
    text_x_shift, text_y_shift = int(1.5*nr), -20
    start_pnt, end_point = (x_px + text_x_shift + 42,y_px+ text_y_shift + 3), (x_px + text_x_shift + 105,y_px+ text_y_shift + 19)
    cv2.putText(output_results, "#{} (Sample #{})".format(for_idx+1,idx), (x_px + text_x_shift, y_px + text_y_shift),cv2.FONT_HERSHEY_SIMPLEX, 0.4, (0, 0, 0), 1)
    cv2.putText(output_results, "Result:", (x_px + text_x_shift, y_px+ text_y_shift + 15),cv2.FONT_HERSHEY_SIMPLEX, 0.4, (0,0,0), 1)
    #cv2.rectangle(output_results, start_pnt, end_point, (10,10,10), -1)
    cv2.putText(output_results, "{}".format("POSITIVE" if result_state else "NEGATIVE"), (x_px + text_x_shift + 45, y_px+ text_y_shift + 15),cv2.FONT_HERSHEY_SIMPLEX, 0.4, (0, 255, 0) if result_state else (255,0,0), 1)
    cv2.putText(output_results, "Viral RNA: {}mM".format(0), (x_px + text_x_shift, y_px+ text_y_shift + 30),cv2.FONT_HERSHEY_SIMPLEX, 0.4, (0, 0, 0), 1)
    #cv2.putText(output_results, "HSV:({:03d},{:03d},{:03d})".format(hue_avr,sat_avr,val_avr), (x_px + text_x_shift, y_px+ text_y_shift + 45),cv2.FONT_HERSHEY_SIMPLEX, 0.4, (0, 0, 0), 1)
    #cv2.putText(output_results, "RGB:({:03d},{:03d},{:03d})".format(red_avr,green_avr,blue_avr), (x_px + text_x_shift, y_px+ text_y_shift +60),cv2.FONT_HERSHEY_SIMPLEX, 0.4, (0, 0, 0), 1)

    

    
# ++++++++++++++++ Write .csv file with all the information ++++++++++++++++ 

samples_data.sort(key=lambda tup: tup[feature_map['index']])  # re-order using some feature
try:
    with open('results.csv','w',newline='') as file:
        fieldnames = ['Sample ID', 'Result', 'Viral RNA conc.(mM)', 'Red','Green','Blue','Hue','Saturation','Value']
        thewriter = csv.DictWriter(file, fieldnames = fieldnames)
        thewriter.writeheader()
        for for_idx, (idx,(x,y,ind_r), result_state, (hue_avr,sat_avr,val_avr),(red_avr,blue_avr,green_avr)) in enumerate(samples_data):
            thewriter.writerow({
                'Sample ID': idx,
                'Result': 'POSITIVE' if result_state else 'NEGATIVE',
                'Viral RNA conc.(mM)': 0,
                'Red': red_avr,
                'Green': green_avr,
                'Blue': blue_avr,
                'Hue': hue_avr,
                'Saturation': sat_avr,
                'Value': val_avr
            })
except:
    print('File could not be opened!')

    
end_time = int(round(time.time() * 1000))
print("[INFO] Execution time: {}ms".format(end_time - start_time))

# Show final results
cv2.imshow("output",img)
cv2.waitKey(0) 
cv2.imshow("output",cv2.cvtColor(output_markers, cv2.COLOR_BGR2RGB))
cv2.waitKey(0)    
cv2.imshow("output",cv2.cvtColor(output_results, cv2.COLOR_BGR2RGB))
cv2.waitKey(0)
cv2.destroyAllWindows()

# # +++++++++ Hough Transform +++++++++ 
# cimg = img
# minCircleRadius = 5
# maxCircleRadius = 20
# circles = cv2.HoughCircles(sample_edges,cv2.HOUGH_GRADIENT,1,minCircleRadius*2,param1=50,param2=10,minRadius=minCircleRadius,maxRadius=maxCircleRadius)
# #circles = cv2.HoughCircles(samples,cv2.HOUGH_GRADIENT,1,2*minCircleRadius,param1=50,param2=30,minRadius=minCircleRadius,maxRadius=maxCircleRadius)
# if circles is not None:
#     circles = np.uint16(np.around(circles))
#     for i in circles[0,:]:
#         # draw the outer circle
#         cv2.circle(cimg,(i[0],i[1]),i[2],(0,255,0),2)
#         # draw the center of the circle
#         cv2.circle(cimg,(i[0],i[1]),2,(0,0,255),3)
#         print(i[2])

#     cv2.imshow('detected circles',cimg)
#     cv2.waitKey(0)
# else:
#     print("No circles detected...")
# cv2.destroyAllWindows()

fig, ax = plt.subplots(figsize=(18, 7))
ax.imshow(eucl_dist, cmap='bwr')
plt.axis('off')
plt.tight_layout()