Rev 83 | Rev 93 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
3 | Werner | 2 | |
3 | #include "tree.h" |
||
83 | Werner | 4 | #include "grid.h" |
30 | Werner | 5 | #include "imagestamp.h" |
3 | Werner | 6 | |
83 | Werner | 7 | #include "stamp.h" |
90 | Werner | 8 | #include "species.h" |
38 | Werner | 9 | |
61 | Werner | 10 | |
39 | Werner | 11 | FloatGrid *Tree::m_grid = 0; |
45 | Werner | 12 | FloatGrid *Tree::m_dominanceGrid = 0; |
40 | Werner | 13 | int Tree::m_statPrint=0; |
48 | Werner | 14 | int Tree::m_statAboveZ=0; |
40 | Werner | 15 | int Tree::m_nextId=0; |
53 | Werner | 16 | float Tree::lafactor = 1.; |
58 | Werner | 17 | int Tree::m_debugid = -1; |
40 | Werner | 18 | |
3 | Werner | 19 | Tree::Tree() |
20 | { |
||
38 | Werner | 21 | m_Dbh = 0; |
22 | m_Height = 0; |
||
23 | m_species = 0; |
||
40 | Werner | 24 | m_id = m_nextId++; |
3 | Werner | 25 | } |
38 | Werner | 26 | |
15 | Werner | 27 | /** get distance and direction between two points. |
38 | Werner | 28 | returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */ |
3 | Werner | 29 | float dist_and_direction(const QPointF &PStart, const QPointF &PEnd, float &rAngle) |
30 | { |
||
31 | float dx = PEnd.x() - PStart.x(); |
||
32 | float dy = PEnd.y() - PStart.y(); |
||
33 | float d = sqrt(dx*dx + dy*dy); |
||
34 | // direction: |
||
35 | rAngle = atan2(dx, dy); |
||
36 | return d; |
||
37 | } |
||
38 | |||
39 | |||
38 | Werner | 40 | void Tree::setup() |
41 | { |
||
42 | if (m_Dbh<=0 || m_Height<=0) |
||
43 | return; |
||
44 | // check stamp |
||
45 | Q_ASSERT_X(m_species!=0, "Tree::setup()", "species is NULL"); |
||
46 | m_stamp = m_species->stamp(m_Dbh, m_Height); |
||
47 | } |
||
39 | Werner | 48 | |
70 | Werner | 49 | #define NOFULLDBG |
77 | Werner | 50 | //#define NOFULLOPT |
51 | /* |
||
39 | Werner | 52 | void Tree::applyStamp() |
53 | { |
||
54 | Q_ASSERT(m_grid!=0); |
||
55 | if (!m_stamp) |
||
56 | return; |
||
57 | |||
40 | Werner | 58 | QPoint pos = m_grid->indexAt(m_Position); |
59 | int offset = m_stamp->offset(); |
||
60 | pos-=QPoint(offset, offset); |
||
61 | QPoint p; |
||
62 | |||
60 | Werner | 63 | float local_dom; // height of Z* on the current position |
40 | Werner | 64 | int x,y; |
58 | Werner | 65 | float value; |
51 | Werner | 66 | QPoint dbg(10,20); |
70 | Werner | 67 | #ifndef NOFULLDBG |
68 | qDebug() <<"indexstampx indexstamy gridx gridy local_dom_height stampvalue 1-value*la/dom grid_before gridvalue_after"; |
||
69 | #endif |
||
40 | Werner | 70 | for (x=0;x<m_stamp->size();++x) { |
71 | for (y=0;y<m_stamp->size(); ++y) { |
||
72 | p = pos + QPoint(x,y); |
||
51 | Werner | 73 | // debug pixel |
61 | Werner | 74 | //if (p==dbg) |
75 | // qDebug() << "#" << id() << "value;"<<(*m_stamp)(x,y)<<"domH"<<dom; |
||
51 | Werner | 76 | |
77 | if (m_grid->isIndexValid(p)) { |
||
78 | // mulitplicative: |
||
53 | Werner | 79 | //m_grid->valueAtIndex(p)*=(*m_stamp)(x,y); |
51 | Werner | 80 | // additiv: |
58 | Werner | 81 | // multiplicatie, v2 |
75 | Werner | 82 | // a optimization for 2m vs 10m grids: local_dom = m_dominanceGrid->valueAtIndex(p.x()/5, p.y()/5); |
83 | // effect is about 20% of the application time |
||
60 | Werner | 84 | local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) ); // todo: could be done more effectively (here 2 transormations are performed)... |
85 | if (local_dom<=0.f) { |
||
61 | Werner | 86 | //qDebug() << "invalid height at " << m_grid->cellCoordinates(p) << "of" << local_dom; |
87 | local_dom = 2.; |
||
60 | Werner | 88 | } |
89 | value = (*m_stamp)(x,y); |
||
90 | value = 1. - value*lafactor / local_dom; |
||
59 | Werner | 91 | value = qMax(value, 0.02f); |
70 | Werner | 92 | #ifndef NOFULLDBG |
93 | qDebug() << x << y << p.x() << p.y() << local_dom << (*m_stamp)(x,y) << 1. - (*m_stamp)(x,y)*lafactor / local_dom << m_grid->valueAtIndex(p) << m_grid->valueAtIndex(p)*value; |
||
94 | #endif |
||
95 | |||
58 | Werner | 96 | m_grid->valueAtIndex(p)*= value; |
51 | Werner | 97 | } |
40 | Werner | 98 | } |
99 | } |
||
58 | Werner | 100 | |
101 | m_statPrint++; // count # of stamp applications... |
||
102 | } |
||
77 | Werner | 103 | */ |
58 | Werner | 104 | |
77 | Werner | 105 | void Tree::applyStamp() |
106 | { |
||
107 | Q_ASSERT(m_grid!=0); |
||
108 | if (!m_stamp) |
||
109 | return; |
||
110 | |||
111 | QPoint pos = m_grid->indexAt(m_Position); |
||
112 | int offset = m_stamp->offset(); |
||
113 | pos-=QPoint(offset, offset); |
||
114 | |||
115 | float local_dom; // height of Z* on the current position |
||
116 | int x,y; |
||
117 | float value; |
||
118 | int gr_stamp = m_stamp->size(); |
||
119 | int grid_x, grid_y; |
||
120 | float *grid_value; |
||
121 | if (!m_grid->isIndexValid(pos) || !m_grid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
||
122 | // todo: in this case we should use another algorithm!!! |
||
123 | return; |
||
124 | } |
||
125 | |||
126 | for (y=0;y<gr_stamp; ++y) { |
||
127 | grid_y = pos.y() + y; |
||
128 | grid_value = m_grid->ptr(pos.x(), grid_y); |
||
129 | for (x=0;x<gr_stamp;++x) { |
||
130 | // suppose there is no stamping outside |
||
131 | grid_x = pos.x() + x; |
||
132 | |||
133 | local_dom = m_dominanceGrid->valueAtIndex(grid_x/5, grid_y/5); |
||
134 | value = (*m_stamp)(x,y); // stampvalue |
||
135 | value = 1. - value*lafactor / local_dom; // calculated value |
||
136 | value = qMax(value, 0.02f); // limit value |
||
137 | |||
138 | *grid_value++ *= value; |
||
139 | } |
||
140 | } |
||
141 | |||
142 | m_statPrint++; // count # of stamp applications... |
||
143 | } |
||
144 | |||
74 | Werner | 145 | /* |
146 | void Tree::heightGrid_old() |
||
58 | Werner | 147 | { |
59 | Werner | 148 | // height of Z* |
74 | Werner | 149 | const float cellsize = m_dominanceGrid->cellsize(); |
62 | Werner | 150 | if (m_Height < cellsize/2.) { |
59 | Werner | 151 | float &dom = m_dominanceGrid->valueAt(m_Position); // modifyable reference |
152 | dom = qMax(dom, m_Height/2.f); |
||
153 | } else { |
||
154 | QPoint p = m_dominanceGrid->indexAt(m_Position); // pos of tree on height grid |
||
62 | Werner | 155 | |
156 | int ringcount = int(floor(m_Height / cellsize)); |
||
59 | Werner | 157 | int ix, iy; |
158 | int ring; |
||
159 | QPoint pos; |
||
62 | Werner | 160 | float hdom; |
161 | float h_out = fmod(m_Height, cellsize) / 2.; |
||
59 | Werner | 162 | for (ix=-ringcount;ix<=ringcount;ix++) |
163 | for (iy=-ringcount; iy<=+ringcount; iy++) { |
||
164 | ring = qMax(abs(ix), abs(iy)); |
||
165 | QPoint pos(ix+p.x(), iy+p.y()); |
||
166 | if (m_dominanceGrid->isIndexValid(pos)) { |
||
167 | // apply value.... |
||
62 | Werner | 168 | if (ring==0) { |
169 | hdom = m_Height- cellsize/4.; |
||
170 | } else if (ring==abs(ringcount)) { |
||
171 | // outermost ring: use again height/2. |
||
172 | hdom = h_out; |
||
173 | } else { |
||
174 | hdom = m_Height- ring*cellsize; |
||
175 | } |
||
176 | m_dominanceGrid->valueAtIndex(pos) = qMax(m_dominanceGrid->valueAtIndex(pos), hdom); |
||
59 | Werner | 177 | } |
45 | Werner | 178 | |
59 | Werner | 179 | } |
180 | } |
||
181 | |||
74 | Werner | 182 | } */ |
183 | |||
184 | /** heightGrid() |
||
185 | This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid. |
||
186 | |||
187 | */ |
||
188 | void Tree::heightGrid() |
||
189 | { |
||
190 | // height of Z* |
||
191 | const float cellsize = m_dominanceGrid->cellsize(); |
||
192 | |||
193 | QPoint p = m_dominanceGrid->indexAt(m_Position); // pos of tree on height grid |
||
194 | QPoint competition_grid = m_grid->indexAt(m_Position); |
||
195 | |||
196 | int index_eastwest = competition_grid.x() % 5; // 4: very west, 0 east edge |
||
197 | int index_northsouth = competition_grid.y() % 5; // 4: northern edge, 0: southern edge |
||
198 | int dist[9]; |
||
199 | dist[3] = index_northsouth * 2 + 1; // south |
||
200 | dist[1] = index_eastwest * 2 + 1; // west |
||
201 | dist[5] = 10 - dist[3]; // north |
||
202 | dist[7] = 10 - dist[1]; // east |
||
203 | dist[8] = qMax(dist[5], dist[7]); // north-east |
||
204 | dist[6] = qMax(dist[3], dist[7]); // south-east |
||
205 | dist[0] = qMax(dist[3], dist[1]); // south-west |
||
206 | dist[2] = qMax(dist[5], dist[1]); // north-west |
||
75 | Werner | 207 | dist[4] = 0; // center cell |
76 | Werner | 208 | /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then: |
209 | index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above. |
||
210 | e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2 |
||
211 | */ |
||
74 | Werner | 212 | |
213 | |||
214 | int ringcount = int(floor(m_Height / cellsize)) + 1; |
||
215 | int ix, iy; |
||
216 | int ring; |
||
217 | QPoint pos; |
||
218 | float hdom; |
||
219 | |||
220 | for (ix=-ringcount;ix<=ringcount;ix++) |
||
221 | for (iy=-ringcount; iy<=+ringcount; iy++) { |
||
222 | ring = qMax(abs(ix), abs(iy)); |
||
223 | QPoint pos(ix+p.x(), iy+p.y()); |
||
224 | if (m_dominanceGrid->isIndexValid(pos)) { |
||
225 | float &rHGrid = m_dominanceGrid->valueAtIndex(pos); |
||
226 | if (rHGrid > m_Height) // skip calculation if grid is higher than tree |
||
227 | continue; |
||
228 | int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y) |
||
229 | hdom = m_Height - dist[direction]; |
||
230 | if (ring>1) |
||
231 | hdom -= (ring-1)*10; |
||
232 | |||
233 | rHGrid = qMax(rHGrid, hdom); // write value |
||
234 | } // is valid |
||
235 | } // for (y) |
||
39 | Werner | 236 | } |
40 | Werner | 237 | |
238 | double Tree::readStamp() |
||
239 | { |
||
51 | Werner | 240 | if (!m_stamp) |
241 | return 0.; |
||
47 | Werner | 242 | const Stamp *stamp = m_stamp->reader(); |
40 | Werner | 243 | if (!stamp) |
244 | return 0.; |
||
245 | QPoint pos = m_grid->indexAt(m_Position); |
||
246 | int offset = stamp->offset(); |
||
247 | pos-=QPoint(offset, offset); |
||
248 | QPoint p; |
||
249 | |||
250 | int x,y; |
||
251 | double sum=0.; |
||
252 | for (x=0;x<stamp->size();++x) { |
||
253 | for (y=0;y<stamp->size(); ++y) { |
||
254 | p = pos + QPoint(x,y); |
||
255 | if (m_grid->isIndexValid(p)) |
||
256 | sum += m_grid->valueAtIndex(p) * (*stamp)(x,y); |
||
257 | } |
||
258 | } |
||
53 | Werner | 259 | float eigenvalue = m_stamp->readSum() * lafactor; |
260 | mImpact = sum - eigenvalue;// additive |
||
261 | float dom_height = (*m_dominanceGrid)[m_Position]; |
||
262 | if (dom_height>0.) |
||
263 | mImpact = mImpact / dom_height; |
||
264 | |||
265 | //mImpact = sum + eigenvalue;// multiplicative |
||
48 | Werner | 266 | // read dominance field... |
53 | Werner | 267 | |
48 | Werner | 268 | if (dom_height < m_Height) { |
269 | // if tree is higher than Z*, the dominance height |
||
270 | // a part of the crown is in "full light". |
||
271 | // total value = zstar/treeheight*value + 1-zstar/height |
||
272 | // reformulated to: |
||
273 | mImpact = mImpact * dom_height/m_Height ; |
||
274 | m_statAboveZ++; |
||
275 | } |
||
47 | Werner | 276 | if (fabs(mImpact < 0.000001)) |
277 | mImpact = 0.f; |
||
278 | qDebug() << "Tree #"<< id() << "value" << sum << "eigenvalue" << eigenvalue << "Impact" << mImpact; |
||
42 | Werner | 279 | return mImpact; |
40 | Werner | 280 | } |
281 | |||
78 | Werner | 282 | /* |
58 | Werner | 283 | double Tree::readStampMul() |
284 | { |
||
285 | if (!m_stamp) |
||
286 | return 0.; |
||
287 | const Stamp *reader = m_stamp->reader(); |
||
288 | if (!reader) |
||
289 | return 0.; |
||
290 | QPoint pos_reader = m_grid->indexAt(m_Position); |
||
291 | |||
292 | int offset_reader = reader->offset(); |
||
293 | int offset_writer = m_stamp->offset(); |
||
294 | int d_offset = offset_writer - offset_reader; |
||
295 | |||
296 | QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer); |
||
297 | pos_reader-=QPoint(offset_reader, offset_reader); |
||
298 | QPoint p; |
||
299 | |||
300 | float dom_height = (*m_dominanceGrid)[m_Position]; |
||
70 | Werner | 301 | float local_dom; |
58 | Werner | 302 | |
303 | int x,y; |
||
304 | double sum=0.; |
||
305 | double value, own_value; |
||
306 | for (x=0;x<reader->size();++x) { |
||
307 | for (y=0;y<reader->size(); ++y) { |
||
308 | p = pos_reader + QPoint(x,y); |
||
309 | if (m_grid->isIndexValid(p)) { |
||
70 | Werner | 310 | local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) ); |
311 | own_value = 1. - m_stamp->offsetValue(x,y,d_offset)*lafactor /local_dom; // old: dom_height; |
||
59 | Werner | 312 | own_value = qMax(own_value, 0.02); |
58 | Werner | 313 | value = m_grid->valueAtIndex(p) / own_value; // remove impact of focal tree |
314 | if (value>0.) |
||
315 | sum += value * (*reader)(x,y); |
||
316 | //value = (1. - m_stamp->offsetValue(x,y,d_offset)/dom_height); |
||
317 | //value = (1 - m_grid->valueAtIndex(p)) / (m_stamp->offsetValue(x,y,d_offset)/dom_height * lafactor); |
||
318 | if (isDebugging()) { |
||
319 | qDebug() << "Tree" << id() << "Coord" << p << "gridvalue" << m_grid->valueAtIndex(p) << "value" << value << "reader" << (*reader)(x,y) << "writer" << m_stamp->offsetValue(x,y,d_offset); |
||
320 | } |
||
321 | |||
322 | } |
||
323 | } |
||
324 | } |
||
325 | mImpact = sum; |
||
326 | // read dominance field... |
||
327 | if (dom_height < m_Height) { |
||
328 | // if tree is higher than Z*, the dominance height |
||
329 | // a part of the crown is in "full light". |
||
330 | m_statAboveZ++; |
||
331 | mImpact = 1. - (1. - mImpact)*dom_height/m_Height; |
||
332 | } |
||
333 | if (mImpact > 1) |
||
334 | mImpact = 1.f; |
||
61 | Werner | 335 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact; |
58 | Werner | 336 | return mImpact; |
78 | Werner | 337 | }*/ |
338 | |||
339 | double Tree::readStampMul() |
||
340 | { |
||
341 | if (!m_stamp) |
||
342 | return 0.; |
||
343 | const Stamp *reader = m_stamp->reader(); |
||
344 | if (!reader) |
||
345 | return 0.; |
||
346 | QPoint pos_reader = m_grid->indexAt(m_Position); |
||
347 | |||
348 | int offset_reader = reader->offset(); |
||
349 | int offset_writer = m_stamp->offset(); |
||
350 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
||
351 | |||
352 | QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer); |
||
353 | pos_reader-=QPoint(offset_reader, offset_reader); |
||
354 | QPoint p; |
||
355 | |||
356 | //float dom_height = (*m_dominanceGrid)[m_Position]; |
||
357 | float local_dom; |
||
358 | |||
359 | int x,y; |
||
360 | double sum=0.; |
||
361 | double value, own_value; |
||
362 | float *grid_value; |
||
363 | int reader_size = reader->size(); |
||
364 | int rx = pos_reader.x(); |
||
365 | int ry = pos_reader.y(); |
||
366 | for (y=0;y<reader_size; ++y, ++ry) { |
||
367 | grid_value = m_grid->ptr(rx, ry); |
||
368 | for (x=0;x<reader_size;++x) { |
||
369 | |||
370 | //p = pos_reader + QPoint(x,y); |
||
371 | //if (m_grid->isIndexValid(p)) { |
||
372 | local_dom = m_dominanceGrid->valueAtIndex((rx+x)/5, ry/5); |
||
373 | //local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) ); |
||
374 | own_value = 1. - m_stamp->offsetValue(x,y,d_offset)*lafactor /local_dom; // old: dom_height; |
||
375 | own_value = qMax(own_value, 0.02); |
||
376 | value = *grid_value++ / own_value; // remove impact of focal tree |
||
377 | //if (value>0.) |
||
378 | sum += value * (*reader)(x,y); |
||
379 | |||
380 | //} // isIndexValid |
||
381 | } |
||
382 | } |
||
383 | mImpact = sum; |
||
384 | // read dominance field... |
||
385 | // this applies only if some trees are potentially *higher* than the dominant height grid |
||
386 | //if (dom_height < m_Height) { |
||
387 | // if tree is higher than Z*, the dominance height |
||
388 | // a part of the crown is in "full light". |
||
389 | // m_statAboveZ++; |
||
390 | // mImpact = 1. - (1. - mImpact)*dom_height/m_Height; |
||
391 | //} |
||
392 | if (mImpact > 1) |
||
393 | mImpact = 1.f; |
||
394 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact; |
||
395 | return mImpact; |
||
58 | Werner | 396 | } |
397 | |||
40 | Werner | 398 | void Tree::resetStatistics() |
399 | { |
||
400 | m_statPrint=0; |
||
48 | Werner | 401 | m_statAboveZ=0; |
40 | Werner | 402 | m_nextId=1; |
403 | } |