【代码备份】pocs.m

超分辨率算法代码

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

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:\Users\chen\Desktop\POCS\code\a_',num2str(num),'.jpg'));  
    else  
      frame = imread(strcat('C:\Users\chen\Desktop\POCS\code\a_',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))

View Code