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-1for j = 1:n1-1if(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];endindext = [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=iedges = [edges;ed];%edges是加上第二圈后的连接边矩阵，一共两列，第一列是加上第二圈后与第i个节点相连接的所有节点（为避免重复运算增加运算速度，列举出第i个节点的所有连接点后，列举第i+1个节点的所有连接点时，去掉第i个节点），第二列是iendend% 计算关联矩阵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.1W=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元素被去除，非零元素及其下标组成矩阵SoptAff = eye(spnum)/(D-alpha*W);%   eye函数返回单位矩阵,Y = eye(n)：返回n*n单位矩阵mz = diag(ones(spnum,1));%将mz设置为方阵，方阵大小为spnum*spnum，内容全为1mz = ~mz;     %将A的对角元素设置为0optAff = optAff.*mz;%将optAff的对角元素归零，因为超像素自身和自身之间的关联性并不考虑，要考虑的话也应该是无限大，没有数值，此矩阵用来表示各个超像素之间的关联性%%-----------------------------显著性检测第一阶段--------------------------%%% 为每个超像素计算显著性值% 作为种子点的上边界% topYt = zeros(spnum,1);%指示向量bst = unique(superpixels(1,1:n));Yt(bst) = 1;%将上边界的超像素赋值为1bsalt = optAff*Yt;%（spnum*spnum）*(spnum*1)，表示从1到spnum的所有超像素对于上边界超像素的关联性和bsalt = (bsalt-min(bsalt(:)))/(max(bsalt(:))-min(bsalt(:)));       %括号内只有一个冒号表示历遍所有元素，将所有元素放进【0,1】中，即正规化数据bsalt = 1-bsalt;                  %补码为显著性度量% downYd = 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; % rightYr = 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;% leftYl = 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;  % combinebsalc = (bsalt.*bsald.*bsall.*bsalr);bsalc = (bsalc-min(bsalc(:)))/(max(bsalc(:))-min(bsalc(:)));% 为每个像素分配显著性值tmapstage1 = zeros(m,n);for i = 1:spnumtmapstage1(inds{i}) = bsalc(i);endtmapstage1 = (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:spnumtmapstage2(inds{i}) = fsal(i);  endtmapstage2 = (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);