-
Notifications
You must be signed in to change notification settings - Fork 1
/
poisson_solver_function.m
37 lines (36 loc) · 1.58 KB
/
poisson_solver_function.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
%% According to http://web.media.mit.edu/~raskar/photo/code.pdf
function [img_direct] = poisson_solver_function(gx,gy,boundary_image);
% function [img_direct] = poisson_solver_function(gx,gy,boundary_image)
% Inputs; Gx and Gy -> Gradients
% Boundary Image -> Boundary image intensities
% Gx Gy and boundary image should be of same size
[H,W] = size(boundary_image);
gxx = zeros(H,W); gyy = zeros(H,W);
f = zeros(H,W); j = 1:H-1; k = 1:W-1;
% Laplacian
gyy(j+1,k) = gy(j+1,k) - gy(j,k); gxx(j,k+1) = gx(j,k+1) - gx(j,k);
f = gxx + gyy; clear j k gxx gyy gyyd gxxd
% boundary image contains image intensities at boundaries
boundary_image(2:end-1,2:end-1) = 0;
disp('Solving Poisson Equation Using DST'); tic
j = 2:H-1; k = 2:W-1; f_bp = zeros(H,W);
f_bp(j,k) = -4*boundary_image(j,k) + boundary_image(j,k+1) +...
boundary_image(j,k-1) + boundary_image(j-1,k) + boundary_image(j+1,k);
clear j k
f1 = f - reshape(f_bp,H,W);% subtract boundary points contribution
clear f_bp f
% DST Sine Transform algo starts here
f2 = f1(2:end-1,2:end-1); clear f1
%compute sine transform
tt = dst(f2); f2sin = dst(tt')'; clear f2
%compute Eigen Values
[x,y] = meshgrid(1:W-2,1:H-2); denom = (2*cos(pi*x/(W-1))-2) + (2*cos(pi*y/(H-1)) - 2) ;
%divide
f3 = f2sin./denom; clear f2sin x y
%compute Inverse Sine Transform
tt = idst(f3); clear f3; img_tt = idst(tt')'; clear tt
time_used = toc; disp(sprintf('Time for Poisson Reconstruction = %f secs',time_used));
% put solution in inner points; outer points obtained from boundary image
img_direct = boundary_image;
img_direct(2:end-1,2:end-1) = 0;
img_direct(2:end-1,2:end-1) = img_tt;