python - 如何获得图像中单个蒙版细胞强度的直方图？

Question

好吧，这里的新手一直在用这里的原始帖子解决一系列家庭作业问题:How do I make a mask from one image and then transfer it to another?.最初的想法是获取 DAPI 图像(灰色图像)并将其作为蒙版应用于 NPM1(绿色)图像。在实现了 HansHirse 推荐的代码(谢谢!)以及我为家庭作业问题制作的其他一些代码之后，我终于得到了图像中所有兼容单元格的工作直方图。 “兼容性”位是任何接触边界的单元格都不应该被计算在内。但是，我仍然需要找到一种方法来获取每个单独单元格的直方图。我也附上了帖子中的原始图片:

为此，我尝试了 blob_doh 和另一种方法来获取每个单元格的分段区域，但不知道如何将这些坐标应用于直方图的图像。

附言。代码有点乱。我对代码进行了分段，使 blob_doh 位于底部附近，而另一种方法也是位于最底部的独立部分。对不起!

%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

您正在寻找的是cv2.connectedComponents。基本上，一旦有了分隔单元格的二元掩码，就可以尝试将掩码的每个连接组件标记为一个单元格:

# 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))

例如第一个和第二个单元格的直方图:

细胞标记物是: