
% Create sharp depth gradient field from segmented ratio images
function CreateDepthGrad

% ----- input: segmented ratio images ------

r_up = imread('input/up.bmp');
r_up = double(r_up); r_up = r_up ./ 255;
r_right = imread('input/right.bmp');
r_right = double(r_right); r_right = r_right ./ 255;
r_down = imread('input/down.bmp');
r_down = double(r_down); r_down = r_down ./ 255;
r_left = imread('input/left.bmp');
r_left = double(r_left); r_left = r_left ./ 255;


% parameters to set:
% ----------------------------------------
thr1 = 10/255; thr2 = 10/255;
max_shadow_width = 100;
filename = 'depthGrad.txt';
% ----------------------------------------

% original images
img_up = imread('input/uphand.jpg');img_right = imread( 'input/righthand.jpg' );img_down = imread( 'input/downhand.jpg' );img_left = imread( 'input/lefthand.jpg' );

% make gray scale images
gimg_up = mean( img_up, 3 ); mean_up = mean(gimg_up(:));
gimg_right = mean( img_right, 3 ); mean_right = mean(gimg_right(:));
gimg_down = mean( img_down, 3 ); mean_down = mean(gimg_down(:));
gimg_left = mean( img_left, 3 ); mean_left = mean(gimg_left(:));

all_mean = mean( [mean_up; mean_right; mean_down; mean_left] );

% normalize intensities
gimg_up = gimg_up *all_mean/mean_up;
gimg_right = gimg_right *all_mean/mean_right;
gimg_down = gimg_down *all_mean/mean_down;
gimg_left = gimg_left *all_mean/mean_left;


depth_grad_x = zeros(size(gimg_left));
depth_grad_y = zeros(size(gimg_left));

[depth_grad_x,shadow_left] = FindEdgesLeftFlash(r_left,img_left,depth_grad_x,thr1,thr2,max_shadow_width);
[depth_grad_x,shadow_right] = FindEdgesRightFlash(r_right,img_right,depth_grad_x,thr1,thr2,max_shadow_width);
[depth_grad_y,shadow_up] = FindEdgesUpFlash(r_up,img_up,depth_grad_y,thr1,thr2,max_shadow_width);
[depth_grad_y,shadow_down] = FindEdgesDownFlash(r_down,img_down,depth_grad_y,thr1,thr2,max_shadow_width);

depth_grad_x = Padding(depth_grad_x,1);
depth_grad_y = Padding(depth_grad_y,2);

SaveDepthGrad(depth_grad_x,depth_grad_y,filename);

edges = abs(depth_grad_x)|abs(depth_grad_y);
edges(find(edges>0)) = 255;
imwrite( edges, 'edges.bmp' );



% input:  r_left: left ratio image
%         img_left: just so overwrite the shadow width in red over it
%         depth_grad: previous depth gradient image which will be updated
%         thr1,thr2: threshold for inner and outer shadow boundary, respectively
%         max_shadow_width
% output: depth_grad: updated depth gradient map
%         shadow_left: left image with shadows in red

function [depth_grad,shadow_left] = FindEdgesLeftFlash(r_left,img_left,depth_grad,thr1,thr2,max_shadow_width)

siz = size(r_left);
height = siz(1);
width  = siz(2);
shadow_left = img_left;

for i = 1:height
    for j = 1:(width-1)
        dif = r_left(i,j+1)- r_left(i,j);
        if (dif < 0) & (abs(dif)>thr1) % found depth edge
            % now walk along epipolar line to compute shadow width
            pos_j = j;
            shadow_width = 0;
            while ((r_left(i,pos_j+shadow_width+1) - r_left(i,pos_j+shadow_width)) < thr2) & (shadow_width < max_shadow_width) & (pos_j+shadow_width+2<width)
                shadow_width = shadow_width + 1;
                shadow_left(i,pos_j+shadow_width,1) = 255;
                shadow_left(i,pos_j+shadow_width,2) = 0;
                shadow_left(i,pos_j+shadow_width,3) = 0;
            end
            if shadow_width == max_shadow_width
                shadow_width = 1;
            end
            depth_grad(i,pos_j) = -shadow_width;    
        end
    end
end
figure(1)
imshow(shadow_left);


%********************************************************************************************************************************

function [depth_grad,shadow_right] = FindEdgesRightFlash(r_right,img_right,depth_grad,thr1,thr2,max_shadow_width)

siz = size(r_right);
height = siz(1);
width  = siz(2);
shadow_right = img_right;

for i = 1:height
    for j = width:-1:2
        dif = r_right(i,j-1)- r_right(i,j);
        if (dif < 0) & (abs(dif)>thr1) % found depth edge
            % now walk along epipolar line to compute shadow width
            pos_j = j;
            shadow_width = 0;
            while ((r_right(i,pos_j-shadow_width-1) - r_right(i,pos_j-shadow_width)) < thr2) & (shadow_width < max_shadow_width) & (pos_j-shadow_width-2>1)
                shadow_width = shadow_width + 1;
                shadow_right(i,pos_j-shadow_width,1) = 255;
                shadow_right(i,pos_j-shadow_width,2) = 0;
                shadow_right(i,pos_j-shadow_width,3) = 0;
            end
            if shadow_width == max_shadow_width
                shadow_width = 1;
            end
            depth_grad(i,pos_j) = shadow_width;    
        end
    end
end
figure(2)
imshow(shadow_right);

%********************************************************************************************************************************

function [depth_grad,shadow_up] = FindEdgesUpFlash(r_up,img_up,depth_grad,thr1,thr2,max_shadow_width)

siz = size(r_up);
height = siz(1);
width  = siz(2);
shadow_up = img_up;

for j = 1:width
    for i = 1:(height-1)
        dif = r_up(i+1,j)- r_up(i,j);
        if (dif < 0) & (abs(dif)>thr1) % found depth edge
            % now walk along epipolar line to compute shadow width
            pos_i = i;
            shadow_width = 0;
            while ((r_up(pos_i+shadow_width+1,j) - r_up(pos_i+shadow_width,j)) < thr2) & (shadow_width < max_shadow_width) & (pos_i+shadow_width+2<height)
                shadow_width = shadow_width + 1;
                shadow_up(pos_i+shadow_width,j,1) = 255;
                shadow_up(pos_i+shadow_width,j,2) = 0;
                shadow_up(pos_i+shadow_width,j,3) = 0;
            end
            if shadow_width == max_shadow_width
                shadow_width = 1;
            end
            depth_grad(pos_i,j) = -shadow_width;    
        end
    end
end
figure(3)
imshow(shadow_up);

%********************************************************************************************************************************

function [depth_grad,shadow_down] = FindEdgesDownFlash(r_down,img_down,depth_grad,thr1,thr2,max_shadow_width)

siz = size(r_down);
height = siz(1);
width  = siz(2);
shadow_down = img_down;

for j = 1:width
    for i = height:-1:2
    
        dif = r_down(i-1,j)- r_down(i,j);
        if (dif < 0) & (abs(dif)>thr1) % found depth edge
            % now walk along epipolar line to compute shadow width
            pos_i = i;
            shadow_width = 0;
            while ((r_down(pos_i-shadow_width-1,j) - r_down(pos_i-shadow_width,j)) < thr2) & (shadow_width < max_shadow_width) & (pos_i-shadow_width-2>1)
                shadow_width = shadow_width + 1;
                shadow_down(pos_i-shadow_width,j,1) = 255;
                shadow_down(pos_i-shadow_width,j,2) = 0;
                shadow_down(pos_i-shadow_width,j,3) = 0;
            end
            if shadow_width == max_shadow_width
                shadow_width = 1;
            end
            depth_grad(pos_i,j) = shadow_width;    
        end
    end
end
figure(4)
imshow(shadow_down);


%********************************************************************************************************************************

%********************************************************************************************************************************


function result = Padding(depth_grad,direction)

siz = size(depth_grad);
maxsiz = max(siz);
n = 1;
while (2^n < maxsiz)
    n = n + 1;
end

new_siz = 2^n+2;
result = zeros(new_siz,new_siz);

pos_x = (new_siz-siz(2))/2;
pos_y = (new_siz-siz(1))/2;
result(pos_y:(pos_y+siz(1)-1), pos_x:(pos_x+siz(2)-1)) = depth_grad(:,:);


  

% Write depth gradient map to be read in VC++ program for Poisson reconstruction
function SaveDepthGrad(depth_grad_x,depth_grad_y,filename)

fp = fopen(filename,'wt');
siz = size(depth_grad_x);
fprintf(fp,'%d %d \n',siz(2),siz(1));
for i=1:siz(1)
  for j=1:siz(2)
      fprintf(fp,'(%d,%d) ',depth_grad_x(i,j),depth_grad_y(i,j));
  end
  fprintf(fp,'\n');
end

fclose(fp);

