MATLAB蛋白质双向电泳图谱分析

USM 锐化

// An highlighted block
A = imread('蛋白质局部2.bmp');
%A = imread('porteinswim.bmp');
A = rgb2gray(A);
n = 3;  % 3*3模板
a = [-1 -1 -1;-1 8 -1;-1 -1 -1];  % USM算子  什么算法这里放什么算子就可以了，但是注意算子必须是3*3的
%A(a:b,c:d)表示A矩阵的第a到b行,第c到d列的所有元素
[p,q] = size(A);   % 获取输入图像的行列数，要求：p>n,q>n
x1 = double(A);
x2 = zeros(size(x1));    % 确保未被赋值的元素取0
for i = 2:p-1 %如果算子模板不是3*3的话那么这里的最后的值就需要更改了
    for j = 2:q-1
        c=x1(i-1:i+1,j-1:j+1).*a;  % 取出x1中从(i,j)开始的3行3列元素与模板a相乘
               
        s = sum(sum(c));                 % 求c矩阵(即模板)中各元素之和
                
        x2(i,j)=s; % 将模板各元素的均值赋给模板中心位置的元素
    end
end
mg = uint8(x2);
%subplot(121);
figure
imshow(A);title('original');
%subplot(122);
figure
imshow(mg);title('USM-sharp');

二值化

level = graythresh(mg);
ori_im2_3 = im2bw(mg,level);
figure
imshow(ori_im2_3,[]);
title('二值化效果');

Canny边缘检测

I2=edge(ori_im2_3,'canny');
figure
imshow(I2);
title('canny算子分割结果');

提取每个封闭区间的面积，若大于阈值，则认为是蛋白质

因为这里提出的图像会有一些小的封闭环，可能是细胞核，也有可能是一些杂质，这时候就需要将他们都删除。

B = bwboundaries(I2, 'noholes');
[L, num] = bwlabel(I2);
j = 1;
for i = (min(min(L))+1):max(max(L))
    
    area(i) = sum(sum(L == i));
    %调整这个阈值来控制识别的蛋白大小
    
    %蛋白质局部阈值45
    %蛋白质局部2阈值55
    if area(i)>= 55
        real_protein_point(j) = i ;
        %后面可能应该将j清0
        j = j + 1;
    end
end
%算logic矩阵的行数和列数
[ m , n ] =size(L);

%结果矩阵清零
result = zeros(m,n);
%计数矩阵清零
%cnt = 0;
for i1 = 1:m
    for j1 = 1:n
        for s = 1:length(real_protein_point)%算符合要求蛋白质矩阵的长度
            if L(i1,j1) == real_protein_point(s) 
            result(i1,j1) = 1;
            %cnt = cnt+1;
            end
        end
    end
end

    
[posr2, posc2] = find(result == 1);
%cnt                                      % compter des coins
disp('蛋白质的个数是')
fprintf('%d\n',length(real_protein_point))
figure

imshow(A);
hold on;
plot(posc2,posr2,'r+');
title('识别蛋白质效果')

代码的原理是
首先获取每一个小区域的面积。考虑先用 bwlabel 对标签的不同区域打上标记，默认按照 8 联通。然后可以根据标记像素值求取每个小区域的面积。
这个面积就存储在area这个矩阵中，可以先
获取每一个闭合小区域的轮廓曲线坐标，直接使用 bwboundaries 函数进行处理
那么此时 B 为一个 n*1 的元胞数组，n 为闭合区域的个数，这个个数和area是对应的。

随后就是提取每个超过阈值的封闭环的对应标签，再返回到标签矩阵L中反向寻找坐标，最后再原图上标出即可。

最后的效果是这样的：同时会得出蛋白质的个数

以上部分整体代码如下：

