关于图像显著性（MR）matlab代码详解

本代码内容是关于Saliency Detection via Graph-Based Manifold Ranking的算法详解，想要运行此代码还需要一系列的文件，单纯此代码无法演示结果的。

可以在网上搜索完整文件，简单来说就是基本的SLIC小程序，其余的不是必要项。

将数据集直接放入test文件夹中，运行sample或者demo都可以吧，我的是师兄传给我的，跟各位下载的结果应该差别不大。

（为了避免大家浪费时间，这里贴出百度云下载链接，链接：

关注并且私信我发链接

clear all;close all;clc;
addpath('./function/');
%%------------------------设置参数---------------------%%
theta = 0.1;                            % 控制边缘的权重
alpha = 0.99;                          % 控制流行排序成本函数两个项的平衡
spnumber = 200;                    % 超像素的数量
imgRoot = './test/';                 % 测试图像的路径
saldir = './saliencymap/';         % 显著性图像的输出路径
supdir = './superpixels/';         % 超像素标签的文件路径
mkdir(supdir);
mkdir(saldir);
imnames = dir([imgRoot '*' 'jpg']);
 
disp(imnames);
imname = [imgRoot imnames.name];
[input_im,w] = removeframe(imname);            %预处理去掉边框
[m,n,k] = size(input_im);
 
%%----------------------生成超像素--------------------%%
  imname = [imname(1:end-4) '.bmp'];       %SLIC软件仅支持bmp格式的图片
  comm = ['SLICSuperpixelSegmentation' ' ' imname ' ' int2str(20) ' ' int2str(spnumber) ' ' supdir]; %设置启用SLIC.exe前参数
  %<filename> <spatial_proximity_weight> <number_of_superpixels> <path_to_save_results>
  system(comm);   %启用SLIC.exe，将每一个像素按照超像素的区域分类，同一区域的超像素赋相同的值，值与此区域的颜色亮度无关
  spname = [supdir imnames.name(1:end-4)  '.dat'];
   
  %超像素标签矩阵
  fid = fopen(spname,'r');%fid是文件代号（句柄）
  A = fread(fid, m * n, 'uint32');     %fread(fid, N, 'str')  N代表读入元素个数， 'str'是格式类型，将此文件顺序读出来，m*n个像素点逐行扫描
  A = A+1;     %把A变成正整数或逻辑值
  B = reshape(A,[n, m]);%将顺序读出的值改为原图像的二维格式
  
  superpixels = B';%B的转置矩阵，转成m*n,从列开始，在同一超像素区域内的像素赋同一值，值仅作为计数用，值为n，表示第n个超像素
  fclose(fid);
  spnum = max(superpixels(:));  %实际的超像素数目 
   
 %%----------------------设计图形模型--------------------------%%
 %计算特征值 (mean color in lab color space)
 %对于每个超像素
    input_vals = reshape(input_im, m*n, k);%将原图像按行扫描转换为（m*n）*3的矩阵，即m*n行,3列的矩阵，分别为rgb
    %input_im的值为0-1之间，是正常的rgb取值范围，全0表示黑，全1表示白
    rgb_vals = zeros(spnum,1,3);
    inds = cell(spnum,1);
    for i = 1:spnum%从1到spnum，找到superpixels值相同的像素，放到一个cell中
        inds{i} = find(superpixels==i);
        rgb_vals(i,1,:) = mean(input_vals(inds{i},:),1);%input_vals每一行代表一个像素，将inds{i}中的所有对应的input_vals中的第几行，也就是
                                                         %第几个像素，取平均值，rgb_vals(i,1,1)也就是第i个超像素的r值
                                                      
    end 
  lab_vals = colorspace('Lab<-', rgb_vals);            %rgb转换成lab空间
  seg_vals = reshape(lab_vals,spnum,3);                 % 每个超像素点的特征，将三维向量lab_vals变为二维向量，
                                                        %也就是以每一个超像素（一行Lab为特征）为单位，seg_vals（1,1）表示第1个超像素的L值
   
 % 求得边界
 %求邻接矩阵
  adjloop = zeros(spnum,spnum);%邻接矩阵，不相连为0，相连为1，默认自身不相连
  [m1,n1] = size(superpixels);
 for i = 1:m1-1
    for j = 1:n1-1
        if(superpixels(i,j)~=superpixels(i,j+1))%从列开始搜索，临近一个像素的超像素标记值不相等，说明这个像素是两个超像素的边界线上点，即这两个超像素相连
            adjloop(superpixels(i,j),superpixels(i,j+1)) = 1;%superpixels(i,j)的值代表第几个超像素
            adjloop(superpixels(i,j+1),superpixels(i,j)) = 1;
        end;%横方向
        if(superpixels(i,j)~=superpixels(i+1,j))
            adjloop(superpixels(i,j),superpixels(i+1,j)) = 1;
            adjloop(superpixels(i+1,j),superpixels(i,j)) = 1;
        end;%竖方向
        if(superpixels(i,j)~=superpixels(i+1,j+1))
            adjloop(superpixels(i,j),superpixels(i+1,j+1)) = 1;
            adjloop(superpixels(i+1,j+1),superpixels(i,j)) = 1;
        end;%捺方向
        if(superpixels(i+1,j)~=superpixels(i,j+1))
            adjloop(superpixels(i+1,j),superpixels(i,j+1)) = 1;
            adjloop(superpixels(i,j+1),superpixels(i+1,j)) = 1;
        end;%撇方向
    end;
