开发流程
1. 基于OpenCV库,对含有细胞的图像进行处理。
2. 统计图像中所含细胞数量。
3. 输出最大和最小细胞的面积,周长,方向,和该细胞的中心位置。(以像素为单位)
4. 输出所有细胞的平均面积。
5. 主要的中间步骤处理结果以cvShowImage的形式输出。
算法具体步骤
1. 图像的读入和转化
读入图像后转化为灰度图像
// Load image to src and Show source image
Mat source = imread("cell1.bmp", CV_LOAD_IMAGE_COLOR);
imshow("Source image", source);
// Get gray graph of source image.
Mat gray;
cvtColor(source, gray, CV_RGB2GRAY);
imshow("Gray image", gray);
2. 分离细胞与背景,并将图像二值化
int thresholdn = Otsu(gray);
Mat bin;
threshold(gray, bin, thresholdn, 255, CV_THRESH_BINARY);
在分离细胞与背景图像的时候,最重要的是要找到一个合适的阈值,来进行图像和背景的二值化分离。这里我阅读了一些参考文献,决定采用otsu算法即最大类间方差法,又称大津算法。它是按图像的灰度特性,将图像分成背景和目标2部分。背景和目标之间的类间方差越大,说明构成图像的2部分的差别越大,当部分目标错分为背景或部分背景错分为目标都会导致2部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。
3. 获取细胞外围轮廓,及其面积、周长等。
这里我们使用findContours()函数来实现对细胞轮廓的获取,在此基础上舍弃一些轮廓,留下最终我们需要的,由于图像中有些和细胞颜色接近而面积远小于细胞的杂质,因此我们还需要对面积进行判断并且舍弃部分杂质。
vector<vector<Point>> contours;
vector<Vec4i> hierarchy;
Mat cont = Mat::zeros(bin.size(), CV_8UC1);
Mat ellp = Mat::zeros(bin.size(), CV_8UC1);
findContours(bin, contours, hierarchy, CV_RETR_TREE, CV_CHAIN_APPROX_SIMPLE, Point(0, 0));
//计算平均面积:
for (int i = 0; i < contours.size(); i++)
{
if (hash[i])
{
area[i] = contourArea(contours[i]);
totalarea += area[i];
}
}
int average = totalarea / counter;
//杂质的判断并记录最大最小细胞的编号
int minarea = cont.rows*cont.cols, maxarea = 0, mini, maxi;
for (int i = 0; i < contours.size(); i++)
{
if (hash[i])
if (area[i]<average * 0.3)
{
hash[i] = false;
counter--;
}
else
{
if (area[i]<minarea)
{
minarea = area[i];
mini = i;
}
if (area[i]>maxarea)
{
maxarea = area[i];
maxi = i;
}
}
}
4. 对细胞进行椭圆拟合,输出最大最小细胞的信息
虽然只需输出最大最小细胞的信息,但为了便于观察分析我对所有细胞都进行了椭圆拟合。
if (hash[i])
{
CvPoint center;
CvSize size;
CvBox2D32f* box;
CvPoint* PointArray;
CvPoint2D32f* PointArray2D32f;
PointArray = (CvPoint*)malloc(count*sizeof(CvPoint));
PointArray2D32f = (CvPoint2D32f*)malloc(count*sizeof(CvPoint2D32f));
// Alloc memory for ellipse data.
box = (CvBox2D32f*)malloc(sizeof(CvBox2D32f));
// Convert CvPoint set to CvBox2D32f set.
for (int j = 0; j<count;j++)
{
PointArray2D32f[j].x = (float)contours[i][j].x;
PointArray2D32f[j].y = (float)contours[i][j].y;
}
// Fits ellipse to current contour.
cvFitEllipse(PointArray2D32f, count, box);
// Convert ellipse data from float to integer representation.
center.x = cvRound(box->center.x);
center.y = cvRound(box->center.y);
size.width = cvRound(box->size.width*0.5);
size.height = cvRound(box->size.height*0.5);
drawContours(cont, contours, i, 255);
ellipse(ellp, center, size, box->angle, 0, 360, 255, 1, CV_AA, 0);
if (i == mini)
cout << "最小细胞面积为" << area[i] << ',' << "周长:" << count << ','<< "短轴与X+轴成:" << 180 - box->angle << "度角," << "中心" << '(' << center.x << ',' << center.y << ");"<<endl;
if (i == maxi)
cout << "最大细胞面积为" << area[i] << ',' << "周长:" << count << ','<< "短轴与X+轴成:" << 180 - box->angle << "度角," << "中心" << '(' << center.x << ',' << center.y << ");" << endl;
// Free memory.
free(PointArray);
free(PointArray2D32f);
free(box);
}
算法实现要点
1. Otsu算法
对于图像I(x,y),前景(即目标)和背景的分割阈值记作T,属于前景的像素点数占整幅图像的比
例记为ω0,其平均灰度μ0;背景像素点数占整幅图像的比例为ω1,其平均灰度为μ1。图像的总平均灰度记为μ,类间方差记为g。
假设图像的背景较暗,并且图像的大小为M×N,图像中像素的灰度值小于阈值T的像素个数记作N0,像素灰度大于阈值T的像素个数记作N1,则有:
ω0=N0/ M×N (1)
ω1=N1/ M×N (2)
N0+N1=M×N (3)
ω0+ω1=1 (4)
μ=ω0*μ0+ω1*μ1 (5)
g=ω0(μ0-μ)^2+ω1(μ1-μ)^2 (6)
将式(5)代入式(6),得到等价公式: g=ω0ω1(μ0-μ1)^2 (7)
采用遍历的方法得到使类间方差最大的阈值T,即为所求。
Otsu算法步骤如下:
设图象包含L个灰度级(0,1…,L-1),灰度值为i的的象素点数为Ni ,图象总的象素点数为N=N0+N1+...+N(L-1)。灰度值为i的点的概为:
P(i) = N(i)/N.
门限t将整幅图象分为暗区c1和亮区c2两类,则类间方差σ是t的函数:
σ=a1*a2(u1-u2)^2 (2)式中,aj为类cj的面积与图象总面积之比,a1 =sum(P(i)) i->t, a2 = 1-a1; uj为类cj的均值,u1 = sum(i*P(i))/a1 0->t,
u2 = sum(i*P(i))/a2, t+1->L-1
该法选择最佳门限t 使类间方差最大,即:令Δu=u1-u2,σb =max{a1(t)*a2(t)Δu^2}
具体实现代码如下:
int Otsu(Mat& src)
{
int height = src.rows;
int width = src.cols;
//histogram
float histogram[256] = { 0 };
for (int i = 0; i < src.rows; i++)
for (int j = 0; j < src.cols; j++)
{
histogram[src.at<uchar>(i, j)]++;
}
//normalize histogram
int size = height * width;
for (int i = 0; i < 256; i++)
{
histogram[i] = histogram[i] / size;
}
//average pixel value
float avgValue = 0;
for (int i = 0; i < 256; i++)
{
avgValue += i * histogram[i]; //整幅图像的平均灰度
}
int threshold;
float maxVariance = 0;
float w = 0, u = 0;
for (int i = 0; i < 256; i++)
{
w += histogram[i]; //假设当前灰度i为阈值, 0~i 灰度的像素(假设像素值在此范围的像素叫做前景像素) 所占整幅图像的比例
u += i * histogram[i]; // 灰度i 之前的像素(0~i)的平均灰度值: 前景像素的平均灰度值
float t = avgValue * w - u;
float variance = t * t / (w * (1 - w));
if (variance > maxVariance)
{
maxVariance = variance;
threshold = i;
}
}
return threshold;
}
Otsu算法二值化结果:
2. 细胞轮廓提取与分离
findContours()函数会提取图像中的所有轮廓,效果如下
此种情况下有两个问题:
1. 由于函数实现中将整幅图的边界像素值当做0来处理,因而会有一个包含了部分图像边界在内的极大轮廓,而且该轮廓并不是我们需要的。
2. 由于细胞内部相比于边界,颜色较浅且较为接近背景颜色,因而利用Otsu算法无法将整个细胞化为一个整体,所以在细胞内会存在“空洞”。
解决方案:
首先,在采用findContours函数的时候,制定MODE为CV_RETR_TREE,这是将所有的轮廓根据包含关系,建成一个树状结构,若大的轮廓包含了小的,则小的轮廓为大的轮廓的儿子节点,且该树可以有多层。
而根据测试,当我们将整幅图像的边界以内宽度为3的边框的像素赋为255的时候,就能够实现一个覆盖整幅图的边框,并将穿越图像边界的细胞(可能为细胞的一部分)完全分离出来。
for (int i = 0; i < bin.rows; i++)
{
for (int j = 1; j < 3; j++)
{
bin.at<uchar>(i, j - 1) = 255;
bin.at<uchar>(i, bin.cols - j) = 255;
}
}
for (int i = 0; i < bin.cols; i++)
{
for (int j = 1; j < 3; j++)
{
bin.at<uchar>(j - 1, i) = 255;
bin.at<uchar>(bin.rows - j, i) = 255;
}
}
于是得到如下算法①将外层轮廓删除②所有父节点是最外层轮廓的(即为细胞最外层)保留,其余轮廓删除。
for (int i = 0; i < contours.size(); i++)
if ((hierarchy[i][3] == -1) || (hierarchy[i][3] != 0))
{
hash[i] = false;
counter --;
}
结果如下:
3. 去除杂质
在观察图像特质之后我们可以得到杂质如下特征:①颜色和细胞颜色较为接近,二值化后为黑色。②面积远小于细胞。
方法一:统计学处理
根据分析很容易得出如下结论,细胞的面积基本相差不大,服从N(μ, σ)的正态分布,其中μ为面积均值,σ为标准差。根据三倍标准差检验,即分布在[μ-3σ,μ+3σ]外的数据为不符合要求数据。
存在问题:
1. 没有总体数据
在没有总体的情况下,如果强行算出所有数据的均值和标准差,会发现混入了杂质后,细胞的均值下降,方差大幅度增加,导致了杂质基本都存在于[μ-3σ,μ+3σ]区间内。
2. 无法迭代操作
如果采用将符合要求的数据逐个插入总体中,进行迭代,则涉及到最刚开始的数据必须是符合要求的细胞的数据,否则如果一开始都是杂质,那么以杂质建立总体样本,就会出现细胞变成了不符合要求的数据,并且被踢出。
方法二:
法二较为简单,利用杂质面积远小于细胞的特性,设置一个参数,为杂质和细胞均值的面积比例,可以通过实验得出,当比例小于0.3的时候就可以很好的区分出细胞和杂质了。但是这一种方法的误差依然存在。