【问题标题】:How to calculate nucleus amount of cell?如何计算细胞核量?
【发布时间】:2019-01-08 08:09:56
【问题描述】:

我正在使用 Python 3.5 和 OpenCV 3 来分析生物学中的细胞图片。我的照片是这样的:

我希望能够计算细胞核面积与整个细胞面积的比率。

在我的幻灯片中,细胞核为深紫色,细胞的其他区域为浅蓝色。还有一些我想完全忽略的棕褐色红细胞。为清楚起见,这是一张带标签的图片:

如何使用图像分割来识别和测量我的感兴趣区域?

我已经尝试关注this guide,但它返回了一个全黑的图像。

【问题讨论】:

  • 您能发布您尝试使用的代码吗?该指南的链接很好,但我们仍在猜测可能做了什么导致看到的问题。

标签: python algorithm opencv image-processing image-segmentation


【解决方案1】:
# light purple color segmentation (to get cells)
cell_hsvmin = (110,40,145)
cell_hsvmax = (150,190,255)

hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)
color_thresh = cv2.inRange(hsv, cell_hsvmin, cell_hsvmax)
# masked = cv2.bitwise_and(frame,frame, mask=color_thresh)
# cv2.imshow('masked0', masked)
ksize = 5
open_thresh = cv2.morphologyEx(color_thresh, cv2.MORPH_OPEN, np.ones((ksize,ksize),'uint8'), iterations=1)
masked = cv2.bitwise_and(frame,frame, mask=open_thresh)
cv2.imshow('masked', masked)

# dark purple color segmentation (to get nucleus)
nucleus_hsvmin = (125,65,160)
nucleus_hsvmax = (150,190,255)

nucleus_color_thresh = cv2.inRange(hsv, nucleus_hsvmin, nucleus_hsvmax)
ksize = 3
nucleus_open_thresh = cv2.morphologyEx(nucleus_color_thresh, cv2.MORPH_OPEN, np.ones((ksize,ksize),'uint8'), iterations=1)
nucleus_masked = cv2.bitwise_and(masked,masked, mask=nucleus_open_thresh)
cv2.imshow('nucleus_masked', nucleus_masked)

"""
HULL APPROXIMATES THE CELLS TO A CIRCLE TO FILL IN GAPS CREATED BY THRESHOLDING AND CLOSING.
FOR NON-CIRCULAR CELLS LIKE IN YOUR SECOND IMAGE, THIS MIGHT CAUSE BAD AREA CALCULATIONS
"""
# doHULL = False
doHULL = True

cells = []
cells_ratio = []
minArea = frame.shape[0]*frame.shape[1]* 0.01
_, contours, _ = cv2.findContours(open_thresh,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)
for cnt in contours:
    area = cv2.contourArea(cnt)
    if area >= minArea:
        cells.append(cnt)
        nucleus_mask = np.zeros(frame.shape[:2], 'uint8')
        if doHULL:
            hull = cv2.convexHull(cnt)
            cv2.drawContours(nucleus_mask, [hull], -1, 255, -1);
        else:
            cv2.drawContours(nucleus_mask, [cnt], -1, 255, -1);
        nucleus_mask = cv2.bitwise_and(nucleus_open_thresh, nucleus_mask)
        nucleus_area = np.count_nonzero(nucleus_mask)

        ratio = nucleus_area / area
        cells_ratio.append(ratio)

        # nucleus_img = cv2.bitwise_and(frame, frame, mask=nucleus_mask)
        # cv2.imshow('nucleus_img', nucleus_img)
        # cv2.waitKey(0)

doDRAWCELLS = False
# doDRAWCELLS = True
if doDRAWCELLS:
    for cell_cnt in cells:
        cells_mask = np.zeros(frame.shape[:2], 'uint8')
        if doHULL:
            hull = cv2.convexHull(cell_cnt)
            cv2.drawContours(cells_mask, [hull], -1, 255, -1);
        else:
            cv2.drawContours(cells_mask, [cell_cnt], -1, 255, -1);
        cells_img = cv2.bitwise_and(frame, frame, mask=cells_mask)
        cv2.imshow('cells_img', cells_img)
        cv2.waitKey(0)

这仅适用于未连接的单元格。您可以以此为基础来使用分水岭算法。 此外,颜色分割参数已根据您发布的 2 张图像进行了调整。其他幻灯片可能会偏离颜色范围,因此您可能需要调整它们。如果调整它们不能让您取得良好的折衷,您可能需要研究 otsu 二值化或自适应阈值来分割颜色。

另一种选择是查看 cv2.MORPH_GRADIENT,它的工作原理类似于边缘检测器。或

gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
kernel = np.array([[1,1,1],[1,-8,1],[1,1,1]],dtype='float32')
laplace = cv2.filter2D(cv2.GaussianBlur(gray,(blur_ksize,blur_ksize),0), -1, kernel)
cv2.imshow('laplace', laplace)

并使用边缘来分割单元格?

【讨论】:

    【解决方案2】:

    首先,我们将在下面使用一些初步代码:

    import numpy as np
    import cv2
    from matplotlib import pyplot as plt
    from skimage.morphology import extrema
    from skimage.morphology import watershed as skwater
    
    def ShowImage(title,img,ctype):
      if ctype=='bgr':
        b,g,r = cv2.split(img)       # get b,g,r
        rgb_img = cv2.merge([r,g,b])     # switch it to rgb
        plt.imshow(rgb_img)
      elif ctype=='hsv':
        rgb = cv2.cvtColor(img,cv2.COLOR_HSV2RGB)
        plt.imshow(rgb)
      elif ctype=='gray':
        plt.imshow(img,cmap='gray')
      elif ctype=='rgb':
        plt.imshow(img)
      else:
        raise Exception("Unknown colour type")
      plt.title(title)
      plt.show()
    

    作为参考,这是您的原始图片:

    #Read in image
    img         = cv2.imread('cells.jpg')
    ShowImage('Original',img,'bgr')
    

    您链接到的文章建议使用Otsu's method 进行颜色分割。该方法假设图像像素的强度可以绘制成双峰直方图,并为该直方图找到最佳分隔符。我应用下面的方法。

    #Convert to a single, grayscale channel
    gray        = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    #Threshold the image to binary using Otsu's method
    ret, thresh = cv2.threshold(gray,0,255,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)
    ShowImage('Grayscale',gray,'gray')
    ShowImage('Applying Otsu',thresh,'gray')
    

    图像的二进制形式不是很好!查看灰度图像,您可以了解原因:Otsu 变换会产生 类像素:深色背景像素、甜甜圈细胞和细胞内部以及细胞核。下面的直方图说明了这一点:

    #Make a histogram of the intensities in the grayscale image
    plt.hist(gray.ravel(),256)
    plt.show()
    

    因此,您已经打破了您所使用算法的假设,因此得到不好的结果也就不足为奇了。通过丢弃颜色信息,我们失去了区分甜甜圈和细胞内部的能力。

    解决此问题的一种方法是根据颜色阈值执行分割。为此,您需要选择一个色彩空间来工作。This guide 对不同空间有出色的图片描述。

    让我们选择 HSV。这样做的好处是,单个通道H 描述了图像的颜色。一旦我们将图像转换到这个空间,我们就可以找到我们感兴趣的颜色的边界。例如,要找到细胞核,我们可以这样做:

    cell_hsvmin  = (110,40,145)  #Lower end of the HSV range defining the nuclei
    cell_hsvmax  = (150,190,255) #Upper end of the HSV range defining the nuclei
    #Transform image to HSV color space
    hsv          = cv2.cvtColor(img,cv2.COLOR_BGR2HSV) 
    #Threshold based on HSV values
    color_thresh = cv2.inRange(hsv, cell_hsvmin, cell_hsvmax) 
    ShowImage('Color Threshold',color_thresh,'gray')
    
    masked = cv2.bitwise_and(img,img, mask=color_thresh)
    ShowImage('Color Threshold Maksed',masked,'bgr')
    

    这看起来好多了!不过,请注意细胞内部的某些部分被标记为细胞核,即使它们不应该这样做。也有人会争辩说,这不是很自动:您仍然必须小心亲自挑选颜色。在 HSV 空间中操作消除了很多猜测,但也许我们可以利用有四种不同颜色的事实来消除对范围的需求!为此,我们通过 k-means clustering algorithm 传递 HSV 像素。

    #Convert pixel space to an array of triplets. These are vectors in 3-space.
    Z = hsv.reshape((-1,3)) 
    #Convert to floating point
    Z = np.float32(Z)
    #Define the K-means criteria, these are not too important
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)
    #Define the number of clusters to find
    K = 4
    #Perform the k-means transformation. What we get back are:
    #*Centers: The coordinates at the center of each 3-space cluster
    #*Labels: Numeric labels for each cluster
    #*Ret: A return code indicating whether the algorithm converged, &c.
    ret,label,center = cv2.kmeans(Z,K,None,criteria,10,cv2.KMEANS_RANDOM_CENTERS)
    
    #Produce an image using only the center colours of the clusters
    center = np.uint8(center)
    khsv   = center[label.flatten()]
    khsv   = khsv.reshape((img.shape))
    ShowImage('K-means',khsv,'hsv')
    
    #Reshape labels for masking
    label = label.reshape(img.shape[0:2])
    ShowImage('K-means Labels',label,'gray')
    

    请注意,这在分离颜色方面做得非常出色,无需手动指定! (除了指定簇的数量。)

    现在,我们需要弄清楚哪些标签对应于单元格的哪些部分。

    为此,我们找到两个像素的颜色:一个明显是细胞核像素,另一个明显是细胞像素。然后我们会找出离这些像素最近的聚类中心。

    #(Distance,Label) pairs
    nucleus_colour = np.array([139, 106, 192])
    cell_colour    = np.array([130, 41,  207])
    nuclei_label  = (np.inf,-1)
    cell_label    = (np.inf,-1)
    for l,c in enumerate(center):
      print(l,c)
      dist_nuc = np.sum(np.square(c-nucleus_colour)) #Euclidean distance between colours
      if dist_nuc<nuclei_label[0]:
            nuclei_label=(dist_nuc,l)
      dist_cell = np.sum(np.square(c-cell_colour)) #Euclidean distance between colours
      if dist_cell<cell_label[0]:
            cell_label=(dist_cell,l)
    nuclei_label = nuclei_label[1]
    cell_label   = cell_label[1]
    print("Nuclei label={0}, cell label={1}".format(nuclei_label,cell_label))
    

    现在,让我们构建我们需要的二元分类器,以便为分水岭算法识别整个单元格:

    #Multiply by 1 to keep image in an integer format suitable for OpenCV
    thresh = cv2.bitwise_or(1*(label==nuclei_label),1*(label==cell_label))
    thresh = np.uint8(thresh)
    ShowImage('Binary',thresh,'gray')
    

    我们现在可以消除单像素噪声:

    #Remove noise by eliminating single-pixel patches
    kernel  = np.ones((3,3),np.uint8)
    opening = cv2.morphologyEx(thresh,cv2.MORPH_OPEN, kernel, iterations = 2)
    ShowImage('Opening',opening,'gray')
    

    我们现在需要识别分水岭的山峰并给它们单独的标签。这样做的目的是生成一组像素,使得每个细胞核+细胞中都有一个像素,并且没有两个细胞核的标识符像素接触。

    为了实现这一点,我们可以执行距离变换,然后过滤掉距离细胞核+细胞中心两个远的距离。

    但是,我们必须小心,因为具有高阈值的长而窄的单元格可能会完全消失。在下图中,我们将右下角的两个相互接触的单元格分开,但完全消除了右上角的长而窄的单元格。

    #Identify areas which are surely foreground
    fraction_foreground = 0.75
    dist         = cv2.distanceTransform(opening,cv2.DIST_L2,5)
    ret, sure_fg = cv2.threshold(dist,fraction_foreground*dist.max(),255,0)
    ShowImage('Distance',dist_transform,'gray')
    ShowImage('Surely Foreground',sure_fg,'gray')
    

    降低阈值会使又长又窄的单元格返回,但会使右下角的单元格保持连接状态。

    我们可以通过使用识别每个局部区域的峰值的自适应方法来解决此问题。这消除了为我们的阈值设置单个全局常量的需要。为此,我们使用h_axima 函数,该函数返回所有大于指定截止值的局部最大值。这与距离函数相反,后者返回大于给定值的所有像素。

    #Identify areas which are surely foreground
    h_fraction = 0.1
    dist     = cv2.distanceTransform(opening,cv2.DIST_L2,5)
    maxima   = extrema.h_maxima(dist, h_fraction*dist.max())
    print("Peaks found: {0}".format(np.sum(maxima)))
    #Dilate the maxima so we can see them
    maxima   = cv2.dilate(maxima, kernel, iterations=2)
    ShowImage('Distance',dist_transform,'gray')
    ShowImage('Surely Foreground',maxima,'gray')
    

    现在我们通过减去最大值来识别未知区域,将由分水岭算法标记的区域:

    # Finding unknown region
    unknown = cv2.subtract(opening,maxima)
    ShowImage('Unknown',unknown,'gray')
    

    接下来,我们为每个最大值赋予唯一标签,然后在最终执行分水岭变换之前标记未知区域:

    # Marker labelling
    ret, markers = cv2.connectedComponents(maxima)
    ShowImage('Connected Components',markers,'rgb')
    
    # 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==np.max(unknown)] = 0
    
    ShowImage('markers',markers,'rgb')
    
    dist    = cv2.distanceTransform(opening,cv2.DIST_L2,5)
    markers = skwater(-dist,markers,watershed_line=True)
    
    ShowImage('Watershed',markers,'rgb')
    imgout = img.copy()
    imgout[markers == 0] = [0,0,255] #Label the watershed_line
    
    ShowImage('img',imgout,'bgr')
    

    这为我们提供了一组代表细胞的标记区域。接下来,我们遍历这些区域,将它们用作标记数据的掩码,并计算分数:

    for l in np.unique(markers):
        if l==0:      #Watershed line
            continue
        if l==1:      #Background
            continue
        #For displaying individual cells
        #temp=khsv.copy()
        #temp[markers!=l]=0
        #ShowImage('out',temp,'hsv')
        temp = label.copy()
        temp[markers!=l]=-1
        nucleus_area = np.sum(temp==nuclei_label)
        cell_area    = np.sum(temp==cell_label)
        print("Nucleus fraction for cell {0} is {1}".format(l,nucleus_area/(cell_area+nucleus_area)))
    

    这给出了:

    Nucleus fraction for cell 2 is 0.9002795899347623
    Nucleus fraction for cell 3 is 0.7953321364452424
    Nucleus fraction for cell 4 is 0.7525925925925926
    Nucleus fraction for cell 5 is 0.8151515151515152
    Nucleus fraction for cell 6 is 0.6808656818962556
    Nucleus fraction for cell 7 is 0.8276481149012568
    Nucleus fraction for cell 8 is 0.878500237304224
    Nucleus fraction for cell 9 is 0.8342518016108521
    Nucleus fraction for cell 10 is 0.9742324561403509
    Nucleus fraction for cell 11 is 0.8728733459357277
    Nucleus fraction for cell 12 is 0.7968570333461096
    Nucleus fraction for cell 13 is 0.8226831716293075
    Nucleus fraction for cell 14 is 0.7491039426523297
    Nucleus fraction for cell 15 is 0.839096357768557
    Nucleus fraction for cell 16 is 0.7589670014347202
    Nucleus fraction for cell 17 is 0.8559168925022583
    Nucleus fraction for cell 18 is 0.7534142640364189
    Nucleus fraction for cell 19 is 0.8036734693877551
    Nucleus fraction for cell 20 is 0.7566037735849057
    

    (请注意,如果您将其用于学术目的,学术诚信需要适当的署名。请联系我了解详情。)

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2010-10-10
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2023-01-26
      • 1970-01-01
      相关资源
      最近更新 更多