超分辨率算法代码
POCS算法，凸集投影法。
pocs.m，没有调用的代码，没看懂。。只有这个函数。。抱歉。
function y = pocs(s,delta_est,factor) % POCS - reconstruct high resolution image using Projection On Convex Sets % y = pocs(s,delta_est,factor) % reconstruct an image with FACTOR times more pixels in both dimensions % using Papoulis Gerchberg algorithm and using the shift and rotation % information from DELTA_EST and PHI_EST % in: % s: images in cell array (s{1}, s{2},...) % delta_est(i,Dy:Dx) estimated shifts in y and x % factor: gives size of reconstructed image %% ----------------------------------------------------------------------- % SUPERRESOLUTION - Graphical User Interface for Super-Resolution Imaging % Copyright (C) 2005-2007 Laboratory of Audiovisual Communications (LCAV), % Ecole Polytechnique Federale de Lausanne (EPFL), % CH-1015 Lausanne, Switzerland % % This program is free software; you can redistribute it and/or modify it % under the terms of the GNU General Public License as published by the % Free Software Foundation; either version 2 of the License, or (at your % option) any later version. This software is distributed in the hope that % it will be useful, but without any warranty; without even the implied % warranty of merchantability or fitness for a particular purpose. % See the GNU General Public License for more details % (enclosed in the file GPL). % % Latest modifications: August 20, 2006, by Karim Krichane max_iter = 50; temp = upsample(upsample(s{1}, factor)', factor)'; y = zeros(size(temp)); coord = find(temp); y(coord) = temp(coord); for i = 2:length(s) temp = upsample(upsample(s{i}, factor)', factor)'; temp = shift(temp, round(delta_est(i, 2)*factor), round(delta_est(i, 1)*factor)); coord = find(temp); y(coord) = temp(coord); end y_prev=y; E=[]; iter=1; blur =[.25 0 1 0 .25;... 0 1 2 1 0;... 1 2 4 2 1;... 0 1 2 1 0;... .25 0 1 0 .25]; blur = blur / sum(blur(:)); wait_handle = waitbar(0, '重构中...', 'Name', '超分辨率重构'); while iter < max_iter waitbar(min(4*iter/max_iter, 1), wait_handle); y = imfilter(y, blur); for i = length(s):-1:1 temp = upsample(upsample(s{i}, factor)', factor)'; temp = shift(temp, round(delta_est(i, 2)*factor), round(delta_est(i, 1)*factor)); coord = find(temp); y(coord) = temp(coord); end delta= norm(y-y_prev)/norm(y); E=[E; iter delta]; iter = iter+1; if iter>3 if abs(E(iter-3,2)-delta) <1e-4 break end end y_prev=y; % if mod(iter,10)==2 % disp(['iteration ' int2str(E(iter-1,1)) ', error ' num2str(E(iter-1,2))]) % end end close(wait_handle);
【其他】貌似这个里面有，可以试一下，没下载过
凸集投影法（POCS）超分辨重建算法MATLAB实现 https://download.csdn.net/download/styyzxjq2009/2312854
POCS 提供了基于POCS算法的超分辨率图像重建的源程序 联合开发网 - pudn.com http://www.pudn.com/Download/item/id/3028355.html
超分辨率的POCS算法–MATLAB中文论坛 http://www.ilovematlab.cn/thread-135641-1-1.html
POCS.m:
close all clear clc t1=clock; NumberOfFrames =3; k = zeros(1,4); %%% 第一帧低分辨率图像与原图 RefImage = imread('a_0.jpg'); %第一帧LW图像 origin=imread('origin.jpg'); %原图 figure(1); imshow(RefImage) RefImageImage =double(RefImage); %%%差值处理，spline，nearest，linear，cubic [x, y] = meshgrid(1:size(RefImage,2), 1:size(RefImage,1)); [X, Y] = meshgrid(1:2.*size(RefImage,2), 1:2.*size(RefImage,1)); upRefImage = interp2(x,y,double(RefImage),X./2,Y./2,'spline'); upRefImage(isnan(upRefImage)) = 0; upRefImage=wiener2(upRefImage); figure(2); imshow(mat2gray(upRefImage)) imwrite(mat2gray(upRefImage),'RefImage_filter_nearest.jpg') %计算信噪比PSNR c=zeros(); [m,n]=size(origin) for i=1:1:m for j=1:1:n minus(i,j)=(origin(i,j)-upRefImage(i,j))^2; end end summ=sum(sum(minus)); PSNR=10*log10(255^2*m*n/summ) %迭代次数 for iter=1:8, disp(iter); for num = 2:NumberOfFrames, %读入其他帧数图像 if (num < 8); frame = imread(strcat('C:UserschenDesktopPOCScodea_',num2str(num),'.jpg')); else frame = imread(strcat('C:UserschenDesktopPOCScodea_',num2str(num),'.jpg')); end frame = double(frame); %%%计算相对第一帧的位置 k = affine(frame,RefImage); u = k(1).*X + k(2).*Y + 2.*k(3); v = -k(2).*X + k(1).*Y + 2.*k(4); mcX = X + u; mcY = Y + v; for m2 = 1:size(frame,2), for m1 = 1:size(frame,1), n1 = 2*m1; n2 = 2*m2; N2 = mcX(n1,n2); N1 = mcY(n1,n2); if ( N1>3 & N1<size(upRefImage,1)-2 & N2>3 & N2<size(upRefImage,2)-2 ) rN1 = round(N1); rN2 = round(N2); windowX = Y(rN1-2:rN1+2,rN2-2:rN2+2); windowY = X(rN1-2:rN1+2,rN2-2:rN2+2); weights = exp(-((N1-windowX).^2+(N2-windowY).^2)./2); weights = weights./sum(sum(weights)); Ihat = sum(sum(weights.*upRefImage(rN1-2:rN1+2,rN2-2:rN2+2))); R = frame(m1,m2) - Ihat; temp = 0; %%% 计算新值 if (R>1) convertedR=double(R); upRefImage(rN1-2:rN1+2,rN2-2:rN2+2) = upRefImage(rN1-2:rN1+2,rN2-2:rN2+2) + ... (weights.*(convertedR-1))./sum(sum(weights.^2)); elseif (R<-1) convertedR=double(R); upRefImage(rN1-2:rN1+2,rN2-2:rN2+2) = upRefImage(rN1-2:rN1+2,rN2-2:rN2+2) + ... (weights.*(convertedR+1))./sum(sum(weights.^2)); end end end end upRefImage(upRefImage<0) = 0; upRefImage(upRefImage>255) = 255; end end %%%展示图像 %%% imwrite(mat2gray(upRefImage),'SRframe_cubic.jpg'); t2=clock; disp(['程序总运行时间：',num2str(etime(t2,t1))]); figure(3); imshow(mat2gray(upRefImage))
另一种POCS算法，myPOCScode.m:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % POCS Image Reconstruction % ------------------------- % AUTHOR: Stephen Rose, Maher Khoury % DATE: March 1, 1999 % PURPOSE: Generates SR frame using the POCS method % % Notes: % -init.m contains the affine transformation parameters ？？？？？？？ % -Assuming a gaussian PSF % -u,v are affine transformation vectors for (x,y) % -mcX,mcY are transformed coordines in SR frame % % Variables: % -ref = LR reference frame % -upref = HR reference frame % -NumberOfFrames = Number of pixel frames to consider % -frame = LR frame currently being examined % -weights = weights based on Gaussian PSF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Initialization 初始化？？？？ %init; clear; close all clc % NumberOfFrames = 4; k = zeros(1,4); wd=1; dlt=5; % max_iter=1; q=2;%放大倍数 % I=imread('E:SRmmreaddiskframe1.bmp'); % [m n]=size(I); % up_ref=zeros(q.*m,q.*n,17); %逐次选择初始图像 %%% Create the high-resolution reference frame % ref=imread(E:SRmmreaddiskframe1.bmp');%低分辨率参考帧 ref=imread('frame1.bmp'); ref=ref(:,:,1); % ref = ref(1:size(ref,1)./2,1:size(ref,2)./2); ref=double(ref); %ref = ref(1:2:size(ref,1),1:2:size(ref,2)); % figure,imshow(ref,[]); % imwrite(mat2gray(ref),'ref.bmp'); % I0=imread('cameraman.bmp');%读入原始清晰图像（计算mse、psnr时，需要用） % mse=zeros(1,max_iter); % psnr=zeros(1,max_iter); % up_ref=zeros(q.*size(ref,1),q.*size(ref,2),iter_max); % for iter_max=1:max_iter % disp(strcat('最大迭代次数:',num2str(iter_max))); % for dlt=1:6 %%%Interpolate values at inbetween points 插值过程 [x, y] = meshgrid(1:size(ref,2), 1:size(ref,1)); [X, Y] = meshgrid(1:q.*size(ref,2), 1:q.*size(ref,1)); upref = interp2(x,y,ref,X./q,Y./q,'bicubic'); %或者linear,bicubic upref1=upref; upref1(isnan(upref1)) = 0; [m,n]=size(upref); % figure,imshow(upref,[]); % imwrite(mat2gray(upref),'upref0.bmp'); % drawnow; %%% Iterate the entire process 迭代过程 % for iter=1:iter_max % disp(strcat('第',num2str(iter),'次迭代')); %%% Iterate over the frames 逐帧迭代 for num = 2:4 frame = imread(strcat('frame',num2str(num),'.bmp')); frame=frame(:,:,1); frame=double(frame); % frame = frame(1:size(frame,1)./q,1:size(frame,2)./q); %%%Calculate the affine motion parameters for this frame %%%计算该帧的仿射系数（估计图像配准参数） k = affine(frame,ref); u = k(1).*X + k(2).*Y + q.*k(3); v = -k(2).*X + k(1).*Y + q.*k(4); %%% Calculate the coordinates of the motion compensated pixels %%% %计算运动补偿像素的坐标？？？？？ mcX = X + u; mcY = Y + v; % Imin=min(min(frame)); % Imax=max(max(frame)); % Rel=zeros(m,n); % for k=1:m % for j=1:n % Rel(k,j)=0.1*(1-2/(Imax-Imin)*abs(upref(k,j)-(Imax-Imin)/2)); % end % end %%% Loop over entire (low-res) frame 逐像素修正 for m2 = 1:size(frame,2) for m1 = 1:size(frame,1) %%% Get high-resolution coordinates n1 = 2*m1; n2 = 2*m2; %%% Get coordinates of the motion compensated pixel 获取运动补偿像素的坐标 N2 = mcX(n1,n2); N1 = mcY(n1,n2); %%% If not a border pixel 排除边缘像素 if ( N1>=wd+1 & N1<=size(upref,1)-wd & N2>=wd+1 & N2<=size(upref,2)-wd ) %？？？？？？原程序为：N1>wd+1 & N1<size(upref,1)-wd %%% Find center of the window where the PSF will be applied %%% 获取PSF作用范围的中心点 rN1 = round(N1); rN2 = round(N2); %%% Calculate the effective window 计算窗口作用范围 windowX = Y(rN1-wd:rN1+wd,rN2-wd:rN2+wd); windowY = X(rN1-wd:rN1+wd,rN2-wd:rN2+wd); %%% Find the value of the gaussian at these points and normalize %%% 计算PSF并归一化 % weights = exp(-1/wd^2*((N1-windowX).^2+(N2-windowY).^2)./2); %%原代码如下计算weights weights = exp(-((N1-windowX).^2+(N2-windowY).^2)./2); weights = weights./sum(sum(weights)); %%% Calculate the value of the estimate Ihat 计算投影像素的估计值 Ihat = sum(sum(weights.*upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd))); %%% Calculate the residual 计算残差 R(m1,m2) = frame(m1,m2) - Ihat; temp = 0; %%% Calculate new values for the reference frame 修正该点的像素值 if (R(m1,m2)>dlt) upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd) = upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd) + (weights.*(R(m1,m2)-dlt))./sum(sum(weights.^2)); % upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd) = upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd) +Rel(rN1-wd:rN1+wd,rN2-wd:rN2+wd).*(R(m1,m2)-dlt); elseif (R(m1,m2)<-dlt) upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd) = upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd) + (weights.*((R(m1,m2)+dlt))./sum(sum(weights.^2))); % upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd) = upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd) +Rel(rN1-wd:rN1+wd,rN2-wd:rN2+wd).*(R(m1,m2)-dlt); % else % upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd) = upref(rN1-wd:rN1+wd,rN2-wd:rN2+wd) + Rel(rN1-wd:rN1+wd,rN2-wd:rN2+wd).*R(m1,m2); end end end end upref(upref<0) = 0; upref(upref>255) = 255; end %upref=255/max(max(upref))*upref; %%% Display the image %%% % up_ref(:,:,start)=uint8(upref); % imwrite(mat2gray(upref),strcat('upref',num2str(start),'.bmp')); % % % % 计算mse与psnr % mse(1,iter_max)=MSE(I0,up_ref(:,:,iter_max)); % psnr(1,iter_max)=PSNR(I0,up_ref(:,:,iter_max)); % imwrite(upref,'SRgirl.tif'); figure,imshow(upref,[]); figure,imshow(upref1,[]); % imwrite(mat2gray(upref),'upref.bmp'); % g=midfilter(upref,1); % figure,imshow(g) % gg=imread('jichang.bmp'); % figure,imshow(gg); % drawnow; % for i=1:10 % up_ref(:,:,i)=double(up_ref(:,:,i)); % imwrite(mat2gray(up_ref(:,:,i),strcat('up_ref',num2str(i),'.bmp'))); % end
网盘文件：
链接：https://pan.baidu.com/s/1qRNjUa93KXKQrFRmwYyWzw
提取码：cmyr