MATLAB实现荧光细胞计数

目录

  1. 算法总体实现过程.......................................1

1.1算法具体实现过程................................................1

1.2算法流程图......................................................1

  1. 算法选择与分析.........................................2

2.1图像读入........................................................2

2.2阈值分割........................................................2

2.2.1直方图法....................................................2

2.2.2最大类间方差法(OTSU法)....................................3

2.2.3基于最大熵的阈值分割法......................................3

2.2.4迭代式阈值分割法............................................5

2.3二值化处理......................................................7

2.4膨胀与腐蚀......................................................7

2.4.1膨胀........................................................7

2.4.2腐蚀........................................................7

2.4.3模板选择....................................................8

2.4.4腐蚀膨胀....................................................9

2.5查找连通域.....................................................10

  1. 整体算法..............................................12

1.算法总体实现过程

1.1算法具体实现过程

算法总体分析:本算法首先通过imread函数读入一副荧光细胞图像,然后通过rgb2gray函数将读入的细胞图像进行灰度转换,将一副彩色荧光细胞图像转换为灰度图,再调用graythresh函数用OSTU阈值分割法将目标荧光细胞与背景进行分割,再用im2bw函数将阈值分割后的图像进行二值化,再用strel函数选定合适的核将图像用imdilate函数进行膨胀处理,然后再用strel函数选定合适的核将膨胀后的图像用imerode函数进行腐蚀处理,进行上述处理后再用bwlabel函数将经过以上步骤处理后的图像进行获取连通域,再用regionprops

函数获取连通域的属性以便于后面能够计算出每个连通域的中心坐标,再用find函数查找出每个连通域像素的坐标,然后计算出每个连通域的中心坐标,再用text函数对每个连通域进行标号。

1.2算法流程图

2.算法选择与分析

2.1 图像读入

调用matlab中的imread函数读入荧光细胞图片,再调用matlab中的rgb2gray函数将荧光细胞图片进行灰度转换

代码如下:

IY = imread('D:\OneDrive\\2757da3f6bed1981ffba2795b2ef01d.png');

figure(1);

I = rgb2gray(IY);

imshow(I);title('灰度图');

灰度图转换示例图

2.2阈值分割

阈值分割的方法主要有直方图法,最大类间方差法(OTSU法),基于最大熵的阈值分割法,迭代式阈值分割法四种方法。

2.2.1直方图法:

概念:阈值分割直方图法是一种常用的基于图像直方图的阈值选取方法。它利用图像的灰度直方图来确定合适的阈值,从而实现图像的分割。

基于直方图的阈值分割方法,阈值的选择对于分割结果至关重要。通常有以下几种阈值选取方法:

单一阈值:直接选择灰度直方图的峰值位置作为阈值,适用于图像灰度分布单峰情况。这种方法简单直接,但对于灰度分布不均的图像效果不佳。

多阈值:在灰度直方图中选择多个峰值位置作为阈值,适用于图像灰度分布多峰情况。这种方法可以分割出多个目标,但需要手动选择阈值数量和位置。

基于图像统计信息:使用图像的一些统计信息(如平均值、方差等)作为阈值,适用于特定的应用场景。这种方法对于不同场景的图像需要根据具体情况选择适当的统计信息作为阈值。

代码:

I=rgb2gray(imread('D:\OneDrive\\2757da3f6bed1981ffba2795b2ef01d.png'));

imhist(I);

