matlab图像分割，从基因芯片荧光图像中提取阴性点（弱）和阳性点（强）

一、要求

1. 芯片具有边框，且不同芯片边框有差异
2. 芯片相对于相机的角度不会严格正对，会略有倾斜
3. 由于一次成像，照明会不均匀
4. 由于工艺或环境影响，芯片上可能会有杂物或划痕，芯片可能会有缺失部分

二、图片3、图片4

三、思路

由于前两张图的图像质量较好，易于进行图像二值化分割计算，后面两张图由于噪点和干扰较多，采用循环+条件来筛选阴性点和阳性点。

1、首先要提取出图片的ROI（芯片区域）

由于原图中有较强的干扰，直接对原图进行二值化分割的效果并不理想，所以利用形态学的方法，去除干扰的部分，只保留芯片区域的ROI.

下面展示了具体的实现细节

chip_edge = edge(chip_bw);%寻找二值图像的边界
se1=strel('square',13);%方型结构元素
BW=imclose(chip_edge,se1);%闭运算去除干扰
A1 = bwmorph(BW,'hbreak',10);%将不同连通域断开连接
BW1 = bwareaopen(A1, 30000);%将连通体面积小于30000的部分删除，去除杂散点
se2=strel('disk',5);
A2=imerode(BW1,se2);%将相连的连通体断开
BW2 = bwareaopen(A2, 10000);%将连通体面积小于10000的部分删除，去除杂散点
se3=strel('disk',6);
A3=imdilate(BW2,se3);%膨胀操作恢复原来联通体大小

在提取出有效区域后，要使阴性点和阳性点的差异变明显，需要对图像进行增强，先进行自适应的直方图增强，然后对图像进行gamma矫正，阳性点和阴性点的差异变明显。

2、对阳性点和阴性点进行计数

注意到图像中有漏液的部分，以及部分点过于明亮，采用二值化的方法再计算边缘并不可行，所以用循环加条件判断来筛选阳性点和阴性点。观察到，阳性点是其九邻域内的最大值，阴性点是其九邻域内的最小值，编程中适当扩大范围。

% 开始计算阳性点数目
[m,n] = size(chip_enhance);
chip_pad = padarray(chip_enhance,[1 1],0,'both');
mean_value = mean(chip_enhance(:));
bright_num = 0;
dark_num = 0;
%防止领域内有其他值重复识别
point_map = zeros(m,n);

figure(4)
subplot(1,2,1),imshow(chip_enhance),title('原图')
subplot(1,2,2),imshow(chip_enhance),title('红色：阳性点，蓝色：阴性点')
hold on 
%利用循环筛选阳性点和阴性点
for i = 3:m-3
    for j = 3:n-3
        %寻找阳性点
        if (chip_enhance(i,j)>34000) && (chip_enhance(i,j)<50000)
            if ismember(1,point_map(i-2:i+2,j-2:j+2)) || ismember(-1,point_map(i-2:i+2,j-2:j+2))
                continue;
            end
            temp = reshape(chip_pad(i-1:i+1,j-1:j+1),1,[]);
            temp = sort(temp,'descend');
            if chip_enhance(i,j) > temp(3) 
                if  j > n/2
                   if (mean(chip_pad(i-2:i+2,j-2:j+2))<1.5*mean_value)
                       plot(j, i,'r.')
                       bright_num = bright_num + 1;
                       point_map(i,j) = 1;
                   end
                else
                    plot(j, i,'r.')
                    bright_num = bright_num + 1;
                    point_map(i,j) = 1;
                    
                end
            end
        end
        %寻找阴性点
        if (chip_enhance(i,j)>10000) && (chip_enhance(i,j)<25000) 
            if ismember(-1,point_map(i-2:i+2,j-2:j+2)) || ismember(1,point_map(i-2:i+2,j-2:j+2))
                continue;
            end
            temp = reshape(chip_pad(i-1:i+1,j-1:j+1),1,[]);
            temp = sort(temp,'descend');
            if chip_enhance(i,j) > temp(3) && ((temp(1)-temp(9))>4000)
                if (mean(chip_pad(i-1:i+1,j-1:j+1))>11000)
                    plot(j, i,'b.')
                    dark_num = dark_num +1;
                    point_map(i,j) = -1;
                end
            end
        end
                    
    end
end
hold off

可以看到对阳性点和阴性点有一个较好的识别，漏液的部分和杂质没有检测，

简单说明：

if ismember(1,point_map(i-2:i+2,j-2:j2)) || ismember(-1,point_map(i-2:i+2,j-2:j+2)) 
%这个限制条件视为防止在一个九连通域内出现好几个检测点

(temp(1)-temp(9))>4000
 % 这个限制条件是为了防止一些不是芯片内的点被识别，如下图情况

3、完整程序

