Как получить гистограмму интенсивности отдельных замаскированных клеток на изображении? - PullRequest
1 голос
/ 06 июня 2019

ОК, так что новичок здесь, который работал над набором домашних заданий с оригинальным постом здесь: Как сделать маску из одного изображения, а затем перенести ее в другое? , Первоначальная идея состояла в том, чтобы взять изображение DAPI (серое изображение) и применить его в качестве маски к изображению NPM1 (зеленое). После реализации предложенного кода от HansHirse (спасибо!) Вместе с другим кодом, который я делал для решения домашней задачи, я, наконец, получил рабочую гистограмму всех совместимых ячеек в изображении. Бит «совместимости» заключается в том, что любые ячейки, соприкасающиеся с границей, не должны учитываться. Однако мне все еще нужно найти способ получить гистограммы для каждой отдельной ячейки. Я приложил оригинальные изображения из поста тоже: DAPI imageNPM1 Image

Для этого я попробовал blob_doh и еще один метод, чтобы получить сегментированные области каждой ячейки, но понятия не имел, как я могу применить эти координаты к изображению для гистограммы.

PS. Код немного грязный. Я сегментировал код так, что 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)

1 Ответ

1 голос
/ 06 июня 2019

То, что вы ищете, это 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))

Например, гистограмма первой и второй ячеек:

enter image description here

enter image description here

И маркеры ячейки:

[enter image description here

...