用POCS方法对图像进行超分辨率重构matlab源码.zip

% 位移、角度估计后利用pocs法求超精度图像 %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% s{1}=double(imread('tu1.jpg')); [row,colum]=size(s{1}); r2=1:2*row;c2=1:2*colum; r2=r2/2;c2=c2/2;c2=c2'; s{1}=interp2(s{1},r2,c2,'cubic') img1=uint8(s{1}); for k=2:4 s{2}=double(imread(strcat('tu',num2str(k),'.jpg'))); [row,colum]=size(s{2}); r2=1:2*row;c2=1:2*colum; r2=r2/2;c2=c2/2;c2=c2'; s{2}=interp2(s{2},r2,c2,'cubic'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %估计角度 img2=uint8(s{2}); [im1, des1, loc1] = sift(img1); [im2, des2, loc2] = sift(img2); distRatio = 0.75; des2t = des2'; for i = 1:size(des1,1) dotprods = des1(i,:) * des2t; [vals,indx] = sort(acos(dotprods)); if (vals(1) < 0.15) && (vals(1) < distRatio * vals(2)) match1(i) = indx(1); else match1(i) = 0; end end cols1 = size(im1,2); j = 1; for i = 1: size(des1,1) if (match1(i) > 0) %line([loc1(i,2) loc2(match1(i),2)+cols1], ... %[loc1(i,1) loc2(match1(i),1)], 'Color', 'r'); point1(j,:) = loc1(i,:); point2(j,:) = loc2(match1(i),:); j = j + 1; end end num = sum(match1 > 0); if num <= 1 %fprintf('too few matches found.\n', num); return; end distpoint = point1 - point2; for i = 1 : size(point1,1) if distpoint(i,4) >= 3.14159 distpoint(i,4) = distpoint(i,4) - 3.14159 * 2; else if distpoint(i,4) <= -3.14159 distpoint(i,4) = distpoint(i,4) + 3.14159 * 2; end end end rotation = 0; count = 0; for i = 2 : size(point1,1) distdistpoint = distpoint(i,4) - distpoint(i-1,4); if abs(distdistpoint) > 0.5*abs(sum(abs(distpoint(:,4)))/size(distpoint,1)) continue; end if abs(distpoint(i,4)) > 1.5*abs(sum(abs(distpoint(:,4)))/size(distpoint,1)) continue; end rotation = rotation + distpoint(i,4); count = count + 1; end if rotation == 0 rotation = sum(distpoint(:,4))/size(distpoint,1) * 180 / 3.14159; else rotation = rotation / count * 180 / 3.14159; end phi_rad(1)=0;phi_rad(2)=rotation; %角度补偿 img2 = imrotate(img2,-rotation,'bicubic','crop'); s{2}=double(img2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 位移估计 delta_est = estimate_shift(s,8); s{2}=double(imread(strcat('tu',num2str(k),'.jpg'))); [row,colum]=size(s{2}); r2=1:2*row;c2=1:2*colum; r2=r2/2;c2=c2/2;c2=c2'; s{2}=interp2(s{2},r2,c2,'cubic'); %求img2补偿后的坐标 factor=1; ss = size(s{1}); center = (ss+1)/2; for i=1:2 r{i} = [1:factor:factor*ss(1)]'*ones(1,ss(2)); % create matrix with row indices c{i} = ones(ss(1),1)*[1:factor:factor*ss(2)]; % create matrix with column indices r{i} = r{i}-factor*center(1); % shift rows to center around 0 c{i} = c{i}-factor*center(2); % shift columns to center around 0 coord{i} = [c{i}(:) r{i}(:)]*[cos(phi_rad(i)) sin(phi_rad(i)); -sin(phi_rad(i)) cos(phi_rad(i))]; % rotate r{i} = coord{i}(:,2)+factor*center(1)+factor*delta_est(i,1); % shift rows back and shift by delta_est c{i} = coord{i}(:,1)+factor*center(2)+factor*delta_est(i,2); % shift columns back and shift by delta_est r{i}=reshape(r{i},400,400); c{i}=reshape(c{i},400,400); end %%%% [m,n]=size(s{1}); for m1=1:m for n1=1:n N1=r{2}(m1,n1); N2=c{2}(m1,n1); if ( N1>3 & N1<size(s{1},1)-2 & N2>3 & N2<size(s{1},2)-2 ) %%% Find center of the window where the PSF will be applied rN1 = round(N1); rN2 = round(N2); %%% Calculate the effective window windowX = r{1}(rN1-2:rN1+2,rN2-2:rN2+2); windowY = c{1}(rN1-2:rN1+2,rN2-2:rN2+2); %%% Find the value of the gaussian at these points and normalize 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.*s{1}(rN1-2:rN1+2,rN2-2:rN2+2))); %%% Calculate the residual R = s{2}(m1,n1) - Ihat; temp = 0; %%% Calculate new values for the reference frame if (R>1) s{1}(rN1-2:rN1+2,rN2-2:rN2+2) = s{1}(rN1-2:rN1+2,rN2-2:rN2+2) + ... (weights.*(R-1))./sum(sum(weights.^2)); elseif (R<-1) s{1}(rN1-2:rN1+2,rN2-2:rN2+2) = s{1}(rN1-2:rN1+2,rN2-2:rN2+2) + ... (weights.*(R+1))./sum(sum(weights.^2)); end end end end clear match1;clear point1;clear point2;clear num;clear distpoint; end