Reference:《Super-Resolution Through Neighbor Embedding》Hong Chang,Dit-Yan Yeung,Yimin Xiong
本文主旨：利用局部邻域嵌入流行学习方法，提出可邻域嵌入的实例学习超分辨方法，该方法假设低分辨图像块与高分辨图像块组成的两个流行具有相似的局部几何结构。根据这假设，将低分辨输入图像块与训练集中的K近邻之间的重构关系，映射到对应的高分辨图像块上，合成需要的高分辨图像。与以往实例学习方法相比，邻域嵌入方法不需要学习大量的训练样本，在重建质量与计算复杂度方面能保持一定的平衡。
基于邻域嵌入的超分辨率重建算法主要分为两大部分：样本训练和图像重建。
（1）样本训练
对自然图像建立图像库，然后以图像块为单位，抽取各图像块的特征信息，生成该图像库的低分辨率样本集及对应的高分辨率样本集；
步骤：
①降质
②特征提取
③分块
（2）图像重建
同样对待重建图像进行分块处理，然后在低分辨率样本集中，为每一个待重建图像块寻找最相似的K 个样本块，同时可以获得这些低分辨率样本块所对应的高分辨率样本块，将K 个低分辨率近邻块的进行线性组合重建当前低分辨率块，使重建误差最小，求得重建权重，以该权重对对应的高分辨率样本块加权求和，即可重建当前图像块，最终可重建整幅图像。
步骤：
①特征提取
②分块
③寻找K个近邻块
④计算权值
⑤重建图像
典型代码附录：
①LLE代码及解析：
LLE算法代码 % LLE ALGORITHM (using K nearest neighbors) % % [Y] = lle(X,K,dmax) % % X = data as D x N matrix (D = dimensionality, N = #points) % K = number of neighbors(领域点的个数) % dmax = max embedding dimensionality(最大嵌入维数) % Y = embedding as dmax x N matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Y] = lle(X,K,d) [D,N] = size(X); %D是矩阵的行数，N是矩阵的列数 fprintf(1,'LLE running on %d points in %d dimensions\n',N,D); % STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS %寻找邻居数据点 fprintf(1,'-->Finding %d nearest neighbours.\n',K); X2 = sum(X.^2,1); %矩阵X中的每个元素以2为指数求幂值，并且竖向相加 %if two point X=(x1,x2),Y=(y1,y2) %than the distance between X and Y is sqtr(x1-y1)+sqtr(x2-y2) distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X; %repmat就是在行方向把X2复制成N份,列方向为1份 [sorted,index] = sort(distance); %sort是对矩阵排序,sorted是返回对每列排序的结果,index是返回排序后矩阵中每个数在矩阵未排序前每列中的位置 neighborhood = index(2:(1+K),:); % STEP2: SOLVE FOR RECONSTRUCTION WEIGHTS %计算重构权 fprintf(1,'-->Solving for reconstruction weights.\n'); if(K>D) fprintf(1,' [note: K>D; regularization will be used]\n'); tol=1e-3; % regularlizer in case constrained fits are ill conditioned else tol=0; end W = zeros(K,N); for ii=1:N z = X(:,neighborhood(:,ii))-repmat(X(:,ii),1,K); % shift ith pt to origin C = z'*z; % local covariance C = C + eye(K,K)*tol*trace(C); % regularlization (K>D) W(:,ii) = C\ones(K,1); % solve Cw=1 W(:,ii) = W(:,ii)/sum(W(:,ii)); % enforce sum(w)=1 end; % STEP 3: COMPUTE EMBEDDING FROM EIGENVECTS OF COST MATRIX M=(I-W)'(I-W) %计算矩阵M=(I-W)'(I-W)的最小d个非零特征值对应的特征向量 fprintf(1,'-->Computing embedding.\n'); % M=eye(N,N); % use a sparse matrix with storage for 4KN nonzero elements M = sparse(1:N,1:N,ones(1,N),N,N,4*K*N); for ii=1:N w = W(:,ii); jj = neighborhood(:,ii); M(ii,jj) = M(ii,jj) - w'; M(jj,ii) = M(jj,ii) - w; M(jj,jj) = M(jj,jj) + w*w'; end; % CALCULATION OF EMBEDDING options.disp = 0; options.isreal = 1; options.issym = 1; [Y,eigenvals] = eigs(M,d+1,0,options); %[Y,eigenvals] = jdqr(M,d+1);%change in using JQDR func Y = Y(:,2:d+1)'*sqrt(N); % bottom evect is [1,1,1,1...] with eval 0 fprintf(1,'Done.\n'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % other possible regularizers for K>D % C = C + tol*diag(diag(C)); % regularlization % C = C + eye(K,K)*tol*trace(C)*K; % regularlization
②训练样本降质
%%% ------ training set high-resolution and low-resolution images -------
RGB_YS = imread('.\pic\train-ab-high.bmp','bmp');
RGB_XS = averaging4(RGB_YS);
RGB_XS = downsize(RGB_XS,4);
function B = averaging4(A);
%average a RGB image A into B
% A: m by n by 3
% B: m by n by 3
[m,n,x] = size(A);
A = double(A);
B = zeros(m,n,3);
for i=1:4:m
for j=1:4:n
temp = round((sum(A(i,j:j+3,:))+sum(A(i+1,j:j+3,:))+sum(A(i+2,j:j+3,:))+sum(A(i+3,j:j+3,:)))/16);
for x=0:3
B(i,j+x,:) = temp;
B(i+1,j+x,:) = temp;
B(i+2,j+x,:) = temp;
B(i+3,j+x,:) = temp;
end
end
end
B = uint8(B);
function B = downsize(A,ratio);
%downsize a RGB image A into B
% A: m by n by ratio
% B: m/ratio by n/ratio by ratio
[m,n,x] = size(A);
newm = floor(m/ratio);
newn = floor(n/ratio);
B = zeros(newm,newn,3);
for i=1:newm
for j=1:newn
B(i,j,:) = A(ratio*i-1,ratio*j-1,:);
end
end
B = uint8(B);
③分块
function [B,row,col] = overlapcut(A,s,overlap);
% cut image A into small overlapping patches
% A : m by n matrix representing the Y component of YIQ
% s : size of the small patches
% overlap : overlapping number of pixels in each dimension
% B: size by size by (row*col) array representing the generated pathes
[m,n] = size(A);
step = s-overlap;
row = floor((m-s+step)/step);%等价于row = floor((m-s+s-overlap)/step) = floor((m-overlap)/step)
col = floor((n-s+step)/step);
for i=1:row
for j=1:col
index = (i-1)*col + j;
B(:,:,index) = A((i-1)*step+1:(i-1)*step+s,(j-1)*step+1:(j-1)*step+s);
end
end
④提取一阶二阶梯度特征和亮度特征
function f = findfeature(A,B,row,col,flag);
% construct the feature vecor for some image
% A : s by s by n matrix, representing patches of Y component of some image
% B : ii by jj matrix of Y component of image
% row: # of rows of patches
% col: # of cols of patches
% f : feature vectors
[s,s,n] = size(A); %n is # of patches
[ii,jj] = size(B);
if (row*col~=n)
fprintf('input parameter error!');
return;
end
f = zeros(s*s,n);
for i=1:n
temp = A(:,:,i);
vector = reshape(temp',s*s,1);
vector = vector - mean(vector);
f(1:s*s,i) = vector; %all piexl values (-mean) in the patches
end
if flag == 1
grad = zeros(ii,jj,4);
for i=1:ii %for each piexl in B matrix
for j=1:jj
if (i+1<=ii) & (i-1>=1)
grad(i,j,1) = B(i+1,j)-B(i-1,j);
end
if (j+1<=jj) & (j-1>=1)
grad(i,j,2) = B(i,j+1)-B(i,j-1);
end
if (i+2<=ii) & (i-2>=1)
grad(i,j,3) = B(i+2,j)-2*B(i,j)+B(i-2,j);
end
if (j+2<=jj) & (j-2>=1)
grad(i,j,4) = B(i,j+2)-2*B(i,j)+B(i,j-2);
end
end
end
for p=1:n
i = ceil(p/col); % left top position of pth patch
j = mod(p-1,col)+1; % in grad matrix is (i,j)
add = s*s;
for si=0:s-1
for sj=0:s-1
for num=1:size(grad,3)
add = add+1;
f(add,p) = grad(i+si,j+sj,num);
end
end
end
end
end
⑤NEE方法 寻找K近邻块，计算权值，重建图像
function [YT,U,neighborhood] = naneighbor(XT,XS,YS,K);
% compute weight matrix U by reconstructing nearest appearance neighbors (nan) of XT in XS
% XT : N by T matrix of T features, each colomn representing a N-D feature vector
% XS : N by S matrix of S features, each colomn representing a N-D feature vector
% YS : M by S matrix of S features, each colomn representing a M-D feature vector (M>N)
% K : # of nearest appearance neighbors, deciding by Euclidean distance
%
% YT : M by T matrix of T features, each colomn representing a M-D feature vector
% U : K by T weight matrix
% neighborhood : K by T matrix indicating the neighbors of each feature vector
[N,T] = size(XT);
[N,S] = size(XS);
[M,SS] = size(YS);
% fprintf('Computing weights U of %d nearest appearance neighbors\n',K);
for i = 1:T
temp = XT(:,i);
distance = dist2(temp',XS');
[sorted,index] = sort(distance');
neighborhood(:,i) = index(2:(K+1));
end
tol=1e-4; % regularlizer in case constrained fits are ill conditioned
U = zeros(K,T);
for ii=1:T
z = XS(:,neighborhood(:,ii))-repmat(XT(:,ii),1,K); % shift ith pt to origin
C = z'*z; % local covariance
if trace(C)==0
C = C + eye(K,K)*tol; % regularlization
else
C = C + eye(K,K)*tol*trace(C);
end
U(:,ii) = C\ones(K,1); % solve C*u=1
U(:,ii) = U(:,ii)/sum(U(:,ii)); % enforce sum(u)=1
end;
for ii = 1:T
YT(:,ii) = YS(:,neighborhood(:,ii))*U(:,ii);
end
⑥图像评价指标代码
function measure(origImg,distImg)
%Program for Image / Picture Quality Measures Calculation
%Program Description
%This program is the main entry of the application.
%This program calculates the difference Image/Picture Quality Measures
%If the input image is rgb, convert it to gray image
noOfDim = ndims(origImg);
if(noOfDim == 3)
origImg = rgb2gray(origImg);
end
noOfDim = ndims(distImg);
if(noOfDim == 3)
distImg = rgb2gray(distImg);
end
%Size Validation
origSiz = size(origImg);
distSiz = size(distImg);
sizErr = isequal(origSiz, distSiz);
if(sizErr == 0)
disp('Error: Original Image & Distorted Image should be of same dimensions');
return;
end
disp('done!!');
PSNR = PeakSignaltoNoiseRatio(origImg, distImg);
disp('PSNR = ');
disp(PSNR);
SS= SSIM(origImg, distImg);
disp('SSIM = ');
disp(SS);
MSS = MSSIM(origImg, distImg);
disp('MSSIM = ');
disp(MSS);
PSNR
%Program for Peak Signal to Noise Ratio Calculation
%Author : Athi Narayanan S
%Student, M.E, Embedded Systems,
%K.S.R College of Engineering
%Erode, Tamil Nadu, India.
%http://www.ksrce.ac.in/
function PSNR = PeakSignaltoNoiseRatio(origImg, distImg)
origImg = double(origImg);
distImg = double(distImg);
[M N] = size(origImg);
error = origImg - distImg;
MSE = sum(sum(error .* error)) / (M * N);
if(MSE > 0)
PSNR = 10*log(255*255/MSE) / log(10);
else
PSNR = 99;
end
SSIM
% % 2009.11.11 % Yanxiong Lu % function result = SSIM(x, y) %% default value C1 = 3; C2 = 30; C3 = 15; alfa = 1; beta = 1; gama = 1; %% x = double(x(:)); y = double(y(:)); N = length(x); mux = mean(x); muy = mean(y); detax = std(x); detay = std(y); tmp = cov(x,y); l = (2*mux*muy+C1) / (mux^2+muy^2+C1); c = (2*detax*detay+C2) / (detax^2+detay^2+C2); s = (tmp(1,2)+C3) / (detax*detay+C3); result = l^alfa * c^beta * s^gama; % if alfa=beta=gama=1, and C3/2=C2; you can use follow expression! % result = ( (2*mux*muy+C1)*(2*tmp(1,2)+C2)) / ((mux^2+muy^2+C1)*(detax^2+detay^2+C2));
MSSIM
%
% 2009.11.11
% Yanxiong Lu
%
function result = MSSIM(x, y)
%% default value
% the size of image x and y must be a multiple of block!!
block = 32;
%%
[M, N] = size(x);
x = double(x);
y = double(y);
sum = 0;
for i = 1:M/block
for j = 1:N/block
tmpx = x( (i-1)*block+1:i*block, (j-1)*block+1:j*block );
tmpy = y( (i-1)*block+1:i*block, (j-1)*block+1:j*block );
sum = sum + SSIM(tmpx, tmpy);
end
end
result = sum/(M*N/block/block);
附录二、超分辨常用训练图像样本和测试图像
附录三、LLE相关的有用PPT
附录四、SR-NE文献和代码
附录五、相关毕业论文----《基于学习的视频超分辨率重建算法研究及实现》