2022/03/06 医学图像处理记录

目标：

1. 将多分类的mask中的细胞核与细胞质分别提取并将背景二值化
2. 生成灰度直方图
3. 根据灰度直方图拟合曲线找到峰谷值，分段计算占比核质的比例

材料：

1. 淋巴细胞原图img
   

2. 淋巴细胞图对应的mask
   

代码实现

1. 背景纯黑

import os
import cv2
import numpy as np

def add_mask2image_binary(images_path, masks_path, masked_path):
# Add binary masks to images
    for img_item in os.listdir(images_path):
        print(img_item)
        img_path = os.path.join(images_path, img_item)
        img = cv2.imread(img_path)
        mask_path = os.path.join(masks_path, img_item[:-4]+'.png')  # mask是.png格式的，image是.jpg格式的
        mask = cv2.imread(mask_path, cv2.IMREAD_GRAYSCALE)  # 将彩色mask以二值图像形式读取
        masked = cv2.add(img, np.zeros(np.shape(img), dtype=np.uint8), mask=mask)  #将image的相素值和mask像素值相加得到结果
        cv2.imwrite(os.path.join(masked_path, img_item), masked)
images_path = '20220116/img/'
masks_path = '20220116/img_out'
masked_path = '20220302/BlackBG/'
add_mask2image_binary(images_path, masks_path, masked_path)

效果

但没有实现核质分离，但是通过mask的灰度去手动获得核与质的灰度，即可分别得到核质与全黑背景

import os
import cv2
import numpy as np

def add_mask2image_binary(images_path, masks_path, masked1_path, masked2_path):
# Add binary masks to images
    for img_item in os.listdir(images_path):
        print(img_item)
        img_path = os.path.join(images_path, img_item)
        img = cv2.imread(img_path)
        mask_path = os.path.join(masks_path, img_item[:-4]+'.png')  # mask是.png格式的，image是.jpg格式的
        # mask = cv2.imread(mask_path, cv2.IMREAD_GRAYSCALE)  # 将彩色mask以二值图像形式读取
        mask = cv2.imread(mask_path, 0)
        # 取核
        upper1 = mask <= 80
        lower1= mask >= 70
        thresh1 = (np.multiply(upper1,lower1).astype(np.float32)*255).astype(np.uint8)
        masked1 = cv2.add(img, np.zeros(np.shape(img), dtype=np.uint8), mask = thresh1)  
        # 取质
        upper2 = mask <= 40
        lower2 = mask >= 30
        thresh2 = (np.multiply(upper2,lower2).astype(np.float32)*255).astype(np.uint8)
        masked2 = cv2.add(img, np.zeros(np.shape(img), dtype=np.uint8), mask = thresh2)
        
        cv2.imwrite(os.path.join(masked1_path, img_item), masked1)
        cv2.imwrite(os.path.join(masked2_path, img_item), masked2)
images_path = '20220116/img/'
masks_path = '20220116/img_out'
masked1_path = '20220302/test/1'
masked2_path = '20220302/test/2'
add_mask2image_binary(images_path, masks_path, masked1_path, masked2_path)

可视化

使用opencv的colormap函数将灰度转换成伪彩图

# https://www.pianshen.com/article/19011178284/
import cv2, csv, os
import numpy as np

imgsrc = '20220302/test/1/'
out_dir = '20220302/test/'
img_lst = []
filelist = os.listdir(imgsrc)
evaluate = [] # 评价指标
for file in filelist:
    filename = os.path.splitext(file)[0]
    # 注意图片后缀
    img = cv2.imread(imgsrc + '{}.jpg'.format(filename), 0)
    # 生成彩图并生成图片
    img_color = cv2.applyColorMap(img, cv2.COLORMAP_JET)
    # 图片输出路径
    theme = 'jet'
    
    cv2.imwrite(out_dir + '{}/{}.png'.format(theme,filename), img_color) 
    # 像素均值、方差
    means, dev = cv2.meanStdDev(img)
    print('means: {}  \n dev: {}'.format(means[0,0], dev[0,0]))
    
    hist = cv2.calcHist([img],[0],None,[256],[0,255])

    means_25 ,devs_25 = 0, 0 
    # 25% 64   /  50% 128 / 75%  192
    for x in range(0,64):
        means_25 += hist[x][0]
        devs_25 += hist[x][0]**2 
    mean_25 = means_25/64
    dev_25 = devs_25/64-mean_25**2
    print(mean_25, dev_25)


    means_50 ,devs_50 = 0, 0 
    for x in range(0,128):
        means_50 += hist[x][0]
        devs_50 += hist[x][0]**2 
    mean_50 = means_50/128
    dev_50 = devs_50/128-mean_50**2
    print(mean_50, dev_50)


    means_75 ,devs_75 = 0, 0 
    for x in range(0,192):
        means_75 += hist[x][0]
        devs_75 += hist[x][0]**2 
    mean_75 = means_75/192
    dev_75 = devs_75/192-mean_75**2
    print(mean_75, dev_75)

    # 存储评价指标
    evaluate.append((mean_25, mean_50, mean_75, means[0,0], dev_25, dev_50, dev_75, dev[0,0]))

header1 = ['means', '', '', '', 'dev' , '', '', '']
header2 = ['25%', '50%', '75%', '100%']
with open(out_dir + '{}/{}.csv'.format(theme, theme),'w') as csvfile:
    writer = csv.writer(csvfile)
    writer.writerow(header1) # 这里要以list形式写入，writer会在新建的csv文件中，一行一行写入
    writer.writerow(header2*2)
    writer.writerows(evaluate)

核质分离计算占比面积

多项式拟合灰度曲线，利用倒数来确定曲线的极值点。

import pandas as pd 
import matplotlib.pyplot as plt 
import numpy as np 
from scipy import signal 
import cv2
from lmfit import Model

imgsrc = 'test/bd/Snipaste_2022-03-07_12-16-36.png'  # 需要图片背景全黑
img = cv2.imread(imgsrc, 0)
rows,cols = img.shape
hist = cv2.calcHist([img],[0],None,[256],[0,255])
CellArea = rows * cols - hist[0][0]
y = []
for value in hist:
    y.append(int(value[0]))
y.pop(0) # 去掉黑色的像素点
y.pop(0)
x = np.arange(0, 254, 1)

# plan A :用n次多项式拟合
n = 7
z1 = np.polyfit(x, y, n)  
p1 = np.poly1d(z1) #多项式系数 print(p1)

yd = np.polyder(p1,1) # 1表示一阶导 
d = []
for i in yd.r:
    d.insert(0, np.round(i))
print(d)
print(y[0], hist[0], hist[1])
# 在屏幕上打印拟合多项式 
yvals=p1(x) 
plt.plot(x, y, '.',label='original values') 
plt.plot(x, yvals, 'r',label='polyfit values') 
plt.xlabel('x axis') 
plt.ylabel('y axis') 
plt.legend(loc=4) 
plt.title('polyfitting {}'.format(n)) 
plt.show()

输出的极值点的横轴值

[23.0, 89.0, 115.0, 175.0, 214.0, 242.0]

添加一下代码计算分段区间占比
步骤：

1. 拟合曲线找到极值点
2. 取两极值点x轴的中点
3. 统计各区段像素个数
4. 计算占比大小

    # 生成分类区段 , 两个极值点取中值
    x_subsection = []
    x_subsection.append(round(d[0]/2))
    for d1 in range(0, len(d) - 1, 1):
        x_subsection.append(round((d[d1] + d[d1+1])/2))
    print(x_subsection)
    
    # 区段像素统计量
    SegmentPixelStatistics= [] 
    for w in range(0, len(x_subsection) - 1, 1):
        SegmentPixel = 0
        for q in range(x_subsection[w], x_subsection[w+1]):
            SegmentPixel += int(hist[q])
        SegmentPixelStatistics.append(SegmentPixel)
    print(SegmentPixelStatistics)
    
    proportion = []
    for p in SegmentPixelStatistics:
        q = float(p/CellArea)
        b = "%.2f%%" % (q * 100)
        proportion.append(b)
    print(proportion)
    
    # 区段名
    section = []

    for s1 in range(len(x_subsection) - 1):
        section.append('{}-{}'.format(x_subsection[s1], x_subsection[s1+1]))
        
    single_dic = dict(zip(section, proportion))
    print(single_dic)
    
    # 在屏幕上打印拟合多项式 
    plt.figure()
    yvals=p1(x) 
    plt.plot(x, y, '.',label='original values') 
    plt.plot(x, yvals, 'r',label='polyfit values') 
    plt.xlabel('x axis') 
    plt.ylabel('y axis') 
    plt.legend(loc=4) 
    plt.title('{} polyfitting {}'.format(filename, n)) 
    plt.show()
    
    with open(out_dir + 'nucleus.csv','a') as csvfile:
        writer = csv.writer(csvfile)
        writer.writerow([filelist[num],'','','',''])
        writer.writerow(section)
        writer.writerow(proportion)
        
    num += 1