Python----opencv计算细胞核质比轮廓描绘椭圆拟合GUI手动调整阈值自定义函数自动检测阈值
Posted 图南zzz
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了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=以上是关于Python----opencv计算细胞核质比轮廓描绘椭圆拟合GUI手动调整阈值自定义函数自动检测阈值的主要内容，如果未能解决你的问题，请参考以下文章
Python + OpenCV分析细胞核细胞质的面积比（自动阈值）