Python----opencv计算细胞核质比、轮廓描绘、椭圆拟合、GUI手动调整阈值、自定义函数自动检测阈值

Python----opencv识别细胞核质、手动调整阈值与阈值的自动检测

一、题目:根据附件 cell.jpg,使用 opencv 库或者 PIL 库计算细胞核与细胞质的面积比。

在这里插入图片描述

二、逐步分解,逐步攻破

1.借助Photoshop初步分析图片
1.1 借助魔棒工具或快速选择工具

1.2 两次分别选择细胞核和整个细胞

在这里插入图片描述

1.3 借助Photoshop的直方图的明度(转化为灰度图)来观察

在这里插入图片描述

1.4初步得出结论

在这里插入图片描述

2. 借助刚刚得到的阈值和opencv库得到细胞的二值图
cell = cv.cvtColor(cell_o, cv.COLOR_BGR2GRAY)

rows, cols = cell_o.shape[:2]  # 获得图片的高和宽
    mapx = np.zeros(cell_o.shape[:2], np.float32)  # mapx是一个图片大小的变换矩阵
    mapy = np.zeros(cell_o.shape[:2], np.float32)  # mapy也是一个图片大小的变换矩阵
    for i in range(rows):  # 遍历,将获得到图片的宽和稿依次赋予两个变换矩阵
        for j in range(cols):
            mapx.itemset((i, j), j)
            mapy.itemset((i, j), i)
            
ret, cell_2v_whole = cv.threshold(cell, x1, 255, cv.THRESH_BINARY_INV)  # 获取整个细胞的二值图
ret, cell_2v_nucleus = cv.threshold(cell, x2, 255, cv.THRESH_BINARY_INV)  # 获取细胞核的二值图
cv.namedWindow("cell_whole", 0)
cv.resizeWindow("cell_whole", 500, 500)
cv.imshow("cell_whole", cell_2v_whole)
cv.namedWindow("cell_nucleus", 0)
cv.resizeWindow("cell_nucleus", 500, 500)
cv.imshow("cell_nucleus", cell_2v_nucleus)

在这里插入图片描述

3.借助opencv描绘出二值图的最大轮廓,并用椭圆拟合

​ 这一步使用opencv处理细胞核质比的关键,通过Photoshop得到的两个阈值,我们已经基本能够得到细胞核和细胞质的轮廓的。但是通过观察上面的二值图可以发现,其实仅仅通过这两个阈值来用opencv提取的话,会有很多细胞质中的深度点被误当成细胞核的部分,因此如果直接用二值图的面积来计算的话,会导致结果及其不准确。

  • 提取最大轮廓

    因为细胞核、细胞质是成块出现的,所以在我们只需要找到最大的一个封闭的轮廓就可以了

  • 椭圆拟合

    ​ 使用椭圆拟合最大的轮廓,是为了让计算更加精确,opencv是有直接的方法来做到这一点的,而椭圆的面积是有直接的公式可以计算的,即:S=Π×a×b。

    ​ 并且使用椭圆拟合后的面积与最大轮廓面积的误差值可以作为一个新的十分重要的数据存在。

import cv2 as cv

# 基于二值图像用外轮廓的模式,通过全点连接轮廓的方法提取轮廓
contours_whole, hierarchy = cv.findContours(cell_2v_whole, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_NONE)
contours_nucleus, hierarchy = cv.findContours(cell_2v_nucleus, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_NONE)

# 获取目标图像的最小矩阵,找到最大的轮廓
maxArea = 0
maxContour = {}
for item in contours_whole:
    x, y, w, h = cv.boundingRect(item)
    if w * h > maxArea:
        maxArea = w * h
        maxContour = item
whole_area = 0
whole_area += cv.contourArea(maxContour)

img = cv.drawContours(cell_o.copy(), maxContour, -1, (0, 255, 0), 1)

cnt = maxContour
ellipse = cv.fitEllipse(cnt)
(x, y), (a, b), angle = cv.fitEllipse(cnt)
cv.ellipse(img, ellipse, (255, 0, 0), 1)
whole_area_ellipse = math.pi * a * b / 4

# 重复,查找细胞核的最大轮廓
maxArea = 0
maxContour = {}
for item in contours_nucleus:
    x, y, w, h = cv.boundingRect(item)
    if w * h > maxArea:
        maxArea = w * h
        maxContour = item
nucleus_area = 0
nucleus_area += cv.contourArea(maxContour)
img = cv.drawContours(img, maxContour, -1, (0, 0, 255), 1)

cnt = maxContour
ellipse = cv.fitEllipse(cnt)
(x, y), (a, b), angle = cv.fitEllipse(cnt)
cv.ellipse(img, ellipse, (255, 0, 0), 1)
nucleus_area_ellipse = math.pi * a * b / 4

在这里插入图片描述

4. 借助Python的GUI界面手动调整阈值,观察数据

在这里插入图片描述

在这里插入图片描述

5.绘制不同阈值下,椭圆拟合面积和最大轮廓面积误差值折线图
for i in range(0, 256):
    scaleX1.set(i)
    scaleX2.set(100)
    try:
        whole_area, nucleus_area, whole_area_ellipse, nucleus_area_ellipse = runDetect(i, 100)
        ls.append(abs(whole_area - whole_area_ellipse))
    except Exception:
        ls.append(0.1)
        
# runDetect函数就是绘制最大轮廓和椭圆拟合的函数,3、4都是其部分代码
# 这里是将阈值从0-255遍历一次,获得各阈值下的核质误差值加入到一个列表中,用于绘制折线图和之后的分析。

在这里插入图片描述

6. 根据误差值折线图特性,自定义函数自动查找两个阈值

  • 分析两个正确阈值的数据特点
    • 误差值很小,趋近于0
    • 处于误差值的上升突变点
    • 处于三个峰的右边两个峰的左侧
  • 转化为数值
    • right > mean * left
    • left != 0
    • abs(ls[i] - ls[i -1]) > 0.01
    • 0 < ls[i] < mean2
def slopeComparison(i,ls,ls1):
    left = abs(ls[i] - ls[i - 1]) + abs(ls[i-1] - ls[i-2])		# 左侧变化率
    right = ls[i+2] - ls[i]		# 右侧变化率
    mean = np.mean(ls1)		# 误差值变化率平均值
    mean2 = np.mean(ls)			# 误差值平均值
    if right > mean * left and left != 0 and (abs(ls[i] - ls[i - 1]) > 0.01) and 0 < ls[i] < mean2:
        return True
    else:
        return False

在这里插入图片描述

三、程序测试

对于大部分单细胞都可以比较准确地进行自动识别阈值并求出核质比,如果有一定误差,也可以通过手动调整来实现准确的结果。

在这里插入图片描述

但是对于部分图片,自动检测阈值会失效,有时候只能检测到一个符合函数的阈值,毕竟这个函数设置的太简单了,还需要进一步加强科学性和准确性。

在这里插入图片描述

四、相关的问题和拓展

曲线平滑处理
  • 使用方差、标准差、最值来逐端平滑曲线
  • 使用部分返回平均值代替单个点特殊值
  • 使用scipy库的savgol_filter方法进行平滑化处理
阈值自动检测
  • 使用平滑后的曲线根据凹凸性找拐点
  • 使用大津法(Otsu)计算最大类间方差判断阈值
canny算子
  • canny算子是一个很强大的库,它有很多很强大的功能,用这个库可以很容易解决这道题,也可以用于优化轮廓,增加准确

总的来说,这篇文章中我们所使用的方法仅仅是在一维的角度上对数据进行处理,肯定存在一定的误差值和失败率,canny算子是在二维矩阵的角度上经过一系列的数学公式计算查找阈值,但是如果要借用他来自动查找的话,需要懂得他的原理和源码,所以对于细胞的识别对于我们目前来说,还是任重道远,之后继续加油吧~~~

五、源码展示

from tkinter import *
from PIL import Image, ImageTk
from tkinter.filedialog import askopenfilename
import cv2 as cv
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter


def slopeComparison(i,ls,ls1):
    left = abs(ls[i] - ls[i - 1]) + abs(ls[i-1] - ls[i-2])
    right = ls[i+2] - ls[i]
    mean = np.mean(ls1)
    if right > mean * left and left != 0 and (abs(ls[i] - ls[i - 1]) > 0.01) and 0 < ls[i] < mean:
        return True
    else:
        return False


def choosepic():
    global scaleX1
    global scaleX2
    global cell_o
    global words
    global rec

    path_ = askopenfilename()
    print(path_)
    rec['text'] = "推荐数值:"
    cell_o = cv.imread(path_, 1)
    changePic(cell_o)

    ls = []

    for i in range(0, 256):
        scaleX1.set(i)
        scaleX2.set(100)
        try:
            whole_area, nucleus_area, whole_area_ellipse, nucleus_area_ellipse = runDetect(i, 100)
            ls.append(abs(whole_area - whole_area_ellipse))
        except Exception:
            ls.append(0.1)

    # ls = list(savgol_filter(ls, 9, 5, mode= 'nearest'))
    flag = False
    ls1 = []
    amount = 0
    for i in range(len(ls)-1):
        ls1.append(abs(ls[i+1] - ls[i]))
    for i in list(range(2, len(ls) - 2))[::-1]:
        if slopeComparison(i,ls,ls1):
            if flag == False:
                flag = True
                print(i-1)
                print("{} {:.2f}".format(i-1, ls[i]))
                rec['text'] = rec['text'] + " " + str(i-1)
                amount = amount + 1
                if amount == 1:
                    scaleX1.set(i-1)
                elif amount == 2:
                    scaleX2.set(i-1)
                    runDetect(scaleX1.get(),scaleX2.get())
        else:
            flag = False

    # for i in range(0, len(ls)):
        # print("{} {:.2f}".format(i, ls[i]))
    plt.plot(range(0, 256), ls, linewidth=1, color='r', marker='o', markersize=0)
    plt.show()


def changePic(cell_o):
    dst = cv.cvtColor(cell_o, cv.COLOR_BGR2RGB)
    up_width = 400
    up_height = 400
    up_points = (up_width, up_height)
    dst = cv.resize(dst, up_points, interpolation=cv.INTER_CUBIC)

    current_image = Image.fromarray(dst)

    imgtk = ImageTk.PhotoImage(image=current_image)
    image_label.config(image=imgtk)
    image_label.image = imgtk  # keep a reference


def runDetect(x1, x2):
    global cell_o
    global words
    cell = cv.cvtColor(cell_o, cv.COLOR_BGR2GRAY)

    rows, cols = cell_o.shape[:2]  # 获得图片的高和宽
    mapx = np.zeros(cell_o.shape[:2], np.float32)  # mapx是一个图片大小的变换矩阵
    mapy = np.zeros(cell_o.shape[:2], np.float32)  # mapy也是一个图片大小的变换矩阵
    for i in range(rows):  # 遍历,将获得到图片的宽和稿依次赋予两个变换矩阵
        for j in range(cols):
            mapx.itemset((i, j), j)
            mapy.itemset((i, j), i)
    # img = cv.remap(cell_o, mapx, mapy, cv.INTER_LINEAR)  # 用仿射函数实现仿射变换

    # 将图像转换为二值图,第二个值为阈值
    ret, cell_2v_whole = cv.threshold(cell, x1, 255, cv.THRESH_BINARY_INV)  # 获取整个细胞的二值图
    ret, cell_2v_nucleus = cv.threshold(cell, x2, 255, cv.THRESH_BINARY_INV)  # 获取细胞核的二值图
    cv.namedWindow("cell_whole", 0)
    cv.resizeWindow("cell_whole", 500, 500)
    cv.imshow("cell_whole", cell_2v_whole)
    cv.namedWindow("cell_nucleus", 0)
    cv.resizeWindow("cell_nucleus", 500, 500)
    cv.imshow("cell_nucleus", cell_2v_nucleus)

    # 基于二值图像用外轮廓的模式,通过全点连接轮廓的方法提取轮廓
    contours_whole, hierarchy = cv.findContours(cell_2v_whole, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_NONE)
    contours_nucleus, hierarchy = cv.findContours(cell_2v_nucleus, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_NONE)

    # 获取目标图像的最小矩阵,找到最大的轮廓
    maxArea = 0
    maxContour = {}
    for item in contours_whole:
        x, y, w, h = cv.boundingRect(item)
        if w * h > maxArea:
            maxArea = w * h
            maxContour = item
    whole_area = 0
    whole_area += cv.contourArea(maxContour)

    img = cv.drawContours(cell_o.copy(), maxContour, -1, (0, 255, 0), 1)

    cnt = maxContour
    ellipse = cv.fitEllipse(cnt)
    (x, y), (a, b), angle = cv.fitEllipse(cnt)
    cv.ellipse(img, ellipse, (255, 0, 0), 1)
    whole_area_ellipse = math.pi * a * b / 4

    # 重复,查找细胞核的最大轮廓
    maxArea = 0
    maxContour = {}
    for item in contours_nucleus:
        x, y, w, h = cv.boundingRect(item)
        if w * h > maxArea:
            maxArea = w * h
            maxContour = item
    nucleus_area = 0
    nucleus_area += cv.contourArea(maxContour)
    img = cv.drawContours(img, maxContour, -1, (0, 0, 255), 1)

    cnt = maxContour
    ellipse = cv.fitEllipse(cnt)
    (x, y), (a, b), angle = cv.fitEllipse(cnt)
    cv.ellipse(img, ellipse, (255, 0, 0), 1)
    nucleus_area_ellipse = math.pi * a * b / 4

    # img = cv.rectangle(img,(x,y),(x+w,y+h),(0,255,0),1)
    changePic(img)
    # 绘制目标框

    words['text'] = "细胞质面积:{} 细胞核面积:{} \n 拟合面积:{:.3f} {:.3f} \n 拟合核质面积比:{:.3f},细胞核质面积比:{:.3f}".format(whole_area,
                                                                                                   nucleus_area,
                                                                                                   whole_area_ellipse,
                                                                                                   nucleus_area_ellipse,
                                                                                                   nucleus_area_ellipse / (whole_area_ellipse - nucleus_area_ellipse),
                                                                                                   nucleus_area / (whole_area - nucleus_area)
                                                                                                   )

    return whole_area, nucleus_area, whole_area_ellipse, nucleus_area_ellipse


def handleLeftButtonReleaseEvent(e):
    global scaleX1
    global scaleX2
    runDetect(scaleX1.get(), scaleX2.get())


root = Tk()
root.geometry('500x600')
root.title('细胞胞质比计算工具')

cell_o = cv.imread("../cell.jpg", 1)

path = StringVar()
Button(root, text='选择图片', command=choosepic).pack()
file_entry = Entry(root, state='readonly', text=path)
# file_entry.pack()
image_label = Label(root)
image_label.pack()
scaleX1 = Scale(root, from_=0, to=255, orient=HORIZONTAL, length=255)
scaleX2 = Scale(root, from_=0, to=255, orient=HORIZONTAL, length=255)
scaleX1.bind("<ButtonRelease-1>", handleLeftButtonReleaseEvent)
scaleX2.bind("<ButtonRelease-1>", handleLeftButtonReleaseEvent)
scaleX1.pack()
scaleX2.pack()
words = Label(root)
words.pack()
rec = Label(root)
rec.pack()

# plt.plot(list(ls),list(range(1,256)),linewidth=1,color='r',marker='o',markersize=0)
# plt.xlabel('Position(2-Theta)')
# plt.ylabel('Intensity')
# plt.legend()
# plt.show()

root.mainloop()

  • 10
    点赞
  • 46
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

图南zzz

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值