function [Sest,lightsourcemask,lightsourcecolor]=getShadowMask(im,depth_data,normal_data,segment)
addpath('utils/matlab');
addpath('utils/npy-matlab/npy-matlab');
%% z->front y->up x->right
%% Real/Synthetic x->front y->up z->left
% path = 'dataset\synthetic\livingroom';
% im=im2double(imread([path '\up.png']));
% normal_data = im2double(exrread([path '\normal.exr']));
% normal_data_x = normal_data(:,:,1);normal_data_y = normal_data(:,:,2);normal_data_z = normal_data(:,:,3);
% normal_data = cat(3, -normal_data_z, normal_data_y, normal_data_x);
% depth_data = exrread([path '\depth.exr']);
% depth_data = depth_data(:,:,1);
%
%% y->front z->up x->right
% im=im2double(imread(['dataset\' path '\data\' path '_color_0_Left_Down_0.0.png']));
% normal_data = im2double(exrread(['dataset\' path '\data\' path '_normal_0_Left_Down_0.0.exr']));
% normal_data_x = normal_data(:,:,1);normal_data_y = normal_data(:,:,2);normal_data_z = normal_data(:,:,3);
% normal_data = cat(3, normal_data_x, normal_data_z, normal_data_y);
% depth_data = im2double(exrread(['dataset\' path '\data\' path '_depth_0_Left_Down_0.0.exr']));
% depth_data = depth_data(:,:,1);
%% Structure3D
% y->up x->front z-> left
% path = 'E:\Python\LRG_360Panoramic-main\data\dataset\72\raw';
% im=im2double(imread([path '\input.png']));
% normal_data = double(imread([path '\output_normal.png'])) ./128 - 1;
% normal_data_x = normal_data(:,:,1);normal_data_y = normal_data(:,:,2);normal_data_z = normal_data(:,:,3);
% normal_data = cat(3, -normal_data_z, normal_data_y, normal_data_x);
% % depth_data = imread([path '\depth.png']);
% depth_data = readNPY([path '\depth.npy']);
% depth_data = double(im2gray(depth_data));
normal_data_x = normal_data(:,:,1);normal_data_y = normal_data(:,:,2);normal_data_z = normal_data(:,:,3);
normal_data = cat(3, -normal_data_z, normal_data_y, normal_data_x);
depth_data = double(im2gray(depth_data));
[h,w,~] = size(depth_data);
gray = im2gray(im);
% bright = gray;bright(bright<0.6)=0;
bright = (im - Shen2009(im));
%bright = im2gray(bright);
ori_bright = bright;
bright = im2gray((bright - min(min(bright)))./(max(max(bright)) - min(min(bright))));
stmp = sort(bright(:),'descend');
threshold = stmp(round(h*w*0.02));
tmp=bright;tmp(tmp>threshold)=1;%0.8
mask = tmp;
% stmp = sort(gray(:),'descend');
% threshold = stmp(round(h*w*0.02));
% tmp=gray;tmp(tmp>threshold)=1;%0.8
% mask = tmp;
[points, color_points] = get_pointcloud(im, depth_data);
normal_data = reshape(pagetranspose(normal_data),h*w,3);
normal_data = normalize(normal_data,2,'norm');
normal_data(isnan(normal_data)) = 0;
[x,y] = find(mask==1);
light = double(points(x*w + y - w,:)) * 0.95;
normal = double(normal_data(x*w + y - w,:));
ls = size(light,1);
tmpp = sort(points(:,2),'descend');
%% Check Mask
count_light = zeros(ls,1);
count_image = zeros(h*w,1);
if ls > 2000
loopind = floor(linspace(1,ls,2000));
else
loopind = 1:ls;
end
for l=loopind
lightpoint = light(l,:);
view = normalize(-light(l,:),'norm');
lightnormal = normal(l,:);
if (norm(lightnormal) <= 1e-6)
continue
end
ref = 2 * (view * lightnormal') * lightnormal - view;
tmpvec = normalize(light - lightpoint,2,'norm');
ind = tmpvec*ref' > 0.99;
count_light(ind) = count_light(ind) + 1;
tmpvec = normalize(points - lightpoint,2,'norm');
ind = tmpvec*ref' > 0.95;
count_image(ind) = count_image(ind) + 1;
% if lightnormal(2) < -0.98 || lightpoint(2) > tmpp(round(h*w*0.2))
% count(l) = 9999;
% %count(x(l) * w + y(l) - w) = 9999;
% end
%
% if lightnormal(2) > 0.98
% count(l) = -999;
% %count(x(l) * w + y(l) - w) = -999;
% end
end
ind = x*w + y - w;
tmp_light = zeros(h*w,1);
tmp_light(ind(count_light > 999)) = 1;
count(count_light >999)=0;
tmp_light(ind(count_light>mean(count_light))) = 1;
% imshow(reshape(tmp,w,h)')
tmp_image = sort(count_image,'descend');
count_image(count_image < tmp_image(round(h*w*0.1)) ) = 0;
count_image = count_image ./ tmp_image(round(h*w*0.0025));
tmp_image = count_image;
[M,~] = max(tmp_image,[],'all');
tmp = reshape(tmp_image,w,h)';
tmp(tmp<M*0.8)=0;
tmp_image = tmp.*reshape(tmp_image,w,h)';
tmp = (reshape(tmp_light,w,h)' + tmp_image)/2;
mask = tmp;mask(mask>0)=1;
mask(normal_data_y>0.5) = 0;
mask = AreaGrowth(im, mask);
mask(normal_data_y>0.8) = 0;
%mask = Dilate(im, mask);
mask = mask.*im;
% for hdr
% tmp = zeros(size(im));tmp(im>max(im,[],'all') * 0.3)=1;
% mask = tmp.*im;
% tmp = mean(tmp,3);
lightsourcecolor = sum(mask,[1 2])./sum(mask>0,[1 2]);
upind = points(:,2) > tmpp(round(h*w*0.2));
normal_data(upind,:) = repmat([0 -1 0],[sum(upind>0),1]);
% mask(floor(h/2):h,:)=0;
% mask(bright==0)=0;
lightsourcemask = mask;
mask = mean(mask,3);
%lightsourcemask(mask<0.5) = 0;
[x,y] = find(mask>0);
light = double(points(x*w + y - w,:))* 0.95;
normal = double(normal_data(x*w + y - w,:));
rowbright = reshape(pagetranspose(mask),h*w,1);
intensity = double(rowbright(x*w+y-w,:));
intensity(intensity>2) = 2;
ls = size(light,1);
shadowmap_h = h;
shadowmap_w = w;
shadow = zeros(h*w,1);
count = 0;
if ls > 20
loopind = floor(linspace(1,ls,20));
else
loopind = 1:ls;
end
for l=loopind
lightpoint = light(l,:);
lightnormal = normal(l,:);
dis = sqrt(sum((points - lightpoint).^2,2));
theta_phi = world_2_angle_coordinate(points - lightpoint);
i_x_y = angle_2_image_coordinate(theta_phi, shadowmap_w, shadowmap_h);
[len,~] = size(points);
new_image = zeros(h*w,1);
norm_dir = normalize(lightpoint - points,2,'norm');
norm_view = normalize(points,2,'norm');
norm_dis = (dis - mean(dis))/(mean(dis) - min(dis));
%Gen shadow map
envmap = zeros(shadowmap_h,shadowmap_w,3);
tmp =zeros(shadowmap_h,shadowmap_w);
%% for point cloud method
shadowmap = ones(shadowmap_h,shadowmap_w) * 999999999;
for i = 1:len
i_x = i_x_y(i,1) + 1;
i_y = i_x_y(i,2) + 1;
if dis(i,1) < shadowmap(i_y,i_x)
shadowmap(i_y, i_x) = dis(i,1);
tmp(i_y,i_x) = i;
envmap(i_y, i_x,:) = color_points(i,:);
end
end
%% Compare depth
for i = 1:len
i_x = i_x_y(i,1) + 1;
i_y = i_x_y(i,2) + 1;
%0.1 control the shadow area
if dis(i,1) > shadowmap(i_y,i_x) + 0.2 && shadowmap(i_y,i_x) > 0.1% && ~(normal_data(i,2) < -0.99)
new_image(i) = 0;
else
vec = norm_dir(i,:);
N = normal_data(i,:);
view_dir = norm_view(i,:);
halfv = (view_dir + vec);
%intensity(l,1)
new_image(i) = intensity(l,1) / ((dis(i,1) * 0.2 + 0.5).^2) * (max(0,vec*N') + max(0,N * halfv').^5) ;
end
end
shadow = shadow + new_image;
count = count + 1;
end
res = reshape(shadow,w,h)'./count;
ori_res = res;
gray = im2gray(im);
[h,w,~] = size(im);
% weight = res ./ mean(res,'all').*0.5;
% res = (weight.*res + (1 - weight).*gray);
weight = zeros([h,w]);
top = floor(h*100/512); bottom = floor(h*410/512);
weight(top:bottom,:) = 1;
weight(1:top,:) = repmat(linspace(0,1,top)',[1,w]);
weight(bottom:end,:) = repmat(linspace(1,0,h - bottom + 1)',[1,w]);
res = res .* weight + gray.*(1-weight);
% 缩小图像范围
res = res - mean(res,'all');
res = res*0.5;
res = res + mean(gray,'all');
res = repmat(res,[1 1 3]).*lightsourcecolor;
Sest = res + ori_bright;
end
% imwrite(res,['shadow_' path '_' method '_auto.png']);
% exrwrite(res,['shadow_' path '_' method '_auto.exr']);
function s = CheckVisibility(p)
global mask;
global points;
global light;
[h,w,~] = size(mask);
s = 0;
count = 0;
ls = size(light,2);
a = sum((points - p).^2,3);
for l=1:10:ls
q = light(1,l,:);
b = sum((points - q).^2,3);
c = sqrt(sum((p - q).^2,3));
dis = sqrt(a.*b - ((a + b - c.^2)/2).^2)./c;
dis = real(dis);
if sum(sum(dis < 0.005)) <= 2
s = s + 1;
end
count = count + 1;
end
s = s / count;
end
function [camera_points, color_points] = get_pointcloud(color_data, depth_data)
[h,w,~] = size(color_data);
[pixel_x, pixel_y] = meshgrid(0:w-1,0:h-1);
theta = pixel_x * 2.0 * pi / w - pi;
phi = (h - pixel_y-0.5) * pi / h - pi / 2;
camera_points_x = depth_data .* cos(phi) .* sin(theta);
camera_points_y = depth_data .* sin(phi);
camera_points_z = depth_data .* cos(phi) .* cos(theta);
camera_points = cat(3,camera_points_x,camera_points_y,camera_points_z);
camera_points = reshape(pagetranspose(camera_points),h*w,3);
color_points = reshape(pagetranspose(color_data),h*w,3);
% global remove_up_down
% if remove_up_down
% depth_data(1:72, :) = 0;
% depth_data(441:end, :) = 0;
% end
% valid_depth_ind = depth_data > 0;
% camera_points = camera_points(valid_depth_ind, :);
% color_points = color_points(valid_depth_ind, :);
end
function [angle] = world_2_angle_coordinate(world_p)
target_p_x = world_p(:,1);
target_p_y = world_p(:,2);
target_p_z = world_p(:,3);
r = sqrt(sum(world_p .^ 2, 2));
phi = asin(target_p_y ./ r);
theta = asin(target_p_x ./ sqrt(target_p_x .^ 2 + target_p_z .^ 2+1e-10));
idx = (target_p_x < 0) & (target_p_z < 0);
theta(idx) = - theta(idx) - pi;
idx = (target_p_x > 0) & (target_p_z < 0);
theta(idx) = - theta(idx) + pi;
angle = cat(3,theta,phi);
end
function [coord] = angle_2_image_coordinate(theta_phi, width, height)
theta = theta_phi(:, 1);
phi = theta_phi(:, 2);
x = (theta + pi) * width / 2 / pi;
y = height - (phi + pi / 2) * height / pi - 0.5;
i_x = mod(int32(round(x)),width);
i_y = mod(int32(round(y)),height);
coord = [i_x i_y];
end
function [inter] = img_inter_func(valid_index, values, method)
[h,w,~] = size(values);
inter = zeros(size(values));
[X, Y] = meshgrid(1:h, 1:w);
[vx,vy] = find(valid_index);
for i = 1:3
chn = values(:, :, i);
Ti = griddata(vx, vy, chn(valid_index), X, Y, method);
inter(:, :, i) = Ti';
end
end
function [C, label, J] = kmeans(I, k)
[m, n, p] = size(I);%图片的大小m*n,p代表RGB三层
X = reshape(double(I), m*n, p);
rng('default');
C = X(randperm(m*n, k), :);%随机选三个聚类中心
J_prev = inf; iter = 0; J = []; tol = 1e-11;%容忍度tol,inf为无穷大
while true,
iter = iter + 1;
dist = sum(X.^2, 2)*ones(1, k) + (sum(C.^2, 2)*ones(1, m*n))' - 2*X*C';%计算图片中各个点到K个聚类中心的距离
[~,label] = min(dist, [], 2) ; %label记录最小值的行数
for i = 1:k,
C(i, :) = mean(X(label == i , :)); %取新的k个聚类中心
end
J_cur = sum(sum((X - C(label, :)).^2, 2));%距离之和
J = [J, J_cur];
display(sprintf('#iteration: %03d, objective fcn: %f', iter, J_cur));
if norm(J_cur-J_prev, 'fro') < tol,% A和A‘的积的对角线和的平方根,即sqrt(sum(diag(A'*A))),本次与上次距离之差
break;
end
if (iter==10),% A和A‘的积的对角线和的平方根,即sqrt(sum(diag(A'*A))),本次与上次距离之差
break;
end
J_prev = J_cur;
end
end
function im = Dilate(image, mask)
[M,I] = max(image.*mask,[],'all');
[h,w] = size(image);
% mask = zeros(h,w);
% mask(mod(I,h)+1,floor(I/h)+1)=1;
mask(mask>0)=1;
flag = true;
% se = strel("disk",5);
% mask = imerode(mask,se);
se = strel("disk",1);
count = 20;
while flag
mask = imdilate(mask,se);
tmpr = image(:,:,1).*mask;
tmpr(mask<1)=[];
tmpg = image(:,:,2).*mask;
tmpg(mask<1)=[];
tmpb = image(:,:,3).*mask;
tmpb(mask<1)=[];
m = (tmpr+tmpg+tmpb) / 3;
if max(mink(m,ceil(size(m,2)*0.25))) < M*0.8 || count==0
flag =false;
end
count = count - 1;
end
if count==0
disp('Early stop dilate.')
end
im = mask;
% tmp = image;
% tmp(tmp<M*2/3)=0;
% im = tmp.*image;
end
function im = AreaGrowth(image, mask)
[M,~] = max(image.*mask,[],'all');
se = strel("disk",1);
I = im2gray(image); % 输入图像
Ir = image(:,:,1);
Ig = image(:,:,2);
Ib = image(:,:,3);
[Gr,~] = imgradient(Ir);
[Gg,~] = imgradient(Ig);
[Gb,~] = imgradient(Ib);
G = max(cat(3,Gr,Gg,Gb),[],3);
[Mr,~] = max(Ir.*mask,[],'all');
[Mg,~] = max(Ig.*mask,[],'all');
[Mb,~] = max(Ib.*mask,[],'all');
Mask = imerode(mask,se); % 原始mask
S = find(Mask); % 种子点集
% S = S(1:20:size(S));
D = S; % 扩散点集
res = false(size(I));
res(S) = true;
step = 0;
while ~isempty(S) && step < 50
step = step + 1;
while ~isempty(S)
p = S(1);
S(1) = [];
[row,col] = ind2sub(size(I),p);
rowMin = max(row-1,1);
rowMax = min(row+1,size(I,1));
colMin = max(col-1,1);
colMax = min(col+1,size(I,2));
[neighsx,neighsy] = meshgrid(rowMin:rowMax, colMin:colMax);
neighsx = neighsx(:);
neighsy = neighsy(:);
for n=1:size(neighsx,1)
q = sub2ind(size(I),neighsx(n),neighsy(n));
if res(q)==false && G(q) < 0.05 && Ir(q) > Mr * 0.8 && Ig(q) > Mg * 0.8 && Ib(q) > Mb * 0.8
D = [D; q];
res(q) = true;
end
end
end
S = D;
D = [];
end
im = res;
end