How to get histogram of intensity of individual masked cells in an image?

Question

OK so newbie here that has been working on a set of homework problems with the original post here: How do I make a mask from one image and then transfer it to another? . The original idea was to take the DAPI image (grey image) and apply it as a mask to the NPM1 (green) image. After implementing the suggested code from HansHirse (thanks!) along with some other code I had been making for the homework problem I finally got a working histogram of all compatible cells in the image. The "compatibility" bit is that any cells touching the border weren't supposed to be counted. However, I still need to find a way to get histograms of each individual cell as well. I've attached the original images from the post too:

To do this, I tried blob_doh and one other method to get segmented regions of each cell but have no idea as to how I can apply these coordinates to an image for the histogram.

PS. The code is a bit messy. I segmented the code such that the blob_doh is near the bottom and the other method is also its own separate piece at the very bottom. Sorry!

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from skimage.feature import blob_dog, blob_log, blob_doh
from skimage.color import rgb2gray
import cv2
import mahotas as mh
import scipy
from scipy import ndimage
import matplotlib.patches as mpatches
from skimage import data
from skimage.filters import threshold_otsu
from skimage.segmentation import clear_border
from skimage.measure import label, regionprops
from skimage.morphology import closing, square
from skimage.color import label2rgb
# Read image into numpy array

image = cv2.imread("NOTREATDAPI.jpg",0)
dna = np.array(image) # must be gray-scale image
plt.show()

# Remove extraneous artifacts from image; set the threshold

dnaf = ndimage.gaussian_filter(dna, 8) #gaussian filter for general image
T = mh.thresholding.otsu(dnaf) # set threshold via mahotas otsu thresholding
theta=np.array(dnaf > T) #setting mask of values in image to calculated otsu threshold
cleared = clear_border(theta) #removes all cells that are in contact with the image border
epsilon = np.array(cleared) #final masked DAPI product

print("DAPI MASK USING GAUSSIAN FILTER AND OTSU THRESHOLDING"); 
plt.imshow(epsilon)
plt.show()



# Load and reset original images

image = cv2.imread("NOTREATDAPI.jpg",0) #The DAPI Image
image1 = cv2.imread("NOTREATNPM1.jpg",0) #The NPM1 Image

print("Original DAPI Image");plt.imshow(image);plt.show() #The DAPI Image
print("Original NPM1 Image");plt.imshow(image1);plt.show() #The NPM1 Image

# Create an array of bool of same shape as image

maskAboveThreshold = epsilon > 0 #Use mask array from above - include only values above non-masked zeros

print("Final Masked Image of NPM1"); plt.imshow(image1 * 
maskAboveThreshold, cmap='gray')
plt.show()

True_NPM1= image1 * maskAboveThreshold # Final masked version of NPM1 set back to grayscale 

# Create a mask using the DAPI image and binary thresholding at 25
_, mask = cv2.threshold(True_NPM1, 1, 255, cv2.THRESH_BINARY)

# Do some morphological opening to get rid of small artifacts
mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, 
cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (15, 15)))

# Calculate the histogram using the NPM1 image and the obtained binary 
mask

hist = cv2.calcHist([image1], [0], mask, [256], [0, 256])

# Show bar plot of calculated histogram

plt.bar(np.arange(256), np.squeeze(hist))
plt.show()

# Show mask image

plt.imshow(mask)
plt.show()




#blob_doh way of segmenting the cells ------

import cv2 as cv
from PIL import Image, ImageDraw

image10 = np.array(Image.open("OXALIDAPI.jpg"))

plt.imshow(image10)

#Convert to gaussian image with thresholds

image10 = cv2.imread("OXALIDAPI.jpg",0)
dna = np.array(image10) # gray-scale image
plt.show()

# Remove extraneous artifacts from image; set the threshold

dnaf = ndimage.gaussian_filter(dna, 8) #gaussian filter for general image
T = mh.thresholding.otsu(dnaf) # set threshold via mahotas otsu thresholding
theta=np.array(dnaf > T) #setting mask of values in image to calculated otsu threshold
cleared = clear_border(theta) #removes all cells that are in contact with the image border
image = np.array(cleared) #final masked DAPI product
#print("DAPI MASK USING GAUSSIAN FILTER AND OTSU THRESHOLDING"); 
plt.imshow(epsilon)
plt.show()

# Convert image to grayscale
image_gray = rgb2gray(image)
plt.imshow(image_gray,cmap="gray")

def plot_blobs(img,blobs):

    fig = plt.figure()
    ax = fig.add_subplot(1,1,1) 

    ax.imshow(img, interpolation='nearest')
    for blob in blobs:
        y, x, r = blob
        c = plt.Circle((x, y), r*1.25, color="red", linewidth=1, fill=False)
        ax.add_patch(c)

# blob_doh

blobs_doh = blob_doh(image_gray, min_sigma=10, max_sigma=256, 
threshold=.025)
plot_blobs(image,blobs_doh)

#get blob coordinates 

def filter_blobs(blobs,r_cutoff=5):

    new_blobs = []
    for b in blobs:
        if b[2] > r_cutoff:
            new_blobs.append(b)

return new_blobs

new_blobs = filter_blobs(blobs_doh)
#plot_blobs(image,new_blobs)

print(new_blobs)




#Other method of segmenting cells. maybe useful?


yeta = cv2.imread("NOTREATDAPI.jpg",0)
image = np.array(yeta)

# apply threshold
dnaf = ndimage.gaussian_filter(image, 8)
T = mh.thresholding.otsu(dnaf) # set threshold
plt.imshow(dnaf > T)
epsilon=np.array(dnaf > T)
plt.show()

# remove artifacts connected to image border
cleared = clear_border(epsilon)

# label image regions

label_image = label(cleared)
image_label_overlay = label2rgb(label_image, image=image)

fig, ax = plt.subplots(figsize=(6, 6))
ax.imshow(image_label_overlay)

for region in regionprops(label_image):
# take regions with large enough areas
    if region.area >= 50:
    # draw rectangle around individual cells
        minr, minc, maxr, maxc = region.bbox
        rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr,
                                  fill=False, edgecolor='red', linewidth=0.5)
    ax.add_patch(rect)

#ax.set_axis_off()
#plt.tight_layout()
plt.show()

howzer=np.array(image_label_overlay)

Answer 1

What you are looking for is cv2.connectedComponents. Basically, once you have the binary mask that separate the cells, you try to label each connected component of the mask as one cell:

# I choose OTSU instead of binary, but they are not much different in this case
_, mask = cv2.threshold(dapi, 25, 255, cv2.THRESH_OTSU)

# compute the connected component
labels, markers = cv2.connectedComponents(mask)

# load 2nd image in grayscale
# as your 2nd image is only green/black
npm1 = cv2.imread('npm1.jpg', cv2.IMREAD_GRAYSCALE)

# for you image (and usually), labels[0] is the background
for label in labels[1:]:
    # compute the histogram over the entire 256 levels of intensity
    hist, bins = np.histogram(npm1[markers==label], bins=range(256))

    # do whatever you like to hist
    # note that bins=range(256) and hist only have 255 values
    plt.bar(bins[1:], hist)
    plt.title('cell number: {:}'.format(label))

So for example the histogram of the first and second cells:

And the cell markers are: