基于稀疏表示的图像检索实现

基于稀疏表示的图像检索实现

稀疏表示：
稀疏表示就是将原始的信号(本文指的是图像特征)在一个变换基独立且不相关的过完备字典上进行分解，使得原始复杂冗余的信号可以由分解得到的很少的一组“基元”（特征）进行完全或近似的线性表示[68]，称该自然信号（本文指的是自然图像）在“基元组合”变换域中具有稀疏性，而所有的“基元”则组成一个过完备的字典矩阵。稀疏表示的字典获取方式有多种，最简单的获取方式可以采用固定的字典即可，而相对比较好的获取方式可以通过不同样本集训练学习自适应得到，因为通过这样训练学习得到的字典不仅更加适用当前的图像数据集，也会更加灵活。
下面详细介绍最基本的稀疏表示模型：
假设一个信号y∈R^N可以由过完备字典D线性稀疏表示，字典D=[d_1,d_2,…,d_K]∈R^N (N≪K)，是包含K列的过完备字典，x是待求的稀疏系数向量，那么原始信号y可由字典D中有限个原子的线性组合近似表示。
稀疏编码一般都分为两个过程：
（1）基向量的学习过程，也可以称为字典的学习。在这个过程中，通过无监督学习方法，利用大量的训练样本，学习获得一组冗余的基向量，这组基向量通常反映了训练样本中一些带有本质特征的基元，如图像中的边界、角点，实验表明，经稀疏编码方法学习获得的基向量在朝向、位置和频率的选择上与人类视觉皮层的V1区具有很大的相似性，换言之就是，字典的学习过程模拟了人类视觉皮层对信息的处理过程。图像中学习得到的字典与人为设定的模板字典有很大的相似性，在一定程度上反映了图像的基本单元。
（2）稀疏系数的求解过程。这个过程根据不同的约束条件，可以得到不同的稀疏系数，然后用该系数向量表示图像特征。

实现步骤：
（1）首先对所有肝脏病变CT图像的肿瘤区域提取SIFT局部特征，并切成
大小为7×7的块（也可尝试其他大小的块，如5×5等）；
（2）使用稀疏表示的方法对这些块进行学习，这个学习过程分为以下两个阶段：
基向量的学习过程，也就是字典训练过程。通过K-SVD算法对这些块进行学习，训练得到一个自适应的过完备字典；
稀疏系数的求解。利用已经训练好的字典通过OMP算法求解稀疏系数，并用该系数向量表示图像特征；
（3）利用相似性度量方法中的欧氏距离来计算待查询图像与数据集中其他图像之间的距离，其实就是求它们特征向量之间的距离；
（4）对求得的这些距离按照从小到大进行排序，距离越小就代表两幅肝脏病变CT图像的差异越小，按照距离从小到大输出检索出的图像，实现最终的检索功能；
（5）利用查准率-查全率、F_1- measure等图像检索的性能评价指标来评估检索的效果。

实现代码：

%% 清空环境变量
close all;
clear;
clc;
%% 读取数据并选取训练集和测试样本
data = importdata('AllTypesART.mat');
data1 = importdata('AllTypeART.mat');
[m,n,s] = size(data);
[mm,nn,ss] = size(data1);
p = pwd;
addpath(fullfile(p, '/functions')) % K-SVD dictionary training algorithm
addpath(fullfile(p, '/ksvdbox')) % K-SVD dictionary training algorithm
addpath(fullfile(p, '/ompbox')) % Orthogonal Matching Pursuit algorithm

% load('oneLiver_ART322.mat')
conf.dict_sizes = 256; % 字典大小
conf.window = [7 7]; % 图像块大小7x7
conf.border = [1 1];
conf.overlap = 1; % 重叠部分为1

imgs = {};
dess = {};
tic;
for ii = 1:ss
    % imgs{ii} = double(data1(:,:,ii));
    im = double(data1(:,:,ii));
    figure; 
    imshow(im,[]); 
    % 采样 
    msgbox('Please separate tumour samples','tumour Samples','help'); 
    pause; 
    [x,y] = ginput(2); 
    hold on; 
    plot(x,y,'r*'); 
    x = uint16(x); 
    y = uint16(y); 
    %imgss{ii} = zeros(49,49,'double');
    imgs{ii} = double(im(y(1):y(2),x(1):x(2))); 
    [im1, des1, loc1] = sift(imgs{ii}); % 提取sift特征
    savecommand=['save ','F:\test\DES\',num2str(ii-1),' des1;'];
    eval(savecommand);
    savecommand=['save ','F:\test\LOC\',num2str(ii-1),' loc1;'];
    eval(savecommand);
    dess{ii} = des1;
end
 patches = collect(conf, dess, 1, {}); % 切成特征块, {}处表示无需提取特征
 % patches = dess;

% Set KSVD configuration
ksvd_conf.iternum = 20; % 迭代次数
ksvd_conf.memusage = 'normal'; 
ksvd_conf.dictsize = conf.dict_sizes; % 字典大小
ksvd_conf.Tdata = 10; % maximal sparsity: 最大稀疏度
ksvd_conf.samples = size(patches,2);
ksvd_conf.data = double(patches);
tic;
fprintf('Training [%d x %d] dictionary on %d vectors using K-SVD\n', ...
    size(ksvd_conf.data, 1), ksvd_conf.dictsize, size(ksvd_conf.data, 2))
[conf.dict, gamma] = ksvd(ksvd_conf); % k-svd算法
D = conf.dict;
T = ksvd_conf.Tdata;
toc;

%ps = randperm(s,1);
%testData = double(data(:,:,ps)); % 测试样本 待查图像
testData = double(data(:,:,36)); % 测试样本 待查图像
%str = sprintf('Im %d:',ps);
%imshow(testData,[]);
%title(str);
[im_test, des_test, loc_test] = sift(testData);
dess_tsets = {};
[x,y,u] = size(des_test);
for ii = 1:u
    dess_tsets{ii} = des_test;
end
patches1 = collect(conf, dess_tsets, 1, {}); % 切成特征块, {}处表示无需提取特征
patches1 = double(patches1);
gamma_test = omp(D,patches1,D'*D,T);% des1 = D*gamma 求解稀疏系数
% gamma_test_mean = mean(gamma_test,2);
gamma_test_mean = mean(full(gamma_test),2); % 池化 平均池化
%if ps >1
%    pss = [1:ps-1,ps+1:s];
%else
%    pss = [ps+1:s];
%end
pss = [1:35 37:100];
trainDatas = double(data(:,:,pss));
% 训练样本个数
Nfiles = s-1;
% 阈值
t = 0.010;
t2 = 0.8;
fprintf('开始查找...\n');
range = 0.0:0.1:1.0;
rangeNew = 0.0:0.05:1.0;
[a,b,c]    = meshgrid(range);
[a2,b2,c2] = meshgrid(rangeNew); % 坐标轴
Similarity = zeros(Nfiles, 1); % Nfiles训练样本数量 Similarity表示与Nfiles个训练样本的相似度量值
nResults = 0;
for ii = 1:Nfiles
    trainData = double(trainDatas(:,:,ii));
    [im_train, des_train, loc_train] = sift(trainData);
    des_trains = {};
    [x2,y2,u2] = size(des_train);
    for k = 1:u2
        des_trains{k} = des_train;
    end
    patches2 = collect(conf, des_trains, 1, {}); % 切成特征块, {}处表示无需提取特征
    patches2 = double(patches2);
    gamma_train = omp(D,patches2,D'*D,T);% des1 = D*gamma 求解稀疏系数
    %gamma_train_mean = mean(gamma_train,2);
    gamma_train_mean = mean(full(gamma_train),2);
    % 计算欧式距离
    %DIFF = abs(gamma_test_mean-gamma_train_mean) ./ gamma_test_mean; 
    DIFF = sum((gamma_test_mean-gamma_train_mean).^2); % 欧式距离
    % 保持与相关查询图像的距离值大于预先设定的阈值
    %DIFF = DIFF(gamma_test_mean>t); 
    % keep error values which are smaller than 1:
    %DIFF2 = DIFF(DIFF<t2);
    %L2 = length(DIFF2); 
    % 计算相似度
    Similarity(ii) = DIFF;  
    % (interface): plot images with small similarity measures:
    plotThres = 0.50 * 10 / length(DIFF);
    if (Similarity(ii) < plotThres)
        % fprintf('%70s %5.2f %5d %5d\n', files{i}, median(DIFF2),
        % length(DIFF), L2);
        nResults = nResults+1;
        subplot(2,2,1);imshow(testData,[]);title('Query image');        
        im_trainData = trainData;
        subplot(2,2,2);imshow(im_trainData,[]);
        title('A similar image ... Still Searching ...');
        subplot(2,2,3);
        plot(DIFF)
        %if (length(DIFF2)>1)
        %    subplot(2,2,4); plot(DIFF2);
        %    axis([1 length(DIFF2) 0.2 1])
        %end
        drawnow
    end
end
% find the nResult "closest" images:
% Similarity = 1./Similarity;
[Sorted, ISorted] = sort(Similarity); % 相似度从小到大排序 升序(直方图相似距离越小 越相似)
nResults = 11;
NRows = ceil((nResults+1) / 3); % ceil(x)大于x的最小整数
% plot query image:
subplot(NRows,3,1); imshow(testData,[]); title('Query Image');
% ... plot similar images:
for (i=1:nResults)
    im_trainData = double(trainDatas(:,:,ISorted(i)));
    str = sprintf('Im %d: %.3f',i,100*Sorted(i));
    subplot(NRows,3,i+1); imshow(im_trainData,[]);  title(str);
end
toc;
disp(['检索时间: ',num2str(toc)]);
fprintf('Done\n');