I. Introduction to iterative step size adaptation

The traditional super-resolution reconstruction algorithm usually uses the gradient descent method to solve, and the step size of iteration is usually determined by experience. And the optimal step sizes of different images are often different. A large step size leads to divergence, while a small step size leads to slow convergence. Based on the realization of regularized super-resolution reconstruction algorithm, the algorithm optimizes the step size selection, deduces the optimal step size of each iteration, and ADAPTS it to improve the convergence of super-resolution algorithm, so as to obtain more accurate reconstruction results in a shorter time. For details, please refer to the corresponding paper: Yingqian Wang, Jungang Yang, Chao Xiao, and Wei An, “Fast convergence strategy for multi-image superresolution via adaptive line search,” IEEE Access, vol. 6, no. 1, pp. 9129-9139.

Two, some source code




clear all
clc

filename = 'Set';
files = dir(fullfile( filename,'*.bmp'));
file_num = 2;       % different number corresponds to defferent test images in 'Set'
reg_term = 1;       %regularization term: 1-BTV, 2-Tikhonov

Image =imread([filename,'\',files(file_num).name]); SZ = size(size(Image)); if (SZ(2)==2) % turn grayscale image to RGB image for qw = 1:3 IMAGE (:,:,qw) = Image; end else IMAGE = Image; End % % Image Degradation D = [1, 1, 2, 1, 1, 3, 3, 2); % Shearing shift Gau = fspecial( 'gaussian'[3 3].1);   % Gaussian bluring kernel
spf = 2;                                              % sampling factor
sigma2 = 1;                                       % variation of noise
LR = ImDegrate(IMAGE,D,Gau,spf,sigma2); % image degradation function

%% Turn RGB to YCbCr, and only SR the Y component
[~, ~, ~, M] = size(LR);
for w = 1:M
    LR(:,:,:,w) = rgb2ycbcr(uint8( squeeze(LR(:,:,:,w))));
end

maxiter = 10;  % maximum number of iteration

y1(:, :, :) =  LR(:,:,1, :); y2(:,:,:) = LR(:,:,2, :); y3(:,:,:) = LR(:,:,3, :); HRitp1 = imresize(y1(:,:,1), spf, 'bicubic');  % bicubic interpolation
HRitp1 = ImWarp(HRitp1, -D(1.1), -D(1.2));  % shift recovering

I1 = Wang_SR(HRitp1, y1, D, Gau, spf, maxiter, reg_term);   %Our proposed SR method

HRitp2 =  imresize(y2(:,:,1), spf, 'bicubic');
HRitp2 = ImWarp(HRitp2, -D(1.1), -D(1.2));
I2 = HRitp2;

HRitp3 =  imresize(y3(:,:,1), spf, 'bicubic');
HRitp3 = ImWarp(HRitp3, -D(1.1), -D(1.2));
I3 = HRitp3;

ImZ(:, :, 1) =  I1;
ImZ(:, :, 2) =  I2;
ImZ(:, :, 3) =  I3;

ImZ = ycbcr2rgb(uint8( ImZ)); % Turn YCbCr to RGB

figure; imshow( uint8( ImZ ) ); title('Wang et al.');
figure; imshow( uint8( IMAGE ) ); title('groundtruth');

%%  Evaluation

If = double(ImZ);   %output image
Is = double(IMAGE); %reference image
[row,col,~]=size(If); 

%RMSE
rmse=0;
for color = 1:3
    Ifc = If(:,:,color);  Isc = Is(:,:,color);
    SSE=sum(sum((Ifc-Isc).^2));
    rmsec=sqrt(SSE/(row*col));
    rmse = rmse+rmsec/3;
end
rmse

%PSNR
psnr=0;
for color = 1:3
    Ifc = If(:,:,color);  Isc = Is(:,:,color);
    maxIs = max(max(Isc));
    minIs = min(min(Isc));
    PSNRc = 10*log10((row*col*(maxIs-minIs)^2)/sum(sum((Ifc-Isc).^2)));
    psnr = psnr+PSNRc/3;
end
psnr

%SSIM
ssim=0;
for color = 1:3
    Ifc = uint8(If(:,:,color));  Isc = uint8(Is(:,:,color));
    ssimc = cal_ssim(Ifc, Isc, 0.0);
    ssim = ssim + ssimc/3;
end
ssim
function ssim = cal_ssim( im1, im2, b_row, b_col)

[h, w, ch] = size( im1 );
ssim = 0;
if (ch == 1)
    ssim = ssim_index ( im1(b_row+1:h-b_row, b_col+1:w-b_col), im2(b_row+1:h-b_row,b_col+1:w-b_col));
else
    for i = 1:ch
        ssim = ssim + ssim_index ( im1(b_row+1:h-b_row, b_col+1:w-b_col, i), im2(b_row+1:h-b_row,b_col+1:w-b_col, i));
    end
    ssim = ssim/3;
end
return

function [mssim, ssim_map] = ssim_index(img1, img2, K, window, L)



if (nargin < 2 || nargin > 5)
   mssim = -Inf;
   ssim_map = -Inf;
   return;
end

if (size(img1) ~= size(img2))
   mssim = -Inf;
   ssim_map = -Inf;
   return;
end

[M N] = size(img1);

if (nargin == 2)
   if ((M < 11) || (N < 11))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   window = fspecial('gaussian'.11.1.5);	%
   K(1) = 0.01;					% default settings
   K(2) = 0.03;					%
   L = 255;                                     %
end

if (nargin == 3)
   if ((M < 11) || (N < 11))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   window = fspecial('gaussian'.11.1.5);
   L = 255;
   if (length(K) == 2)
      if (K(1) < 0 || K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

if (nargin == 4)
   [H W] = size(window);
   if ((H*W) < 4 || (H > M) || (W > N))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   L = 255;
   if (length(K) == 2)
      if (K(1) < 0 || K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end

if (nargin == 5)
   [H W] = size(window);
   if ((H*W) < 4 || (H > M) || (W > N))
	   mssim = -Inf;
	   ssim_map = -Inf;
      return
   end
   if (length(K) == 2)
      if (K(1) < 0 || K(2) < 0)
		   mssim = -Inf;
   		ssim_map = -Inf;
	   	return;
      end
   else
	   mssim = -Inf;
   	ssim_map = -Inf;
	   return;
   end
end


img1 = double(img1);
img2 = double(img2);

% automatic downsampling
f = max(1,round(min(M,N)/256));
%downsampling by f
%use a simple low-pass filter 
if(f>1)
    lpf = ones(f,f);
    lpf = lpf/sum(lpf(:));
    img1 = imfilter(img1,lpf,'symmetric'.'same');
    img2 = imfilter(img2,lpf,'symmetric'.'same');

    img1 = img1(1:f:end,1:f:end);
    img2 = img2(1:f:end,1:f:end);
end

C1 = (K(1)*L)^2;
C2 = (K(2)*L)^2;
window = window/sum(sum(window));

mu1   = filter2(window, img1, 'valid');
mu2   = filter2(window, img2, 'valid');
mu1_sq = mu1.*mu1;
mu2_sq = mu2.*mu2;
mu1_mu2 = mu1.*mu2;
sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;
sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;
sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;

if (C1 > 0 && C2 > 0)
   ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
else
   numerator1 = 2*mu1_mu2 + C1;
   numerator2 = 2*sigma12 + C2;
	denominator1 = mu1_sq + mu2_sq + C1;
   denominator2 = sigma1_sq + sigma2_sq + C2;
   ssim_map = ones(size(mu1));
   index = (denominator1.*denominator2 > 0);
   ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
   index = (denominator1 ~= 0) & (denominator2 == 0);
   ssim_map(index) = numerator1(index)./denominator1(index);
end

mssim = mean2(ssim_map);

return


Copy the code

3. Operation results

Matlab version and references

1 matlab version 2014A

2 Reference [1] CAI Limei. MATLAB Image Processing — Theory, Algorithm and Case Analysis [M]. Tsinghua University Press, 2020. [2] Yang Dan, ZHAO Haibin, LONG Zhe. Examples of MATLAB Image Processing In detail [M]. Tsinghua University Press, 2013. [3] Zhou Pin. MATLAB Image Processing and Graphical User Interface Design [M]. Tsinghua University Press, 2013. [4] LIU Chenglong. Proficient in MATLAB Image Processing [M]. Tsinghua University Press, 2015.