First commit

blur.m

@@ -0,0 +1,17 @@
+function [a]=blur(Hr)
+[row col]=size(Hr(:,:,1));
+%operator
+n = row*col;
+e = ones(n,1);
+a1 = spdiags([e e e], -1:1, n, n);
+a2 = spdiags([e e e], -row-1:-row+1, n, n);
+a3 = spdiags([e e e],col-1:col+1, n, n);%blur matrix
+a=a1+a2+a3;
+a(1:row:n-row,col:col:n)=0;
+a(row:row:n-row,2*col+1:col:n-col)=0;
+a(row:row:n,col+1:col:n-col+1)=0;
+a(row+1:row:n-row+1,col:col:n)=0;
+a(row:row:n,1:col:n-col)=0;
+a(2:row:n-row,col:col:n-col)=0;
+a=a/9;
+end

blur2.m

@@ -0,0 +1,19 @@
+function [B]=blur2(Hr)
+
+test_mask=(1/16)*[1,2,1;
+           1,4,1;
+           1,2,1];
+
+[row_tmat col_tmat band]=size(Hr);
+Blur_matrix = sparse(row_tmat*col_tmat,row_tmat*col_tmat);
+
+template = single(toeplitz_generator(1,2,3,row_tmat));
+for m = 1 : row_tmat
+    for n = 1 : row_tmat
+        if template(m,n) ~= 0
+            t = template(m,n);
+            Blur_matrix(col_tmat*(m-1)+1:col_tmat*(m), col_tmat*(n-1)+1:col_tmat*(n)) = toeplitz_generator(test_mask(t,1),test_mask(t,2),test_mask(t,3),col_tmat);
+        end
+    end
+B=Blur_matrix;
+end

build.m

@@ -0,0 +1,33 @@
+clear all; close all; clc;
+
+test_mask = single(-1 * ones(3,3));
+test_mask(2,2) = 8;
+test_matrix = single(reshape(1:24,6,4)')
+
+[row_tmask,col_tmask] = size(test_mask);
+[row_tmat,col_tmat] = size(test_matrix);
+
+
+Blur_matrix = single(zeros(row_tmat*col_tmat,row_tmat*col_tmat));
+
+template = single(toeplitz_generator(1,2,3,row_tmat))
+
+for m = 1 : row_tmat
+    for n = 1 : row_tmat
+        if template(m,n) ~= 0
+            t = template(m,n);
+            Blur_matrix(col_tmat*(m-1)+1:col_tmat*(m), col_tmat*(n-1)+1:col_tmat*(n)) = toeplitz_generator(test_mask(t,1),test_mask(t,2),test_mask(t,3),col_tmat);
+        end
+    end
+end
+
+
+%% test scheme
+
+result1 = filter2(test_mask, test_matrix)
+
+test_vector = reshape(test_matrix',row_tmat*col_tmat,1);
+
+result2_vector = Blur_matrix * test_vector;
+
+result2 = reshape(result2_vector,col_tmat,row_tmat)'

compute_psnr.m

@@ -0,0 +1,16 @@
+function psnr=compute_psnr(im1,im2)
+if size(im1, 3) == 3,
+    im1 = rgb2ycbcr(im1);
+    im1 = im1(:, :, 1);
+end
+
+if size(im2, 3) == 3,
+    im2 = rgb2ycbcr(im2);
+    im2 = im2(:, :, 1);
+end
+
+imdff = double(im1) - double(im2);
+imdff = imdff(:);
+
+rmse = sqrt(mean(imdff.^2));
+psnr = 20*log10(255/rmse);

compute_rmse.m

@@ -0,0 +1,17 @@
+function [rmse] = compute_rmse(im1, im2)
+
+if size(im1, 3) == 3,
+    im1 = rgb2ycbcr(im1);
+    im1 = im1(:, :, 1);
+end
+
+if size(im2, 3) == 3,
+    im2 = rgb2ycbcr(im2);
+    im2 = im2(:, :, 1);
+end
+
+imdff = double(im1) - double(im2);
+imdff = imdff(:);
+
+rmse = sqrt(mean(imdff.^2));
+end

create_low.m

@@ -0,0 +1,14 @@
+function [s]=create_low(first_image)
+[row col]=size(first_image);
+%downsampled low images
+low1=first_image(1:2:end,1:2:end);
+%transform low images seperately
+%low1reference
+%low2
+low2 = imtranslate(low1,[0.8, 0.3]);
+%low3
+low3= imtranslate(low1,[0.25, 0.4]);
+%low4
+low4= imtranslate(low1,[0.6, 0.5]);
+s={low1,low2,low3,low4};
+end

deepconvolutional.m

@@ -0,0 +1,20 @@
+clear all;
+close all;
+addpath('C:\Users\PROCOMP7\Desktop\Superresolution\Matlab\superresolution_matlab');
+[FileName,PathName] = uigetfile('*.tif','Select GeoTIFF file');
+LR1= imread(strcat(PathName,FileName));
+[FileName,PathName] = uigetfile('*.tif','Select GeoTIFF file');
+LR2 = imread(strcat(PathName,FileName));
+[FileName,PathName] = uigetfile('*.tif','Select GeoTIFF file');
+LR3 = imread(strcat(PathName,FileName));
+
+b=zeros(5,5,3);
+filter=[0,0,1,0,0;0,1,2,1,0;1,2,-16,2,1;0,1,2,1,0;0,0,1,0,0];
+b(:,:,1)=filter;
+b(:,:,2)=filter;
+b(:,:,3)=filter;
+LR1up=imresize(LR1,2,'bicubic');
+LR1conv=convn(LR1up,b);
+
+figure;imshow(uint16(30*LR1conv(:,:,1)));
+figure;imshow(uint16(30*LR1(:,:,1)));

downfunction.m

@@ -0,0 +1,15 @@
+function D=downfunction(HR)
+HR=double(HR);
+[rowX, colX] = size(HR(:,:,1));
+
+D = sparse((rowX/2)*(colX/2),(rowX)*(colX));
+rowD = 0;
+for i = 1:2:rowX;
+    for j = 1:2:colX;
+        colD = colX*(i-1)+j;
+        rowD = rowD+1;
+        D(rowD,colD) = 1;
+    end
+end
+end
+

estimate_shift.m

@@ -0,0 +1,45 @@
+function [shiftx shifty] = estimate_shift(s,n)
+
+h = waitbar(0, 'Shift Estimation');
+set(h, 'Name', 'Please wait...');
+
+nr = length(s);
+delta_est=zeros(nr,2);
+p = [n n]; % only the central (aliasing-free) part of NxN pixels is used for shift estimation
+sz = size(s{1});
+sz=sz(:,1:2);
+S1 = fftshift(fft2(s{1})); % Fourier transform of the reference image
+for i=2:nr
+  waitbar(i/nr, h, 'Shift Estimation');
+  S2 = fftshift(fft2(s{i})); % Fourier transform of the image to be registered
+  S2(S2==0)=1e-10;
+  Q = S1./S2;
+  A = angle(Q); % phase difference between the two images
+
+  % determine the central part of the frequency spectrum to be used
+  beginy = floor(sz(1)/2)-p(1)+1;
+  endy = floor(sz(1)/2)+p(1)+1;
+  beginx = floor(sz(2)/2)-p(2)+1;
+  endx = floor(sz(2)/2)+p(2)+1;
+
+  % compute x and y coordinates of the pixels
+  x = ones(endy-beginy+1,1)*[beginx:endx];
+  x = x(:);
+  y = [beginy:endy]'*ones(1,endx-beginx+1);
+  y = y(:);
+  v = A(beginy:endy,beginx:endx);
+  v = v(:);
+
+  % compute the least squares solution for the slopes of the phase difference plane
+  M_A = [x y ones(length(x),1)];
+  r = M_A\v;
+  delta_est(i,:) = -[r(2) r(1)].*sz/2/pi;
+%   shiftx={abs(delta_est(2,2)),abs(delta_est(3,2)),abs(delta_est(4,2))};
+%   shifty={abs(delta_est(2,1)),abs(delta_est(3,1)),abs(delta_est(4,1))};
+for a=2:nr
+shiftx{a-1}=(delta_est(a,2));
+shifty{a-1}=(delta_est(a,1));
+end
+end
+
+close(h);

geotiff.m

@@ -0,0 +1,20 @@
+clear all;
+close all;
+addpath spotpan
+[a R]=geotiffread('spotpan1.tif');
+info = geotiffinfo('spotpan1.tif');
+int=imresize(a,2,'bicubic');
+Rint=R;
+Rint.RasterSize = size(int);
+info.GeoTIFFTags.GeoKeyDirectoryTag.ProjectedCSTypeGeoKey=32635;
+subImage = int;
+ xi = [1 - .5, 4092 + .5];
+ yi = [1 - .5, 2762 + .5];
+[xlimits, ylimits] = intrinsicToWorld(Rint, xi, yi);
+subR = Rint;
+subR.RasterSize = size(subImage);
+subR.XLimWorld = sort(xlimits);
+subR.YLimWorld = sort(ylimits);
+filename = 'int2.tif';
+geotiffwrite(filename, subImage, subR,  ...
+       'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);      

geowrite.m

@@ -0,0 +1,17 @@
+function result=geowrite(im,row,col);
+[FileName,PathName] = uigetfile('*.tif','Select GeoTIFF file');
+[a R]= geotiffread(strcat(PathName,FileName));
+info = geotiffinfo('strcat(PathName,FileName');
+xi = [1 - .5, col + .5];
+yi = [1 - .5, row + .5];
+[xlimits, ylimits] = intrinsicToWorld(R, xi, yi);
+subR = R;
+subR.RasterSize =[2*row 2*col];
+subR.XLimWorld = sort(xlimits);
+subR.YLimWorld = sort(ylimits);
+[FileName,PathName] = uigetfile('*.tif','Select GeoTIFF file');
+info= geotiffinfo(strcat(PathName,FileName));
+info.GeoTIFFTags.GeoKeyDirectoryTag.ProjectedCSTypeGeoKey=32635;
+im=uint16(im);
+geotiffwrite('im.tif', im, subR,'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
+end

interp.m

@@ -0,0 +1,9 @@
+function int=interp(s,nrofim);
+[row col band]=size(s{1});
+
+for i=1:nrofim
+  for l=1:band
+  int{i}(:,:,l)=imresize(s{i}(:,:,l),2,'bicubic');
+  end
+end
+ end

l1btv.m

@@ -0,0 +1,22 @@
+function [result]=l1btv(W,s,nrofim);
+initial=imresize(s{1},2,'bicubic');
+[row col band]=size(initial);
+[initialvector lr]=preparetoreconstruction(s,nrofim,initial);
+for i=1:band
+lowresvec=lr(:,i);
+iter=1;
+max_iter=3;
+x={0.5,-0.5};
+y={0.5,-0.5};
+I=speye(row*col,col*row);
+S=shiftmatrix(x,y,row,col,3);
+while iter < max_iter
+     initialvector(:,i)=initialvector(:,i)-0.001*W'*sign(W*initialvector(:,i)-lowresvec)+(I-S{2})*sign(initialvector(:,i)-S{3}*initialvector(:,i));
+     iter=iter+1;
+end
+sr(:,:,i)=reshape(real(initialvector(:,i)),col,row)';
+huberresult2(:,:,i) = medfilt2(sr(:,:,i),[3 3]);
+result(:,:,i)=huberresult2(2:end-1,2:end-1,i);
+end
+
+end

l1normlaplacien2.m

@@ -0,0 +1,19 @@
+function [l1norm2]=l1normlaplacien2(s,W,nrofim);
+initial=imresize(s{1},2,'bicubic');
+[row col band]=size(initial);
+[initialvector lr]=preparetoreconstruction(s,nrofim,initial);
+for i=1:band
+lowresvec=lr(:,i);
+iter=0;
+max_iter=7;
+operatormask=[-1,-1,-1;-1,8,-1;-1,-1,-1];
+ L=laplacienoperator(initial(:,:,i),operatormask);
+while iter < max_iter
+     initialvector(:,i)=initialvector(:,i)-0.01*W'*sign(W*initialvector(:,i)-lowresvec)+0.01*L*initialvector(:,i);
+     iter=iter+1;
+end
+sr(:,:,i)=reshape(real(initialvector(:,i)),col,row)';
+l1norm22(:,:,i) = medfilt2(sr(:,:,i),[3 3]);
+l1norm2(:,:,i)=l1norm22(2:end-1,2:end-1,i);
+end
+end

l2normlaplacien.m

@@ -0,0 +1,19 @@
+function [l2norm]=l2normlaplacien(s,W,nrofim);
+initial=imresize(s{1},2,'bicubic');
+[row col band]=size(initial);
+[initialvector lr]=preparetoreconstruction(s,nrofim,initial);
+for i=1:band
+lowresvec=lr(:,i);
+iter=0;
+max_iter=50;
+operatormask=[-1,-1,-1;-1,8,-1;-1,-1,-1];
+ L=laplacienoperator(initial(:,:,i),operatormask);
+while iter < max_iter
+     initialvector(:,i)=initialvector(:,i)-0.0001*W'*(W*initialvector(:,i)-lowresvec)+0.001*L*initialvector(:,i);
+     iter=iter+1;
+end
+sr(:,:,i)=reshape(real(initialvector(:,i)),col,row)';
+l2norm2(:,:,i) = medfilt2(sr(:,:,i),[3 3]);
+l2norm(:,:,i)=l2norm2(2:end-1,2:end-1,i);
+end
+end

l2normlaplacien2.m

@@ -0,0 +1,19 @@
+function [l2norm2]=l2normlaplacien2(s,W,nrofim);
+initial=imresize(s{1},2,'bicubic');
+[row col band]=size(initial);
+[initialvector lr]=preparetoreconstruction(s,nrofim,initial);
+for i=1:band
+lowresvec=lr(:,i);
+iter=0;
+max_iter=30;
+operatormask=[1,-2,1;-2,5,-2;1,-2,1];
+ L=laplacienoperator(initial(:,:,i),operatormask);
+while iter < max_iter
+     initialvector(:,i)=initialvector(:,i)+0.01*W'*(W*initialvector(:,i)-lowresvec)+0.001*L*initialvector(:,i);
+     iter=iter+1;
+end
+sr(:,:,i)=reshape(real(initialvector(:,i)),col,row)';
+l2norm22(:,:,i) = medfilt2(sr(:,:,i),[3 3]);
+l2norm2(:,:,i)=l2norm22(2:end-1,2:end-1,i);
+end
+end

l2normlaplacien3.m

@@ -0,0 +1,19 @@
+function [l2norm3]=l2normlaplacien3(s,W,nrofim);
+initial=imresize(s{1},2,'bicubic');
+[row col band]=size(initial);
+[initialvector lr]=preparetoreconstruction(s,nrofim,initial);
+for i=1:band
+lowresvec=lr(:,i);
+iter=0;
+max_iter=6;
+operatormask=[0,-1,0;-1,4,-1;0,-1,0];
+ L=laplacienoperator(initial(:,:,i),operatormask);
+while iter < max_iter
+     initialvector(:,i)=initialvector(:,i)+0.003*W'*(W*initialvector(:,i)-lowresvec)+0.01*L*initialvector(:,i);
+     iter=iter+1;
+end
+sr(:,:,i)=reshape(real(initialvector(:,i)),col,row)';
+l2norm32(:,:,i) = medfilt2(sr(:,:,i),[3 3]);
+l2norm3(:,:,i)=l2norm32(2:end-1,2:end-1,i);
+end
+end

laplacienoperator.m

@@ -0,0 +1,15 @@
+function [L]=laplacienoperator(Hr,mask) 
+
+[row_tmat col_tmat]=size(Hr);
+Laplacien = sparse(row_tmat*col_tmat,row_tmat*col_tmat);
+
+template = single(toeplitz_generator(1,2,3,row_tmat));
+for m = 1 : row_tmat
+    for n = 1 : row_tmat
+        if template(m,n) ~= 0
+            t = template(m,n);
+            Laplacien(col_tmat*(m-1)+1:col_tmat*(m), col_tmat*(n-1)+1:col_tmat*(n)) = toeplitz_generator(mask(t,1),mask(t,2),mask(t,3),col_tmat);
+        end
+    end
+L=Laplacien;
+end