基于matlab的传统算法分割

需求：运用形态学方法先分割出细胞膜的结构，然后运用圆形膨胀法扩展整个细胞范围

分割

 Orig=imread('F:\ProgramFiles\MATLAB\R2016a\bin\1G.bmp');
 I=rgb2gray(Orig);%转为灰度图像

原图：

 

自定阈值分割

 FG1=I>125    %大于125的地方为白（1），小于为黑（0）

OtSU阈值分割

Otsu阈值分割，可以借助函数graythresh来自动获取阈值，然后根据这个阈值进行分割；

 M=graythresh(I)%自动获取阈值
 FG2=im2bw(I,M);%根据M的大小来进行分割

 

迭代阈值分割

首先是读取灰度图片，设置一个精度值（这里设的是S），在设置一个初值变量S1等于灰度图最大值和最小值之和的二分之一，再设置两个变量r1和r2，r1是所有大于初值变量S1的元素和，r2是所有小于等于S1的元素和；在设置变量S2等于（r1的平均数+r2的平均数）/2；在进入while循环，经过转换，让S1等于S2，再重复进行前几步，直到精度小于设定的精度S1；这时得到的S2为阈值，再通过im2bw来显示二值图；

 I=im2double(I);
 S=0.1;
 S1=(max(I(:))+min(I(:)))/2;
 r1=find(I>S1);
 r2=find(I<=S1);
 S2=(mean(I(r1))+mean(I(r2)))/2;
 while abs(S2-S1)>S
    S1=S2;
    r1=find(I>S1);
    r2=find(I<=S1);
    S2=(mean(I(r1))+mean(I(r2)))/2;
 end
 FG3=im2bw(I,S2);

 

区域生长法（有缺陷）

1. 鼠标点击选取生长点

2. 3*3 8-邻域生长

3. 生长准则：图像与生长点灰度差小于阈值

4. 生长点压入栈中，直到栈中没有元素时停止

I= imread('F:\Program Files\MATLAB\R2016a\bin\1G.bmp');
 %imshow(I)
 ​
 %I=rgb2gray(Orig)
 ​
 ​
 if isinteger(I)
     I=im2double(I);
 end
 I = rgb2gray(I);
 ​
 figure 
 imshow(I)
 [M,N]=size(I);
 [y,x]=getpts; %单击取点后，按enter结束
 x1=round(x);
 y1=round(y);
 seed=I(x1,y1); %获取中心像素灰度值
 ​
 J=zeros(M,N);
 J(x1,y1)=1;
 ​
 count=1; %待处理点个数
 threshold=0.15;
 while count>0
     count=0;
     for i=1:M %遍历整幅图像
     for j=1:N
         if J(i,j)==1 %点在“栈”内
         if (i-1)>1&(i+1)<M&(j-1)>1&(j+1)<N %3*3邻域在图像范围内
             for u=-1:1 %8-邻域生长
             for v=-1:1
                 if J(i+u,j+v)==0&abs(I(i+u,j+v)-seed)<=threshold
                     J(i+u,j+v)=1;
                     count=count+1;  %记录此次新生长的点个数
                 end
             end
             end
         end
         end
     end
     end
 end
 ​
 subplot(1,2,1),imshow(I);
 title('original image')
 subplot(1,2,2),imshow(J);
 title('segmented image')

 ​
 ​

 

缺点：只会生长鼠标所选位置周围

分水岭分割算法⭐️⭐️⭐️

原博客地址：MATLAB图像分割——标记分水岭分割算法matlab分水岭算法计算机视觉从零学的博客-CSDN博客

1.步骤：

➊读取图像; ❷求取图像的边界，在此基础上可直接应用分水岭分割算法，但效果不佳; ❸对图像的前景和背景进行标记，其中每个对象内部的前景像素值都是相连的，背景里面的每个像素值都不属于任何目标物体; ❹计算分割函数，应用分水岭变换。

2.读取图像并求其边界

rgb = imread('F:\Program Files\MATLAB\R2016a\bin\1G.bmp');%读取原图像
 I = rgb2gray (rgb) ; %转化为灰度图像
 figure; subplot (121) %显示灰度图像
 imshow (I)
 text (732, 501, ' Image courtesy of Corel', ...
 'FontSize', 7, 'HorizontalAlignment' , ' right')
 hy = fspecial('sobel') ;%sobel算子
 hx = hy';
 Iy = imfilter(double(I),hy,'replicate');%滤波求y方向边缘
 Ix = imfilter(double(I),hx,'replicate');%滤波求x方向边缘
 gradmag = sqrt(Ix.^2 + Iy.^2);%求模
 subplot(122) ; imshow (gradmag, []),%显示梯度
 title ('Gradient magnitude (gradmag)')

 