end;   
bd = unique([superpixels(1,:),superpixels(m,:),superpixels(:,1)',superpixels(:,n)']);%边界的超像素的标记值
for i = 1:length(bd)
    for j = i+1:length(bd)
        adjloop(bd(i),bd(j)) = 1;
        adjloop(bd(j),bd(i)) = 1;
    end
end        %将边界超像素赋值1，是将边界超像素相连
    
    edges = [];
   for i = 1:spnum;         %真实超像素总数
        indext = [];
        ind = find(adjloop(i,:)==1);%邻接矩阵行i为1的项，也就是找到所有与第i个超像素相连接的超像素 
        for j = 1:length(ind)
            indj = find(adjloop(ind(j),:)==1);%indj仅仅起到过渡作用，将与第i个超像素相连接的所有像素的连接情况一行一行的传递到indext中
            indext = [indext,indj];
        end
        
        indext = [indext,ind];%indext是第i个超像素连接的超像素的连接情况，ind是第i个超像素的连接情况
        indext = indext((indext>i));%只保留indext中值大于i的项
        indext = unique(indext);%合并重复项，只保留不重复的项，即获取矩阵indext的不同元素构成的向量
                                %其目的是将第i个超像素的两圈节点与i连接，例如1超像素周围一圈再加边界一共52个与1相连，再加上这52个节点的周围节点也与1相连接，一共92个节点与1相连
        if(~isempty(indext))
            ed = ones(length(indext),2);%
            ed(:,2) = i*ed(:,2);
            ed(:,1) = indext;%第一列是列举了与i节点相连接的节点，第二列的值就全为i，因为1*i=i
            edges = [edges;ed];%edges是加上第二圈后的连接边矩阵，一共两列，第一列是加上第二圈后与第i个节点相连接的所有节点（为避免重复运算增加运算速度，列举出第i个节点的所有连接点后，列举第i+1个节点的所有连接点时，去掉第i个节点），第二列是i
        end
   end
  
  % 计算关联矩阵
   valDistances = sqrt(sum((seg_vals(edges(:,1),:)-seg_vals(edges(:,2),:)).^2,2));
   %b=sum(a,dim); a表示矩阵；dim等于1或者2，1表示每一列进行求和，2表示每一行进行求和
   %对于valDistances来说,valDistances(1)=sqrt(([1连接的超像素列的第一个的L]-[1的L])^2+[1连接的超像素列的第一个的a]-[1的a])^2+[1连接的超像素列的第一个的b]-[1的b])^2)
   valDistances = normalize(valDistances); %Normalize to [0,1]
   weights = exp(-valDistances/theta);%theta=0.1
   W=sparse([edges(:,1);edges(:,2)],[edges(:,2);edges(:,1)], ...%将weight和edges合并，更加直观
       [weights;weights],spnum,spnum);  
