-
Notifications
You must be signed in to change notification settings - Fork 57
/
ProjectDEM2Plane.m
108 lines (94 loc) · 4.04 KB
/
ProjectDEM2Plane.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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
% verison
% change to project the pixels from the potential shadow layer. Mar., 3, 2018 by Shi Qiu
% sometimes the projection will generate Nan values. fixed this bug by
% Feb., 13, 2018 by Shi Qiu
% fix the bug some pixels out of image. Shi 9/12/2017
%
function [dem_plane_r,dem_plane_c]= ProjectDEM2Plane(dim,pshadow,dem,dem_base,...
sun_elevation_rad,sun_tazi_rad,sun_azimuth_deg,dim_expd,resolu)
% DEM is vertically obervated by scene.
% assume resolu.x=resolu.y
spatial_resol=resolu(1);
% get the location of orginal pixel
dem_loc=zeros(dim,'uint8')+1;
[loc_i,loc_j]=find(dem_loc);
% must project all pixels.
% [loc_i,loc_j]=find(pshadow==1);
% loc_i=uint16(loc_i);
% loc_j=uint16(loc_j);
% if isequal(dem,0) % when no DEMs, no cloud shadow shape correction.
if isempty(dem) % when no DEMs, no cloud shadow shape correction.
loc_i_plane=loc_i;
loc_j_plane=loc_j;
else % when having DEMs
% dem_dif=dem(pshadow==1)-dem_base;
dem_dif=dem-dem_base;
dem_dif=dem_dif(:);
clear dem_loc dem dem_base;
% sun angle geometry
% sun_elevation_deg=90-sun_zenith_deg;
% sun_elevation_rad=deg2rad(sun_elevation_deg);
d_ij=dem_dif./(spatial_resol*tan(sun_elevation_rad));
% Sun_tazi=sun_azimuth_deg-90;
% sun_tazi_rad=deg2rad(Sun_tazi);
% if sun_azimuth_deg< 180
if sun_azimuth_deg< 180
di=0-d_ij.*sin(sun_tazi_rad); % i for row, Y
dj=0-d_ij.*cos(sun_tazi_rad); % j for col, X
else
di=d_ij.*sin(sun_tazi_rad); % i for row, Y
dj=d_ij.*cos(sun_tazi_rad); % j for col, X
end
di=double(di);
dj=double(dj);
clear dem_dif d_ij Sun_tazi sun_elevation_rad sun_zenith_deg Sun_tazi sun_tazi_rad;
% original location adds the bias.
loc_i_plane=round(loc_i+di);
loc_j_plane=round(loc_j+dj);
clear di dj;
end
% dim_expd=1000;% 1000 buffer
dim_plane_expd=dim+(2*dim_expd);
% expand the range of observations
loc_i_plane_expd=loc_i_plane+dim_expd;
loc_j_plane_expd=loc_j_plane+dim_expd;
% remove the pixels out of box.
% % out_ids=loc_i_plane_expd(:)<1|loc_j_plane_expd(:)<1|...
% % loc_i_plane_expd(:)>dim_plane_expd(1)|loc_j_plane_expd(:)>dim_plane_expd(2);
% % loc_i_plane_expd(out_ids)=[];
% % loc_j_plane_expd(out_ids)=[];
% % loc_i(out_ids)=[];
% % loc_j(out_ids)=[];
loc_i_plane_expd(loc_i_plane_expd(:)<1)=1;
loc_j_plane_expd(loc_j_plane_expd(:)<1)=1;
loc_i_plane_expd(loc_i_plane_expd(:)>dim_plane_expd(1))=dim_plane_expd(1);
loc_j_plane_expd(loc_j_plane_expd(:)>dim_plane_expd(2))=dim_plane_expd(2);
clear loc_i_plane loc_j_plane dj dim_expd;
dem_plane_r=zeros(dim_plane_expd);
dem_plane_c=zeros(dim_plane_expd);
% dem_plane_r = NaN(dim_plane_expd,'double');
% dem_plane_c = NaN(dim_plane_expd,'double');
% recorders.
tmp_id_plane=sub2ind(dim_plane_expd,double(loc_i_plane_expd),double(loc_j_plane_expd));
clear loc_i_plane_expd loc_j_plane_expd;
% sometimes there may be some Nan values!
nan_ids=isnan(tmp_id_plane);
if sum(nan_ids(:))>0
tmp_id_plane(nan_ids)=[];
loc_i(nan_ids)=[];
loc_j(nan_ids)=[];
end
dem_plane_r(tmp_id_plane)=loc_i;
dem_plane_c(tmp_id_plane)=loc_j;
clear loc_i loc_j tmp_id_plane;
% filled the holes due to the dicretization projection.
% dem_plane_r=fillmissing(dem_plane_r,'nearest',1);% by column.
% dem_plane_c=fillmissing(dem_plane_c,'nearest',1);% by column.
% fill remainings as 0 value.
% dem_plane_r=fillmissing(dem_plane_r,'constant',0);% by column.
% dem_plane_c=fillmissing(dem_plane_c,'constant',0);% by column.
% dem_plane_r=round(dem_plane_r);
% dem_plane_c=round(dem_plane_c);
% dem_plane_r=imfill(dem_plane_r,'holes');
% dem_plane_c=imfill(dem_plane_c,'holes');
end