3.分别对前景和背景进行标记

使用形态学重建技术对前景对象进行标记，首先使用imopen函数对图像进行开操作，如图3中左边图像所示，使用半径为20的圆形结构元素，开操作是膨胀和腐蚀操作的结合。另外一种方法是先对图像进行腐蚀，然后对图像进行形态学重建，处理后的图像如图右边图像所示。

 se = strel('disk', 20) ; %圆形结构元素
 Io = imopen(I,se) ;%形态学开操作
 figure; subplot(121);
 imshow(Io),% 显示执行开操作后的图像
 title('形态学开操作')
 Ie = imerode(I, se);%对图像进行腐蚀
 Iobr = imreconstruct(Ie,I) ;%形态学重建
 subplot (122); imshow(Iobr), %显示重建后的图像
 title ('形态学腐蚀重建之后 ')

在开操作之后进行关操作可以去除一些很小的目标，例如苹果上的茎干和一些污点等，如下图中左边图像所示。另外一种方法是先对图像进行腐蚀，然后对图像进行形态学重建，重建后的图像如图右边图像所示，注意在重建之前需要先对图像求反，然后在重建之后再进行一次求反。

Ioc=imclose(Io,se);%形态学关操作
 figure; subplot(121)
 imshow(Ioc),%显示关操作后的图像
 title('关操作 ')
 Iobrd = imdilate(Iobr, se) ; %对图像进行膨胀
 Iobrcbr = imreconstruct(imcomplement (Iobrd),imcomplement(Iobr)) ;%形态学重建
 Iobrcbr = imcomplement(Iobrcbr);%图像求反
 subplot(122); imshow(Iobrcbr),%显示重建求反后的图像
 title ('腐蚀重建之后求反 ')

对比两幅图像l，以重建为基础的开关操作比一般的开关操作（纵向两幅）在去除小的污点时会更有效，并且不会影响这个图像的轮廓

求腐蚀重建之后的极大值：

fgm = imregionalmax(Iobrcbr) ; %局部极大值
 figure; imshow(fgm), %显示重建后局部极大值图像
 title ('Regional maxima of opening-closing by reconstruction (局部极大值) ')

 

4.求原图的极大值：

为了更好地理解这个结果，我们可以在原来图像的基础上，显示局部极大值，对前景图像进行标记

I2=I;
 I2(fgm) = 255;%局部极大值处像素值设置为255
 figure; imshow(I2),%在原图上显示极大值区域
 title (' Regional maxima superimposed on original image (在原图求极大值) ')

 

注意到图像中还有少部分目标物体，未被正确地标记出，如果这些目标物体不能被正确地进行标记，则不能正确地进行分割。并且,少部分前景目标物体已经扩展到边缘，因此应该收缩一下边缘，可以先对图像进行关操作，然后再进行腐蚀来达到这样的效果。   这个过程会产生一些孤立的像素点，可以使用bwareaopen函数来达到这样的效果，这个过程将像素点数量较少的孤立像素点去除，如下图左边图像所示。 ​ 将图像Iobrcbr 使用合适的阈值转化成二值图像，其中淡颜色的值为背景，转化成的二值图像如图右边图像所示

se2 = strel (ones(5,5));%结构元素
 fgm2 = imclose(fgm,se2);%关操作.
 fgm3 = imerode (fgm2, se2) ; %腐蚀
 fgm4 = bwareaopen (fgm3, 20) ;%开操作
 I3=I;
 I3 (fgm4) = 255;%前景处设置为255
 figure; subplot (121)
 imshow (I3) %显示修改后的极大值区域
 title ( 'Modified regional maxima')
 bw = im2bw (Iobrcbr, graythresh (Iobrcbr)) ; %转化为二值图像
 subplot (122); imshow(bw), %显示二值图像
 title ('Thresholded opening-closing by reconstruction')

求蛋白质含量Demo