clc,clear,close all;
%USM 锐化
A = imread('蛋白质局部2.bmp');
%A = imread('porteinswim.bmp');
A = rgb2gray(A);
n = 3;  % 3*3模板
a = [-1 -1 -1;-1 8 -1;-1 -1 -1];  % USM算子  什么算法这里放什么算子就可以了，但是注意算子必须是3*3的
%A(a:b,c:d)表示A矩阵的第a到b行,第c到d列的所有元素
[p,q] = size(A);   % 获取输入图像的行列数，要求：p>n,q>n
x1 = double(A);
x2 = zeros(size(x1));    % 确保未被赋值的元素取0
for i = 2:p-1 %如果算子模板不是3*3的话那么这里的最后的值就需要更改了
    for j = 2:q-1
        c=x1(i-1:i+1,j-1:j+1).*a;  % 取出x1中从(i,j)开始的3行3列元素与模板a相乘
               
        s = sum(sum(c));                 % 求c矩阵(即模板)中各元素之和
                
        x2(i,j)=s; % 将模板各元素的均值赋给模板中心位置的元素
    end
end
mg = uint8(x2);
%subplot(121);
figure
imshow(A);title('original');
%subplot(122);
figure
imshow(mg);title('USM-sharp');



%二值化
level = graythresh(mg);
ori_im2_3 = im2bw(mg,level);
figure
imshow(ori_im2_3,[]);
title('二值化效果');



I2=edge(ori_im2_3,'canny');
figure
imshow(I2);
title('canny算子分割结果');


B = bwboundaries(I2, 'noholes');


[L, num] = bwlabel(I2);
%area1 = sum(sum(L == 1));
%area2 = sum(sum(L == 2));
j = 1;
for i = (min(min(L))+1):max(max(L))
    
    area(i) = sum(sum(L == i));
    %调整这个阈值来控制识别的蛋白大小
    
    %蛋白质局部阈值45
    %蛋白质局部2阈值55
    if area(i)>= 55
        real_protein_point(j) = i ;
        %后面可能应该将j清0
        j = j + 1;
    end
end
%算logic矩阵的行数和列数
[ m , n ] =size(L);

%结果矩阵清零
result = zeros(m,n);
%计数矩阵清零
%cnt = 0;
for i1 = 1:m
    for j1 = 1:n
        for s = 1:length(real_protein_point)%算符合要求蛋白质矩阵的长度
            if L(i1,j1) == real_protein_point(s) 
            result(i1,j1) = 1;
            %cnt = cnt+1;
            end
        end
    end
end

    
[posr2, posc2] = find(result == 1);
%cnt                                      % compter des coins
disp('蛋白质的个数是')
fprintf('%d\n',length(real_protein_point))
figure

imshow(A);
hold on;
plot(posc2,posr2,'r+');
title('识别蛋白质效果')

蛋白质匹配

这里可以直接使用matlab的surf函数

/% --------------------------------------
% -------- matlab自带SURF函数示例 --------
% --------------------------------------
clc
close all
clear 
tic;

% -------- 读取图像 --------
proportion = 0.9; %定义缩放比例
I1= imread('mathch1.bmp');
I1=imresize(I1,proportion);
I1=rgb2gray(I1);
I2= imread('match2.jpg');
I2=imresize(I2,proportion);
I2=rgb2gray(I2);

% -------- 寻找特征点 --------
points1 = detectSURFFeatures(I1);
points2 = detectSURFFeatures(I2);

% -------- 显示特征点 --------
figure
imshow(I1);
hold on
plot(points1);
figure
imshow(I2);
hold on
plot(points2);

% -------- 计算描述向量 --------
[f1, vpts1] = extractFeatures(I1, points1);
[f2, vpts2] = extractFeatures(I2, points2);

% -------- 进行匹配 --------
indexPairs = matchFeatures(f1, f2, 'Prenormalized', true);
matched_pts1 = vpts1(indexPairs(:, 1));
matched_pts2 = vpts2(indexPairs(:, 2));

% -------- 显示匹配 --------
figure('name','Surf匹配后的图像'); 
showMatchedFeatures(I1,I2,matched_pts1,matched_pts2,'montage'); 
legend('matched points 1','matched points 2');
toc;