POCS算法.zip_POCS_POCS算法_POCS算法原理_pocs算法代码_slowly2bc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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:\SR\mmread\disk\frame1.bmp'); % [m n]=size(I); % up_ref=zeros(q.*m,q.*n,17); %逐次选择初始图像 %%% Create the high-resolution reference frame % ref=imread(E:\SR\mmread\disk\frame1.bmp');%低分辨率参考帧 ref=imread('frame1.bmp'); % ref=imread('lena.jpg'); 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