-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathFcn_StemBases.m
executable file
·177 lines (167 loc) · 7.18 KB
/
Fcn_StemBases.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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
function Fcn_StemBases(StemLinePathName, StemLineFileName, ...
DEMPathName, DEMFileName, ...
StemBasePathName, StemBaseFileName)
% Registration, step 3, calculate the intersection of stem line and DEM
% that is used as the base point of stems. And output the intersection
% point to a text file.
% SYNTAX: Fcn_StemBases(StemLinePathName, StemLineFileName, ...
% DEMPathName, DEMFileName, ...
% StemBasePathName, StemBaseFileName)
% INPUT:
% 1. StemLinePathName, StemLineFileName: Text files of fitted lines of stems. (FittedStemLine_)
% Format:
% m, n, p, x0, y0, z0, TreeNO, number_of_stem_centers_in_fitting_line_finally
% ([m, n, p] is direction of line, [x0, y0, z0] is coordinates of one
% point on line)
% 2. DEMPathName, DEMFileName: Text file of DEM derived from the lowest
% points.
% File format: (EXAMPLE)
% "ncols 181
% nrows 181
% xllcorner -90.496391
% yllcorner -86.518433
% cellsize 1
% NODATA_value -9999
% [nrows by ncols matrix]"
% 3. StemBasePathName, StemBaseFileName: Text file for outputting stem
% bases.
% File format:
% #1: 'x, y, z, TreeNO'
% #2-#end: x, y, z, TreeNO
% OUTPUT:
% REQUIRED ROUTINES:
% METHODS AND PROCESS:
ParameterFileName=['Parameters_Fcn_StemBases_', StemBaseFileName];
fid=fopen(fullfile(StemBasePathName, ParameterFileName), 'w');
fprintf(fid, [...
'the input stem line file: %s\r\n', ...
'the DEM file: %s\r\n' ...
'the output file for stem base: %s' ...
], ...
fullfile(StemLinePathName, StemLineFileName), ...
fullfile(DEMPathName, DEMFileName), ...
fullfile(StemBasePathName, StemBaseFileName));
fclose(fid);
% ---read from stem line file---
fid=fopen(fullfile(StemLinePathName, StemLineFileName), 'r');
stemline=textscan(fid, '%f %f %f %f %f %f %f %f');
stemline=cell2mat(stemline);
fclose(fid);
% ---read from DEM file.---
demfid=fopen(fullfile(DEMPathName, DEMFileName), 'r');
data=textscan(demfid, '%s %f', 6);
deminfo=data{1, 2};
formatstr=ones(1, deminfo(1)*3);formatstr=char(formatstr);
for i=1:deminfo(1)
formatstr((i*3-2):(i*3))='%f ';
end
data=textscan(demfid,formatstr);
dem=cell2mat(data);
dem=flipud(dem);
fclose(demfid);
numstem = size(stemline, 1);
StemBase = zeros(numstem, 4);
for count=1:numstem
m=stemline(count,1);n=stemline(count,2);p=stemline(count,3);
x0=stemline(count,4);y0=stemline(count,5);z0=stemline(count,6);
if m==0
if n==0 % m=0 and n=0 means that the line is perpendicular to the xoy plane.
xpix = fix( ( x0 - deminfo(3) ) / deminfo(5) ) + 1;
ypix = fix( ( y0 - deminfo(4) ) / deminfo(5) ) + 1;
if dem(ypix, xpix)==deminfo(6) % the grid has no height value.
% delete this line due to no accurate ground position.
StemBase(count, :) = [ NaN, NaN, NaN, stemline(count, 7) ];
% debug
fprintf('No accurate ground position\n');
% debug
else
StemBase(count, :) = [ x0, y0, dem(ypix, xpix), stemline(count, 7) ];
end
continue;
else
ypix=1:1:deminfo(2);
xpix = ( fix( ( x0 - deminfo(3) ) / deminfo(5) ) + 1 ) * ones(size(ypix));
ymap= (ypix-ones(size(ypix)) + 0.5 )*deminfo(5) + deminfo(4)*ones(size(ypix));
t = (ymap - y0*ones(size(ymap)))/n;
zmap = z0*ones(size(xmap)) + p*t;
tempindex = true(size(ypix));
end
else
if abs(n/m)<=1
xpix=1:1:deminfo(1);
xmap= (xpix-ones(size(xpix)) + 0.5 )*deminfo(5) + deminfo(3)*ones(size(xpix));
t = (xmap - x0*ones(size(xmap)))/m;
ymap = y0*ones(size(xmap)) + n*t;
zmap = z0*ones(size(xmap)) + p*t;
ypix = fix( ( ymap - deminfo(4)*ones(size(ymap)) ) / deminfo(5) ) + ones(size(ymap));
tempindex = ypix >=1 & ypix <= deminfo(2);
else
ypix=1:1:deminfo(2);
ymap= (ypix-ones(size(ypix)) + 0.5 )*deminfo(5) + deminfo(4)*ones(size(ypix));
t = (ymap - y0*ones(size(ymap)))/n;
xmap = x0*ones(size(ymap)) + m*t;
zmap = z0*ones(size(ymap)) + p*t;
xpix = fix( ( xmap - deminfo(3)*ones(size(xmap)) ) / deminfo(5) ) + ones(size(xmap));
tempindex = xpix >=1 & xpix <= deminfo(1);
end
end
zdem = dem(sub2ind([deminfo(2), deminfo(1)], ypix(tempindex), xpix(tempindex)));
validzindex = zdem~=deminfo(6);
zdem = zdem(validzindex);
if length(zdem)<1 % no valid height value
% delete this line due to no accurate ground position.
StemBase(count, :) = [ NaN, NaN, NaN, stemline(count, 7) ];
% debug
fprintf('No accurate ground position\n');
% debug
else
% -----------------------------------------------------------------
% code of original version
% zmap = zmap(tempindex); zmap = zmap(validzindex);
% tempdiff = abs( zdem - zmap );
% tempsub = find( tempdiff==min(tempdiff) );
% if length(tempsub)>1
% tempz = mean(zdem(tempsub));
% t = ( tempz - z0 ) / p;
% StemBase(count, :) = [ x0+m*t, y0+n*t, tempz, stemline(count, 7) ];
% else
% t = ( zdem(tempsub) - z0 ) / p;
% StemBase(count, :) = [ x0+m*t, y0+n*t, zdem(tempsub), stemline(count, 7) ];
% end
% -----------------------------------------------------------------
% -----------------------------------------------------------------
% revision on June 30, 2011, by Zhan Li
zmap = zmap(tempindex); zmap = zmap(validzindex);
tempdiff = zmap - zdem;
% the intersection must be found from those zmap greater than zdem,
% i.e. above the ground rather than below the ground (zmap < zdem)
tempdiff = tempdiff( tempdiff>=0 );
if isempty(tempdiff)
% delete this line due to the whole stem line below the ground.
StemBase(count, :) = [ NaN, NaN, NaN, stemline(count, 7) ];
% debug
fprintf('The whole stem line is below the ground in the DEM extent, possibly due to inaccurate DEM\n');
% debug
else
tempsub = find( tempdiff==min(tempdiff) );
if length(tempsub)>1
tempz = mean(zdem(tempsub));
t = ( tempz - z0 ) / p;
StemBase(count, :) = [ x0+m*t, y0+n*t, tempz, stemline(count, 7) ];
else
t = ( zdem(tempsub) - z0 ) / p;
StemBase(count, :) = [ x0+m*t, y0+n*t, zdem(tempsub), stemline(count, 7) ];
end
end
% -----------------------------------------------------------------
end
end
% delete those rows with NaN, i.e. stem line without finding accurate
% ground position.
tempindex = ~isnan(StemBase(:, 3));
StemBase = StemBase(tempindex, :);
fid=fopen(fullfile(StemBasePathName, StemBaseFileName), 'w');
fprintf(fid, '%s\r\n', 'x, y, z, TreeNO');
fclose(fid);
dlmwrite(fullfile(StemBasePathName, StemBaseFileName), StemBase, '-append', 'delimiter', ',', 'precision', 6, 'newline', 'pc');
end