%    sparse函数用法：
%    例如：
% 0, 0, 0, 0;
% 0, 0, 1, 0;
% 0, 0, 0, 0;
% 0, 1, 0, 2;
% 计算机存储稀疏矩阵可以有两种思路：
% 1.按照存储一个普通矩阵一样存储一个稀疏矩阵，比如上面这个稀疏矩阵中总共十六个元素（三个非零元素），把这些元素全部放入存储空间中。这种存储方式，在matlab就叫做full storage organization。
% 2.只存储非零元素，那么怎么存储呢？
% (4,2) 1
% (2,3) 1
% (4,4) 2
 
   % 最优化关联矩阵 (公式3)   f∗ =  y/(D − αW) .
   dd = sum(W); %sum函数不写第二个参数默认为列求和
   D = sparse(1:spnum,1:spnum,dd); 
   clear dd;      %S = sparse(i,j,s,m,n,nzmax)由向量i,j,s生成一个m*n的含有nzmax个非零元素的稀疏矩阵S;即矩阵A中任何0元素被去除，非零元素及其下标组成矩阵S
   optAff = eye(spnum)/(D-alpha*W);%   eye函数返回单位矩阵,Y = eye(n)：返回n*n单位矩阵

   mz = diag(ones(spnum,1));%将mz设置为方阵，方阵大小为spnum*spnum，内容全为1
   mz = ~mz;     %将A的对角元素设置为0
   optAff = optAff.*mz;%将optAff的对角元素归零，因为超像素自身和自身之间的关联性并不考虑，要考虑的话也应该是无限大，没有数值，此矩阵用来表示各个超像素之间的关联性
   
   %%-----------------------------显著性检测第一阶段--------------------------%%
  % 为每个超像素计算显著性值
  % 作为种子点的上边界
  % top
    Yt = zeros(spnum,1);%指示向量
    bst = unique(superpixels(1,1:n));
    Yt(bst) = 1;%将上边界的超像素赋值为1
    bsalt = optAff*Yt;%（spnum*spnum）*(spnum*1)，表示从1到spnum的所有超像素对于上边界超像素的关联性和
    bsalt = (bsalt-min(bsalt(:)))/(max(bsalt(:))-min(bsalt(:)));       %括号内只有一个冒号表示历遍所有元素，将所有元素放进【0,1】中，即正规化数据
    bsalt = 1-bsalt;                  %补码为显著性度量
   
  % down
    Yd = zeros(spnum,1);
    bsd = unique(superpixels(m,1:n));
    Yd(bsd) = 1;
    bsald = optAff*Yd;       %f*（i） 此向量中的每个元素表示节点与背景种子点的相关性
    bsald = (bsald-min(bsald(:)))/(max(bsald(:))-min(bsald(:))); 
    bsald = 1-bsald; 
    
% right
    Yr = zeros(spnum,1);
    bsr = unique(superpixels(1:m,1));
    Yr(bsr) = 1;
    bsalr = optAff*Yr;
    bsalr = (bsalr-min(bsalr(:)))/(max(bsalr(:))-min(bsalr(:)));
    bsalr = 1-bsalr;
   
% left
    Yl = zeros(spnum,1);
    bsl = unique(superpixels(1:m,n));
    Yl(bsl) = 1;
    bsall = optAff*Yl;
    bsall = (bsall-min(bsall(:)))/(max(bsall(:))-min(bsall(:)));
    bsall = 1-bsall;  
    
% combine
    bsalc = (bsalt.*bsald.*bsall.*bsalr);
    bsalc = (bsalc-min(bsalc(:)))/(max(bsalc(:))-min(bsalc(:)));
   
% 为每个像素分配显著性值
  tmapstage1 = zeros(m,n);
  for i = 1:spnum
      tmapstage1(inds{i}) = bsalc(i);
  end
  tmapstage1 = (tmapstage1-min(tmapstage1(:)))/(max(tmapstage1(:))-min(tmapstage1(:)));%正规化
   mapstage1 = zeros(w(1),w(2));
   mapstage1(w(3):w(4),w(5):w(6)) = tmapstage1;
   mapstage1 = uint8(mapstage1*255); %显著性值的高低俩分配灰度值，越显著值越大
    
  outname = [saldir imnames.name(1:end-4)  '_stage1' '.png'];
  imwrite(mapstage1,outname);
   
  %%----------------------显著性检测第二阶段-------------------------%%
  % 自适应阈值二值化
  th = mean(bsalc);                    %阈值被设置为整个显著图上的平均显著性
  bsalc(bsalc<th)=0;%背景
  bsalc(bsalc>=th)=1;%前景
   
  % 为每个超像素计算显著性值
  fsal = optAff*bsalc;%optAff是所有超像素对所有超像素的关联性，只选取其中的前景部分，即所有前景的显著性值
   
  % 为每个像素分配显著性值
    tmapstage2 = zeros(m,n);
    for i = 1:spnum
        tmapstage2(inds{i}) = fsal(i);  
    end
    tmapstage2 = (tmapstage2-min(tmapstage2(:)))/(max(tmapstage2(:))-min(tmapstage2(:)));
     
    mapstage2 = zeros(w(1),w(2));
    mapstage2(w(3):w(4),w(5):w(6)) = tmapstage2;
    mapstage2 = uint8(mapstage2*255);
    outname = [saldir imnames.name(1:end-4)  '_stage2' '.png'];  
    imwrite(mapstage2,outname);

根据代码内容我有感而发，f∗ =  y/(D − αW)此函数算出所有超像素的关联性后，稍微整理一下就可以当做显著值赋予灰度值，至于论文中提到的两次流形排序，
其实在代码中只有一次用到了此函数计算了所有超像素之间的关联性（127-135）。第一次和第二次都只是分别用了类似于指示向量的向量提取出来而已，本质上没有进行两次流形排序，
而是算了一次总的，而后再分别提取，算是简化了过程吧，至于模型的化简，仍然还是比较复杂的，需要摸索的，本文仅仅对代码详解注释，希望对大家的理解过程有帮助。