%imshow(red_gray)bgm_ = imread('F:\Program Files\MATLAB\R2016a\bin\1G.bmp');%读取原图像
 I = rgb2gray (bgm_) ; %转化为灰度图像
 ​
 se = strel('disk', 20) ; %圆形结构元素
 ​
 Io = imopen(I,se) ;%形态学开操作
 ​
 Ie = imerode(I, se);%对图像进行腐蚀
 ​
 Iobr = imreconstruct(Ie,I) ;%形态学重建
 ​
 Ioc=imclose(Io,se);%形态学关操作
 ​
 Iobrd = imdilate(Iobr, se) ; %对图像进行膨胀
 Iobrcbr = imreconstruct(imcomplement (Iobrd),imcomplement(Iobr)) ;%形态学重建
 Iobrcbr = imcomplement(Iobrcbr);%图像求反
 fgm = imregionalmax(Iobrcbr) ; %局部极大值
 %imshow(fgm)效果还行
 ​
 I2=I;
 I2(fgm) = 255;%局部极大值处像素值设置为255
 se2 = strel (ones(5,5));%结构元素
 fgm2 = imclose(fgm,se2);%关操作.
 fgm3 = imerode (fgm2, se2) ; %腐蚀
 fgm4 = bwareaopen (fgm3, 20) ;%开操作
 %imshow(fgm4) 效果一般
 ​
 I3=I;
 I3 (fgm4) = 255;%前景处设置为255 I3最好
 %imshow(I3)
 % bw = im2bw (Iobrcbr, graythresh (Iobrcbr)) ; %转化为二值图像
 bgm=im2bw(I3, 0.05)
 bgm_mask = uint8(bgm.*1);
 imshow(bgm_mask)
 ​
 red=imread('F:\Program Files\MATLAB\R2016a\bin\1R.bmp')
 red_gray= rgb2gray(red)%转灰度图像
 %imshow(red_gray)
 ​
 pm=sum(sum(red_gray.*bgm_mask)')%图像矩阵.*按位置相乘然后计算在图像G中R整幅图像的灰度值（荧光强度）
 pm_total = sum(sum(red_gray)')%计算R整幅图像的荧光强度
 r=pm/pm_total;%求比例

r为比例，pm为图R在图G中的荧光强度，pm_total是R整幅图像的荧光强度

 

膨胀扩展像素

相关函数及原理

1. imdilate 功能：进行膨胀操作 用法： IM2 = imdilate(IM,SE)

   对灰度图像或二值图像IM进行膨胀操作，返回结果图像IM2。

   SE为由strel函数生成的结构元素对象。

2. Strel 功能：构建形态学运算中的结构元素 用法： strel(shape,parameters)

   shape为形状参数，即设置什么样的结构元素；parameters为控制形状参数大小方向的参数。

   常用shape种类：

   'disk' 平面圆盘型

   'rectangle' 矩形

   'line' 线型

   'square' 正方形

3. 用法：这里根据需要自己设定膨胀的大小，膨胀太大会引起失真

   s2=strel('disk',2)%圆盘半价为2的膨胀
    s20=strel('disk',20)%圆盘半价为20的膨胀
    s10=strel('disk',10)%圆盘半价为10的膨胀
    dilated_mask2=imdilate(I3,s2)
    dilated_mask10=imdilate(I3,s10)
    dilated_mask20=imdilate(I3,s20)

4. 结果：

        

 

加入膨胀后求蛋白质含量Demo

bgm_ = imread('F:\Program Files\MATLAB\R2016a\bin\1G.bmp');%读取原图像
 I = rgb2gray (bgm_) ; %转化为灰度图像
 ​
 se = strel('disk', 20) ; %圆形结构元素
 ​
 Io = imopen(I,se) ;%形态学开操作
 ​
 Ie = imerode(I, se);%对图像进行腐蚀
 ​
 Iobr = imreconstruct(Ie,I) ;%形态学重建
 ​
 Ioc=imclose(Io,se);%形态学关操作
 ​
 Iobrd = imdilate(Iobr, se) ; %对图像进行膨胀
 Iobrcbr = imreconstruct(imcomplement (Iobrd),imcomplement(Iobr)) ;%形态学重建
 Iobrcbr = imcomplement(Iobrcbr);%图像求反
 fgm = imregionalmax(Iobrcbr) ; %局部极大值
 %imshow(fgm)效果还行
 ​
 I2=I;
 I2(fgm) = 255;%局部极大值处像素值设置为255
 se2 = strel (ones(5,5));%结构元素
 fgm2 = imclose(fgm,se2);%关操作.
 fgm3 = imerode (fgm2, se2) ; %腐蚀
 fgm4 = bwareaopen (fgm3, 20) ;%开操作
 %imshow(fgm4) 效果一般
 ​
 I3=I;
 I3 (fgm4) = 255;%前景处设置为255 I3效果最好
 ​
 s2=strel('disk',2)
 I3=imdilate(I3,s2)%将细胞膜外扩2个像素
 ​
 bgm=im2bw(I3, 0.05)
 bgm_mask = uint8(bgm.*1);
 imshow(bgm_mask)
 ​
 red=imread('F:\Program Files\MATLAB\R2016a\bin\1R.bmp')
 red_gray= rgb2gray(red)
 ​
 pm=sum(sum(red_gray.*bgm_mask)')
 pm_total = sum(sum(red_gray)')
 r=pm/pm_total;

 

与原始Demo相比，发现r值（比例值上升），表明膨胀操作有效

批量化处理

将所需要处理的.bmp后缀图片放到一个路径下，每次读取两张图片进行蛋白质含量比例测量，最后r值输出每次结果，out数组可看得见所有比对结果

file_list = dir(strcat('*.bmp'));     
 ​
 file_num = length(file_list);
 data_length = file_num/2;
 out = zeros(2, data_length); % 生成一个 2 * data_length 大小的矩阵
 ​
 for i = 1:2:file_num % 初值 步长 终值
     first = file_list(i).name; % 第一张
     second = file_list(i+1).name; % 第二张
     f= imread(strcat(first));
     s= imread(strcat(second));% 读取图片
     I= rgb2gray(f);
    
     s_g=rgb2gray(s);
     imshow(s_g)
     
     se = strel('disk', 20) ; %圆形结构元素
 ​
     Io = imopen(I,se) ;%形态学开操作
     Ie = imerode(I, se);%对图像进行腐蚀
     Iobr = imreconstruct(Ie,I) ;%形态学重建
     Ioc=imclose(Io,se);%形态学关操作
     Iobrd = imdilate(Iobr, se) ; %对图像进行膨胀
     Iobrcbr = imreconstruct(imcomplement (Iobrd),imcomplement(Iobr)) ;%形态学重建
     Iobrcbr = imcomplement(Iobrcbr);%图像求反
     fgm = imregionalmax(Iobrcbr) ; %局部极大值
 ​
     I2=I;
     I2(fgm) = 255;%局部极大值处像素值设置为255
     se2 = strel (ones(5,5));%结构元素
     fgm2 = imclose(fgm,se2);%关操作.
     fgm3 = imerode (fgm2, se2) ; %腐蚀
     fgm4 = bwareaopen (fgm3, 20) ;%开操作
     %imshow(fgm4) 效果一般
     I3=I;
     I3 (fgm4) = 255;%前景处设置为255 I3最好
 ​
     s2=strel('disk',2)
     I3=imdilate(I3,s2)
     bgm=im2bw(I3, 0.05)
     %imshow(bgm)
     bgm_mask = uint8(bgm.*1);%背景
     
     pm=sum(sum(s_g.*bgm_mask)')
     pm_total = sum(sum(s_g)')
     r=pm/pm_total;
     
     k = uint8(i/2)+1
     out(1, k) = r;
 end

 

结果与上面一致，程序运行无异常

背景去噪

博客：matlab——imadjust函数_一千种风的味道的博客-CSDN博客

t1=cputime
bgm_ = imread('F:\Program Files\MATLAB\R2016a\bin\4g.jpg');%读取原图像
I = rgb2gray (bgm_) ; %转化为灰度图像

se = strel('disk', 20) ; %圆形结构元素

Io = imopen(I,se) ;%形态学开操作

Ie = imerode(I, se);%对图像进行腐蚀

Iobr = imreconstruct(Ie,I) ;%形态学重建

Ioc=imclose(Io,se);%形态学关操作

Iobrd = imdilate(Iobr, se) ; %对图像进行膨胀
Iobrcbr = imreconstruct(imcomplement (Iobrd),imcomplement(Iobr)) ;%形态学重建
Iobrcbr = imcomplement(Iobrcbr);%图像求反
fgm = imregionalmax(Iobrcbr) ; %局部极大值
%imshow(fgm)效果还行

I2=I;
I2(fgm) = 255;%局部极大值处像素值设置为255
se2 = strel (ones(5,5));%结构元素
fgm2 = imclose(fgm,se2);%关操作.
fgm3 = imerode (fgm2, se2) ; %腐蚀
fgm4 = bwareaopen (fgm3, 20) ;%开操作
%imshow(fgm4) 效果一般

I3=I;
I3 (fgm4) = 255;%前景处设置为255 I3最好
%imshow(I3)
% bw = im2bw (Iobrcbr, graythresh (Iobrcbr)) ; %转化为二值图像
s2=strel('disk',2)
I3=imdilate(I3,s2)
bgm=im2bw(I3, 0.05)
bgm_mask = uint8(bgm.*1);
%imshow(bgm_mask)

red=imread('F:\Program Files\MATLAB\R2016a\bin\4r.jpg')
red_gray= rgb2gray(red)
%imshow(red_gray)

% pm=sum(sum(red_gray.*I3)')
% pm_total = sum(sum(pm)')
% r=pm/pm_total;

r_img= red_gray.*bgm_mask
subplot(1,3,1),imshow(r_img),title('开始的图像');
r_r_img=r_img
r_r_img=imadjust(r_r_img,[0.05;1],[0;1])%将低于0.03*25像素的值变为0 这里乘的系数0.03控制去噪的程度
r_r_img_=imadjust(r_r_img,[0.09;1],[0;1])%将低于0.09*25像素的值变为0
subplot(1,3,2),imshow(r_r_img),title('去噪的图像');  
subplot(1,3,3),imshow(r_r_img_),title('更去噪的图像'); 
%imshow(r_img)
pm=sum(sum(red_gray.*bgm_mask)')
pm_=sum(sum(r_r_img)')%计算去噪之后的灰度值

pm_total = sum(sum(red_gray)')
r=pm/pm_total;

t2=cputime
t=t2-t1%计算程序运行的时间

 

最终可膨胀可去噪后的求图像荧光强度值DEMo

file_list = dir(strcat('*.jpg'));     % 这里改后缀用于处理不同格式图片

file_num = length(file_list);
data_length = file_num/2;
out = zeros(2, data_length); % 生成一个 2 * data_length 大小的矩阵

for i = 1:2:file_num % 初值 步长 终值
    first = file_list(i).name; % 第一张
	second = file_list(i+1).name; % 第二张
    f= imread(strcat(first));
    s= imread(strcat(second));% 读取图片
    I= rgb2gray(f);
   
    s_g=rgb2gray(s);
    %imshow(s_g)
    
    se = strel('disk', 20) ; %圆形结构元素

    Io = imopen(I,se) ;%形态学开操作
    Ie = imerode(I, se);%对图像进行腐蚀
    Iobr = imreconstruct(Ie,I) ;%形态学重建
    Ioc=imclose(Io,se);%形态学关操作
    Iobrd = imdilate(Iobr, se) ; %对图像进行膨胀
    Iobrcbr = imreconstruct(imcomplement (Iobrd),imcomplement(Iobr)) ;%形态学重建
    Iobrcbr = imcomplement(Iobrcbr);%图像求反
    fgm = imregionalmax(Iobrcbr) ; %局部极大值

    I2=I;
    I2(fgm) = 255;%局部极大值处像素值设置为255
    se2 = strel (ones(5,5));%结构元素
    fgm2 = imclose(fgm,se2);%关操作.
    fgm3 = imerode (fgm2, se2) ; %腐蚀
    fgm4 = bwareaopen (fgm3, 20) ;%开操作
    %imshow(fgm4) 效果一般
    I3=I;
    I3 (fgm4) = 255;%前景处设置为255 I3最好
    %imshow(I3)
    % bw = im2bw (Iobrcbr, graythresh (Iobrcbr)) ; %转化为二值图像
    s2=strel('disk',2)
    I3=imdilate(I3,s2)
    bgm=im2bw(I3, 0.05)
    %imshow(bgm)
    bgm_mask = uint8(bgm.*1);%背景
    
    
    
    pm=sum(sum(s_g.*bgm_mask)')%pm值为图像的荧光强度灰度值
    r_img=s_g.*bgm_mask
    subplot(1,2,1),imshow(r_img),title('去噪前')
    r_r_img=imadjust(r_img,[0.08;1],[0;1])%将低于0.05*255像素的值变为0 这里乘的系数0.05控制去噪的程度
    subplot(1,2,2),imshow(r_r_img),title('去噪后')
    pm_=sum(sum(r_r_img)')%pm_值为去噪后图像的荧光强度灰度值
    pm_total = sum(sum(s_g)')
    r=pm/pm_total;
    
    k = uint8(i/2)+1
    out(1, k) = pm_;%out记录的就是所有的值
end