R=I>125;%大于125的地方为白(1),小于为黑(0

figure(1);

subplot(1,2,1),imshow(R),title('分割图');

subplot(1,2,2),imshow(I),title('原始图像');

直方图法示例图

2.2.2最大类间方差法(OTSU法)

OTSU(又称大律法)按照图像的灰度特征,将图像分成背景和前景两部分。方差是灰度分布均匀性的一种度量,背景和前景之间的类间方差越大,说明构成图像的两部分的差别越大。当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。因此,当类间方差最大时分割错分概率小。在MATLAB中, graythresh函数使用最大类间方差法获得图像的阈值。

代码:

I=rgb2gray(imread('D:\OneDrive\\2757da3f6bed1981ffba2795b2ef01d.png'));

M=graythresh(I);%自动获取阈值

J=im2bw(I,M);%根据M的大小来进行分割

subplot(121),imshow(I);title('原始图像')

subplot(122),imshow(J);title('OTSU阈值分割')

最大类间方差法示例图

2.2.3基于最大熵的阈值分割法

最大熵阈值分割是一种常用的图像分割方法,通过最大化图像熵来确定图像的阈值,从而将图像分为两个部分:背景和前景。其基本思想是找到一个阈值,将图像分为两个部分,使得两个部分的熵之和最大,即最大化图像的熵。

 在图像处理中,熵通常用于衡量图像的复杂度,也常被用作图像分割算法中的评价指标。对于一个二维灰度图像,其熵的计算通常分为以下步骤:

统计图像的灰度直方图,即对每个灰度值出现的次数进行计数。通常采用256个灰度级(0~255)进行离散化。

根据灰度直方图计算每个灰度级出现的概率,即将每个灰度值的出现次数除以总像素数。可以得到一个概率分布向量P,其中P[i]表示灰度值为i的像素出现的概率。

使用熵的定义式计算图像的熵,即:

H = - Σ P[i] * log2(P[i])

其中,Σ是对所有可能的灰度值i求和。

如果需要计算局部熵,则可以将图像划分为若干个小区域,对每个区域计算局部熵。通常将区域大小设置为一个固定值。

在计算熵时,如果像素值分布较为均匀,熵就会比较大,说明图像比较复杂;如果像素值分布较为集中,熵就会比较小,说明图像比较简单。因此,熵可以作为图像复杂度的一种度量方式。同时,熵也可以用于图像分割算法中的阈值选择,通常采用最大熵原理来寻找最佳分割阈值。图像中灰度大于此阈值的像素作为前景,否则作为背景。

代码:

Imag = imread('D:\OneDrive\桌面\2757da3f6bed1981ffba2795b2ef01d.png');

[X, Y] = size(Imag);

imhist(Imag); % 计算图像直方图

hist = imhist(Imag);

p = hist/(X*Y); % 各灰度概率

sumP = cumsum(p);

sumQ = 1-sumP;%256个灰度作为256个分割阈值,分别计算各阈值下的概率密度函数

c0 = zeros(256,256);

c1 = zeros(256,256);

for i = 1:256

                 for j = 1:i

                      if sumP(i) > 0

                       c0(i,j) = p(j)/sumP(i); %计算各个阈值下的前景概率密度函数

        else

            c0(i,j) = 0;

        end

        for k = i+1:256

            if sumQ(i) > 0;

                c1(i,k) = p(k)/sumQ(i); %计算各个阈值下的背景概率密度函数

            else

                c1(i,k) = 0;

                       end

                    end

                 end

end

%计算各个阈值下的前景和背景像素的累计熵

H0 = zeros(256,256);

H1 = zeros(256,256);

for i = 1:256

           for j = 1:i

                if c0(i,j) ~=0

                 H0(i,j) =  - c0(i,j).*log10(c0(i,j));  %计算各个阈值下的前景熵

                     end

                   for k = i+1:256

                  if c1(i,k) ~=0

               H1(i,k) =  -c1(i,k).*log10(c1(i,k));  %计算各个阈值下的背景熵

                  end

                end

             end 

end

HH0 = sum(H0,2);

HH1 = sum(H1,2);

H = HH0 + HH1;

[value, Threshold] = max(H);

BW = imbinarize(Imag, Threshold/255);

figure (1)

subplot(1,2,1),imshow(Imag);title('原始图像')

subplot(1,2,2),imshow(BW);title('最大熵')

基于最大熵的阈值分割法示例图

2.2.4迭代式阈值分割法

迭代式阈值分割法是一种用于图像分割的方法,它将图像分为背景和前景两个部分。该方法基于灰度值的统计特性,通过迭代的方式自适应地选择一个合适的阈值来将图像分割成两个部分。

迭代式阈值分割法的基本思路如下:

(1).选择一个初始阈值T0和预定义两阈值之差d;

(2).利用阈值T0把给定的图像分割成两组图像,记为R1和R2;

(3).计算R1、R2的均值μ1和μ2;

(4).选择新的阈值T1,T1=(μ1+μ2)/2

(5).重复第2至4步,直到连续两次T的差值小于预设值d为止。

迭代式阈值分割法可以用于处理不同的图像,因为它是自适应的,即它能够根据图像的灰度值分布来选择合适的阈值。此外,该方法也易于实现,并且计算速度较快,因此在实际应用中得到了广泛的应用。

代码:

I=imread('D:\OneDrive\桌面\2757da3f6bed1981ffba2795b2ef01d.png');

I=im2double(I);

d=0.01;        %预定义两阈值之差

T0=(min(min(I))+max(max(I)))/2; %初始阈值

R1=find(I>T0);  %分组R1

R2=find(I<=T0); %分组R2

m1=mean(I(R1)); %平均灰度值

m2=mean(I(R2));

T1=(m1+m2)/2;   %新阈值

while abs(T1-T0)>d

            T0=T1;

                  R1=find(I>T0);

            R2=find(I<=T0);

                  T2=(mean(I(R1))+mean(I(R2)))/2;

end

J=im2bw(I,T1);  %把图像转换为二值图,进行阈值分割

set(0,'defaultFigurePosition',[100,100,1000,500]);

set(0,'defaultFigureColor',[1 1 1]);

figure(1);

subplot(1,2,1),imshow(I),title('原图')     

subplot(1,2,2),imshow(J),title('迭代式阈值分割法')

迭代式阈值分割法示例图

根据以上代码以及结果图分析,选择的是OTSU阈值分割方法,原因主要有两个,一是OTSU阈值分割能够自动选择合适的阈值将图像中的目标荧光细胞图像与背景进行分割,而其余的三种方法需要自己设置阈值,如果阈值选择不当就会使目标荧光细胞与背景分割时使某些荧光细胞被当成背景被分割出去,造成连通域减少,在查找连通域时有的连通域不能被查找到,造成荧光细胞计数的误差,并且通过以上结果图也能看出OTSU阈值分割法较其他三种方法效果更好。二是OTSU法代码简洁便于实现,其他三种方法代码繁多并且效果不好。

2.3二值化处理

图像二值化就是将图像上的像素点的灰度值设置为0或255,也就是将整个图像呈现出明显的黑白效果的过程。图像的二值化使图像中数据量大为减少,从而能凸显出目标的轮廓。要得到二值化图像,首先要把图像灰度化,然后将256个亮度等级的灰度图像通过适当的阈值选取而获得仍然可以反映图像整体和局部特征的二值化图像。所有灰度大于或等于阈值的像素被判定为属于特定物体,其灰度值为255,否则这些像素点被排除在物体区域以外,灰度值为0,表示背景或者例外的物体区域。

将上述经过OTSU处理后的图像用im2bw函数转化为二值图像

代码:

I=imread('D:\OneDrive\桌面\2757da3f6bed1981ffba2795b2ef01d.png');

M=graythresh(I);%自动获取阈值

J=im2bw(I,M);%根据M的大小来进行分割

subplot(121),imshow(I);title('原始图像')

subplot(122),imshow(J);title('二值化后图像')

2.4膨胀与腐蚀

2.4.1膨胀:

膨胀是将目标区域接触的背景点合并到该目标物中,使目标边界向外扩张的处理。膨胀可以用来填补目标区域存在的某些空洞,也可以用来消除包含在目标区域中的小颗粒噪声。膨胀处理是腐蚀处理的对偶。

算法步骤:

(1).将预先设定好形状以及原点位置的结构元素移动到图像中可包容(可包容的意思是说:结构

(2).元素的覆盖范围必须完全落在图像内)的所有像素点

(3).判断该结构元素覆盖范围内的像素值是否存在至少一个为1的点,如果存在,则膨胀后图像

(4).中与结构元素原点相同位置的像素值为1,如果该覆盖范围内所有的像素值为0,则膨胀后图像相同位置上的像素值为0

(5).对原图中的所有像素重复进行步骤2

对于此算法先对图像进行膨胀处理填补目标荧光细胞残缺的部分,防止在查找连通域时将残缺的部分当成是一个连通域,造成计数不准确。

2.4.2腐蚀:

腐蚀是一种消除边界点,使边界向内部收缩的过程。可以用来消除小且无意义的目标物。如果两目标物间有细小的连通,可以选取足够大的结构元素,将细小连通腐蚀掉。

算法的具体步骤:

(1).构建结构元素

(2).扫描原图,找到第一个像素值为1(因为是二值图,非0即1)的点。

(3).将预先设定好形状以及原点位置的结构元素的原点移动到该点

(4).判断该结构元素覆盖范围内的图像的像素值是否全部为1,如果是,则腐蚀后图像相同位置上的像素值为1,如果至少有一个像素的值为0,则腐蚀后图像的相同位置上的像素值为0

(5).对原图中所有像素值为1的点重复进行步骤3和步骤4

2.4.3模板选择:

膨胀:

选用一个尺寸为4*4的方形结构体元素,用于膨胀,选用尺寸为4*4大小的原因是荧光细胞的残缺部分比较小因此尺寸可以选小一些,而此结构的形状没有要求,可以选用方形也可以选用圆形。

腐蚀:

选用一个尺寸为半径为9的圆形形结构体元素,用于腐蚀,选用尺寸大小为半径为9的原因是荧光细胞有的是连在一起的而且连接部分较大,因此尺寸要选大一些,而此结构的形状没有要求,可以选用方形也可以选用方形。

模板为4乘4的方形结构体元素与模板为16乘16的方形结构体元素进行膨胀的对比

代码:

IY = imread('D:\OneDrive\桌面\2757da3f6bed1981ffba2795b2ef01d.png');

I = rgb2gray(IY);

I1=graythresh(I);

bw = im2bw(I,I1);

SE1 = strel('square',4);

I2 = imdilate(bw,SE1);

SE2 = strel('square',16);

I3 = imdilate(bw,SE2);

subplot(1,2,1),imshow(I2),title('模板为4的膨胀处理');

subplot(1,2,2),imshow(I3),title('模板为16的膨胀处理');

膨胀处理示例图

模板为半径为9的圆形体元素与模板为3的方形结构体元素进行腐蚀的对比

代码:

IY = imread('D:\OneDrive\桌面\2757da3f6bed1981ffba2795b2ef01d.png');

I = rgb2gray(IY);

I1=graythresh(I);

bw = im2bw(I,I1);

SE1 = strel('disk',9);

I2 = imerode(bw,SE1);

SE2 = strel('disk',3);

I3 = imerode(bw,SE2);

subplot(1,2,1),imshow(I2),title('模板为9的腐蚀处理');

subplot(1,2,2),imshow(I3),title('模板为3的腐蚀处理');

腐蚀处理示例图

2.4.4腐蚀膨胀

先膨胀再腐蚀与先腐蚀再膨胀进行比较

先腐蚀再膨胀

代码:

I = rgb2gray(IY);

I1=graythresh(I);

bw = im2bw(I,I1);

SE1 = strel('square',4);

I2 = imdilate(bw,SE1);

SE2 = strel('disk',9);

I3 = imerode(I2,SE2);

SE3 = strel('disk',9);

I4 = imdilate(bw,SE3);

SE4 = strel('square',4);

I5 = imerode(I4,SE4);

subplot(1,2,1),imshow(I3),title('先膨胀再腐蚀处理');

subplot(1,2,2),imshow(I5),title('先腐蚀再处理');

膨胀腐蚀示意图

根据以上代码以及结果图分析选择模板选用一个尺寸为4*4的方形结构体元素,用于膨胀,再选用一个尺寸为半径为9的圆形形结构体元素,用于腐蚀。进行先膨胀再腐蚀。

原因如下:

选用一个尺寸为4*4的方形结构体元素,用于膨胀,选用尺寸为4*4大小的原因是荧光细胞的残缺部分比较小因此尺寸可以选小一些,而此结构的形状没有要求,可以选用方形也可以选用圆形。

选用一个尺寸为半径为9的圆形形结构体元素,用于膨胀,选用尺寸大小为半径为9的原因是荧光细胞有的是连在一起的而且连接部分较大,因此尺寸要选大一些,而此结构的形状没有要求,可以选用方形也可以选用方形。

先对图像进行膨胀再进行腐蚀的原因是:如果先将荧光细胞图像的二值图进行腐蚀的话会将不完整的细胞中残缺的部分扩大,再进行膨胀有的荧光细胞残缺部分不能被补全,而且会使某些原本分离的细胞连在一块,这样误差就有两个一个是进行腐蚀后使细胞不完整多出现连通域,另一个是膨胀使细胞连在一起而少了的连通域。这样就会使误差大大增大。

而先对图像进行膨胀再进行腐蚀的话能很大程度上减小连通域的误差。原因是荧光细胞残缺部分十分的小,因此先对荧光细胞选用一个较小的正方形模板对其进行膨胀这样能使荧光细胞图像完整,且不会出现多的连通域,再用一个较大的圆形模板对荧光细胞进行腐蚀,这样能使连接在一起的荧光细胞分割开,不会使连通域变小因此大大减少了连通域误差。

2.5查找连通域

      

[L,num] = bwlabel(I3,4);

对腐蚀处理后的图像I3进行连通区域的标记,返回标记矩阵L和连通域的数目

STATS1=regionprops(L,'Area');

计算标记矩阵L中每个连通域的区域属性,这里只计算区域面积。Regionprops在matlab中是用于获取连通域的区域属性的如面积和长度等属性这里只获取面积属性。

ahe=size(STATS1);

获取刚刚用Regionprops函数获取的连通域面积的尺寸即面积的数量,一个面积元素代表一个数量

figure(5);imshow(IY);

创建一个窗口,并显示原始彩色图像

m1=ahe(1,1);

获取连通域的数量并储存在变量m1中

m=zeros(2,m1);

初始化一个2行m1列的零矩阵m,用于储存每个连通域的中心坐标

for i=1:m1

开始一个新的循环,遍历每个连通域以获得每个连通域的x坐标和y坐标

[p,q]=find(L==i);

用find函数查找图像L中每个连通域中所有像素的坐标

temp=[p,q];   

把找到的坐标合成一个temp矩阵

[x,y]=size(temp);

用size函数获取temp的大小,即连通域的像素点个数

m(1,i)=sum(p)/x;

计算每个连通域中每个像素点X坐标的平均值

m(2,i)=sum(q)/x;

计算每个连通域中每个像素点Y坐标的平均值

End

结束循环

for i=1:m1

开始另一个循环遍历每个连通域,以刚获得的每个连通域的中心坐标为坐标标出每个连通域的编号

  figure(5);

  text(m(2,i),m(1,i),int2str(i),'FontSize',12,'color','red','LineWidth',10)

在原图上用text函数在每个连通域的中心位置进行编号。编号I被转换为字符串,字体大小为12,颜色为红色,线条宽度为10。

End

结束循环

3.整体算法

  

       代码:

IY = imread('D:\OneDrive\桌面\2757da3f6bed1981ffba2795b2ef01d.png');

figure(1);

I = rgb2gray(IY);

imshow(I);title('灰度图');

I1=graythresh(I);

bw = im2bw(I,I1);

figure(2);

imshow(bw),title('二值化后图像');

SE1 = strel('square',4);

I2 = imdilate(bw,SE1);

SE2 = strel('disk',9);

I3 = imerode(I2,SE2);

figure(4);

subplot(1,2,1),imshow(I2),title('膨胀处理');

subplot(1,2,2),imshow(I3),title('腐蚀处理');

[L,num] = bwlabel(I3,4);

STATS1=regionprops(L,'Area');

ahe=size(STATS1);

figure(5);imshow(IY);

m1=ahe(1,1);

m=zeros(2,m1);

for i=1:m1

  % 计算目标区域中心,用于显示编号的位置

  [p,q]=find(L==i);

  temp=[p,q];  

  [x,y]=size(temp);

  m(1,i)=sum(p)/x;

  m(2,i)=sum(q)/x;

end

for i=1:m1

  figure(5);

  text(m(2,i),m(1,i),int2str(i),'FontSize',12,'color','red','LineWidth',10)

end

num

  • 7
    点赞
  • 15
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值