demo.m

Matlab Code for Poisson Image Reconstruction from Image Gradients

%%%%%%%% Read in and store  images %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear; 
p = 2; 

addpath('../')
load face_adyoron_1_GS.mat
frame = face_adyoron_1_GS;

M = size(frame(:,:,1),1); 
N = size(frame(:,:,1),2); 
[m,n] = meshgrid(1:M,1:N); 

%%%%%%%%% Find Translation Parameters %%%%%%%%%%%%%%%%%%%%% 
K(1,:) = zeros(1,4); 

for num= 2:p, 
  K(num,:) = affine(frame(:,:,num),frame(:,:,1)); 
end 

deltaK = K(:,3:4); 

%%%%%%%%%% Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Lu = 1; 
Lv = 1; 
Tx = 2; 
Ty = 2; 

wsx = 1./Tx; 
wsy = 1./Ty; 

F1 = fft2(frame(:,:,1)); 
F2 = fft2(frame(:,:,2)); 

F = zeros(2*Lu*M,2*Lv*N); 

%%%%%%%%%%% Calculating values in matrix F  %%%%%%%%%%%%%%%%%% 
for n = 1:N, 
  disp(n); 
  for m = 1:M, 
    Gmn = [F1(m,n);F2(m,n)]; 

    for k=1:p, 
      for r = 1:4 
 PHImn(k,r) = exp(j*2*pi.*(deltaK(k,1)*(m./M*Tx)-(Lu/Tx)+deltaK(k,2)*(n./N*Ty)-(Lv/Ty)))*... 
     exp(j*2*pi.*((deltaK(k,1)./Tx)*rem(r-1,2*Lu)+(deltaK(k,2)*floor((r-1)./2*Lv)))); 
      end 
    end 
  
    Fmn = pinv(PHImn)*Gmn; 
  
    F(m,n) = Fmn(4,1); 
    F(m+M,n) = Fmn(2,1); 
    F(m,n+N) = Fmn(3,1); 
    F(m+M,n+N) = Fmn(1,1); 
  
  end 
end 
  

%%%%%%%%%%% Get f from F and display  %%%%%%%%%%%%%%%%%%%%%%% 
image = abs(ifft2(fftshift(F))); 
figure; 
%showIm(image); 
imagesc(image); 
colormap('gray')