% 图片4完整程序
clc;clear;
%% 提取ROI
%读取图像
chip = imread('作业4图像4.tif');
figure(1);
subplot(221);imshow(chip);title('原始图像');
%直方图均衡化(显示原图高对比度)，
J = histeq(chip);
subplot(222);imshow(J);title('直方图均衡化后');
%对图像进行二值化
chip_bw1 = imbinarize(chip(1:1040,1:696));%对图像左半部分进行二值化
chip_bw2 = imbinarize(chip(1:1040,697:1392));%对图像右半部分进行二值化
chip_bw = [chip_bw1,chip_bw2];%将二值化图像合成
subplot(223);imshow(chip_bw);title('二值化后图像');

chip_edge = edge(chip_bw);%寻找二值图像的边界
se1=strel('square',13);%方型结构元素
BW=imclose(chip_edge,se1);%闭运算去除干扰
A1 = bwmorph(BW,'hbreak',10);%将不同连通域断开连接
BW1 = bwareaopen(A1, 30000);%将连通体面积小于30000的部分删除，去除杂散点
se2=strel('disk',5);
A2=imerode(BW1,se2);%将相连的连通体断开
BW2 = bwareaopen(A2, 10000);%将连通体面积小于10000的部分删除，去除杂散点
se3=strel('disk',6);
A3=imdilate(BW2,se3);%膨胀操作恢复原来联通体大小
subplot(224),imshow(A3);title('提取出的芯片区域');
%将芯片区域用矩形画出并显示
st = regionprops(A3, 'BoundingBox', 'Area' );
[maxArea, indexOfMax] = max([st.Area]);
rectangle('Position',.....
    [st(indexOfMax).BoundingBox(1),st(indexOfMax).BoundingBox(2),st(indexOfMax).BoundingBox(3),st(indexOfMax).BoundingBox(4)],......
    'EdgeColor','r','LineWidth',2)

figure(2)
subplot(221),imshow(chip);title('原图提升对比度后，芯片所在位置');

rectangle('Position',.....
    [st(indexOfMax).BoundingBox(1),st(indexOfMax).BoundingBox(2),st(indexOfMax).BoundingBox(3),st(indexOfMax).BoundingBox(4)],......
    'EdgeColor','r','LineWidth',2)

crop_chip = imcrop(chip, [st.BoundingBox]);

se4=strel('disk',20);
A4=imdilate(BW2,se4);%膨胀操作恢复原来联通体大小
crop_A4 = imcrop(A4, [st.BoundingBox]);

subplot(222),imshow(crop_chip);title('原图ROI裁剪');

chip_enhance = adapthisteq(crop_chip);
subplot(223),imshow(chip_enhance);title('裁剪ROI增强对比度');



chip_enhance = imadjust(chip_enhance,[],[],0.7);

subplot(224);imshow(chip_enhance);title('增强对比度后Gama矫正');
%% 寻找阳性点、阴性点
% 开始计算阳性点数目
[m,n] = size(chip_enhance);
chip_pad = padarray(chip_enhance,[1 1],0,'both');
mean_value = mean(chip_enhance(:));
bright_num = 0;
dark_num = 0;
%防止领域内有其他值重复识别
point_map = zeros(m,n);

figure(4)
subplot(1,2,1),imshow(chip_enhance),title('原图')
subplot(1,2,2),imshow(chip_enhance),title('红色：阳性点，蓝色：阴性点')
hold on 
%利用循环筛选阳性点和阴性点
for i = 3:m-3
    for j = 3:n-3
        %寻找阳性点
        if (chip_enhance(i,j)>34000) && (chip_enhance(i,j)<50000)
            if ismember(1,point_map(i-2:i+2,j-2:j+2)) || ismember(-1,point_map(i-2:i+2,j-2:j+2))
                continue;
            end
            temp = reshape(chip_pad(i-1:i+1,j-1:j+1),1,[]);
            temp = sort(temp,'descend');
            if chip_enhance(i,j) > temp(3) 
                if  j > n/2
                   if (mean(chip_pad(i-2:i+2,j-2:j+2))<1.5*mean_value)
                       plot(j, i,'r.')
                       bright_num = bright_num + 1;
                       point_map(i,j) = 1;
                   end
                else
                    plot(j, i,'r.')
                    bright_num = bright_num + 1;
                    point_map(i,j) = 1;
                    
                end
            end
        end
        %寻找阴性点
        if (chip_enhance(i,j)>10000) && (chip_enhance(i,j)<25000) 
            if ismember(-1,point_map(i-2:i+2,j-2:j+2)) || ismember(1,point_map(i-2:i+2,j-2:j+2))
                continue;
            end
            temp = reshape(chip_pad(i-1:i+1,j-1:j+1),1,[]);
            temp = sort(temp,'descend');
            if chip_enhance(i,j) > temp(3) && ((temp(1)-temp(9))>4000)
                if (mean(chip_pad(i-1:i+1,j-1:j+1))>11000)
                    plot(j, i,'b.')
                    dark_num = dark_num +1;
                    point_map(i,j) = -1;
                end
            end
        end
                    
    end
end
hold off

str1=['阳性点个数=' num2str(bright_num)];
disp(str1);
str2=['阴性点个数=' num2str(dark_num)];
disp(str2);

不足之处：少部分阳性点没检测出来，部分芯片外的阴性点也检测到了。

参考
https://blog.csdn.net/wolfcsharp/article/details/85945491