Rev 70 | Rev 75 | 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" |
||
16 | Werner | 4 | #include "core/grid.h" |
30 | Werner | 5 | #include "imagestamp.h" |
3 | Werner | 6 | |
38 | Werner | 7 | #include "core/stamp.h" |
8 | #include "treespecies.h" |
||
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 |
39 | Werner | 50 | void Tree::applyStamp() |
51 | { |
||
52 | Q_ASSERT(m_grid!=0); |
||
53 | if (!m_stamp) |
||
54 | return; |
||
55 | |||
40 | Werner | 56 | QPoint pos = m_grid->indexAt(m_Position); |
57 | int offset = m_stamp->offset(); |
||
58 | pos-=QPoint(offset, offset); |
||
59 | QPoint p; |
||
60 | |||
60 | Werner | 61 | float local_dom; // height of Z* on the current position |
40 | Werner | 62 | int x,y; |
58 | Werner | 63 | float value; |
51 | Werner | 64 | QPoint dbg(10,20); |
70 | Werner | 65 | #ifndef NOFULLDBG |
66 | qDebug() <<"indexstampx indexstamy gridx gridy local_dom_height stampvalue 1-value*la/dom grid_before gridvalue_after"; |
||
67 | #endif |
||
40 | Werner | 68 | for (x=0;x<m_stamp->size();++x) { |
69 | for (y=0;y<m_stamp->size(); ++y) { |
||
70 | p = pos + QPoint(x,y); |
||
51 | Werner | 71 | // debug pixel |
61 | Werner | 72 | //if (p==dbg) |
73 | // qDebug() << "#" << id() << "value;"<<(*m_stamp)(x,y)<<"domH"<<dom; |
||
51 | Werner | 74 | |
75 | if (m_grid->isIndexValid(p)) { |
||
76 | // mulitplicative: |
||
53 | Werner | 77 | //m_grid->valueAtIndex(p)*=(*m_stamp)(x,y); |
51 | Werner | 78 | // additiv: |
58 | Werner | 79 | // multiplicatie, v2 |
60 | Werner | 80 | local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) ); // todo: could be done more effectively (here 2 transormations are performed)... |
81 | if (local_dom<=0.f) { |
||
61 | Werner | 82 | //qDebug() << "invalid height at " << m_grid->cellCoordinates(p) << "of" << local_dom; |
83 | local_dom = 2.; |
||
60 | Werner | 84 | } |
85 | value = (*m_stamp)(x,y); |
||
86 | value = 1. - value*lafactor / local_dom; |
||
59 | Werner | 87 | value = qMax(value, 0.02f); |
70 | Werner | 88 | #ifndef NOFULLDBG |
89 | 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; |
||
90 | #endif |
||
91 | |||
58 | Werner | 92 | m_grid->valueAtIndex(p)*= value; |
51 | Werner | 93 | } |
40 | Werner | 94 | } |
95 | } |
||
58 | Werner | 96 | |
97 | m_statPrint++; // count # of stamp applications... |
||
98 | } |
||
99 | |||
74 | Werner | 100 | /* |
101 | void Tree::heightGrid_old() |
||
58 | Werner | 102 | { |
59 | Werner | 103 | // height of Z* |
74 | Werner | 104 | const float cellsize = m_dominanceGrid->cellsize(); |
62 | Werner | 105 | if (m_Height < cellsize/2.) { |
59 | Werner | 106 | float &dom = m_dominanceGrid->valueAt(m_Position); // modifyable reference |
107 | dom = qMax(dom, m_Height/2.f); |
||
108 | } else { |
||
109 | QPoint p = m_dominanceGrid->indexAt(m_Position); // pos of tree on height grid |
||
62 | Werner | 110 | |
111 | int ringcount = int(floor(m_Height / cellsize)); |
||
59 | Werner | 112 | int ix, iy; |
113 | int ring; |
||
114 | QPoint pos; |
||
62 | Werner | 115 | float hdom; |
116 | float h_out = fmod(m_Height, cellsize) / 2.; |
||
59 | Werner | 117 | for (ix=-ringcount;ix<=ringcount;ix++) |
118 | for (iy=-ringcount; iy<=+ringcount; iy++) { |
||
119 | ring = qMax(abs(ix), abs(iy)); |
||
120 | QPoint pos(ix+p.x(), iy+p.y()); |
||
121 | if (m_dominanceGrid->isIndexValid(pos)) { |
||
122 | // apply value.... |
||
62 | Werner | 123 | if (ring==0) { |
124 | hdom = m_Height- cellsize/4.; |
||
125 | } else if (ring==abs(ringcount)) { |
||
126 | // outermost ring: use again height/2. |
||
127 | hdom = h_out; |
||
128 | } else { |
||
129 | hdom = m_Height- ring*cellsize; |
||
130 | } |
||
131 | m_dominanceGrid->valueAtIndex(pos) = qMax(m_dominanceGrid->valueAtIndex(pos), hdom); |
||
59 | Werner | 132 | } |
45 | Werner | 133 | |
59 | Werner | 134 | } |
135 | } |
||
136 | |||
74 | Werner | 137 | } */ |
138 | |||
139 | /** heightGrid() |
||
140 | This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid. |
||
141 | |||
142 | */ |
||
143 | void Tree::heightGrid() |
||
144 | { |
||
145 | // height of Z* |
||
146 | const float cellsize = m_dominanceGrid->cellsize(); |
||
147 | |||
148 | QPoint p = m_dominanceGrid->indexAt(m_Position); // pos of tree on height grid |
||
149 | QPoint competition_grid = m_grid->indexAt(m_Position); |
||
150 | |||
151 | int index_eastwest = competition_grid.x() % 5; // 4: very west, 0 east edge |
||
152 | int index_northsouth = competition_grid.y() % 5; // 4: northern edge, 0: southern edge |
||
153 | int dist[9]; |
||
154 | dist[3] = index_northsouth * 2 + 1; // south |
||
155 | dist[1] = index_eastwest * 2 + 1; // west |
||
156 | dist[5] = 10 - dist[3]; // north |
||
157 | dist[7] = 10 - dist[1]; // east |
||
158 | dist[8] = qMax(dist[5], dist[7]); // north-east |
||
159 | dist[6] = qMax(dist[3], dist[7]); // south-east |
||
160 | dist[0] = qMax(dist[3], dist[1]); // south-west |
||
161 | dist[2] = qMax(dist[5], dist[1]); // north-west |
||
162 | dist[4] = 0; |
||
163 | |||
164 | |||
165 | int ringcount = int(floor(m_Height / cellsize)) + 1; |
||
166 | int ix, iy; |
||
167 | int ring; |
||
168 | QPoint pos; |
||
169 | float hdom; |
||
170 | |||
171 | for (ix=-ringcount;ix<=ringcount;ix++) |
||
172 | for (iy=-ringcount; iy<=+ringcount; iy++) { |
||
173 | ring = qMax(abs(ix), abs(iy)); |
||
174 | QPoint pos(ix+p.x(), iy+p.y()); |
||
175 | if (m_dominanceGrid->isIndexValid(pos)) { |
||
176 | float &rHGrid = m_dominanceGrid->valueAtIndex(pos); |
||
177 | if (rHGrid > m_Height) // skip calculation if grid is higher than tree |
||
178 | continue; |
||
179 | int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y) |
||
180 | hdom = m_Height - dist[direction]; |
||
181 | if (ring>1) |
||
182 | hdom -= (ring-1)*10; |
||
183 | |||
184 | rHGrid = qMax(rHGrid, hdom); // write value |
||
185 | } // is valid |
||
186 | } // for (y) |
||
39 | Werner | 187 | } |
40 | Werner | 188 | |
189 | double Tree::readStamp() |
||
190 | { |
||
51 | Werner | 191 | if (!m_stamp) |
192 | return 0.; |
||
47 | Werner | 193 | const Stamp *stamp = m_stamp->reader(); |
40 | Werner | 194 | if (!stamp) |
195 | return 0.; |
||
196 | QPoint pos = m_grid->indexAt(m_Position); |
||
197 | int offset = stamp->offset(); |
||
198 | pos-=QPoint(offset, offset); |
||
199 | QPoint p; |
||
200 | |||
201 | int x,y; |
||
202 | double sum=0.; |
||
203 | for (x=0;x<stamp->size();++x) { |
||
204 | for (y=0;y<stamp->size(); ++y) { |
||
205 | p = pos + QPoint(x,y); |
||
206 | if (m_grid->isIndexValid(p)) |
||
207 | sum += m_grid->valueAtIndex(p) * (*stamp)(x,y); |
||
208 | } |
||
209 | } |
||
53 | Werner | 210 | float eigenvalue = m_stamp->readSum() * lafactor; |
211 | mImpact = sum - eigenvalue;// additive |
||
212 | float dom_height = (*m_dominanceGrid)[m_Position]; |
||
213 | if (dom_height>0.) |
||
214 | mImpact = mImpact / dom_height; |
||
215 | |||
216 | //mImpact = sum + eigenvalue;// multiplicative |
||
48 | Werner | 217 | // read dominance field... |
53 | Werner | 218 | |
48 | Werner | 219 | if (dom_height < m_Height) { |
220 | // if tree is higher than Z*, the dominance height |
||
221 | // a part of the crown is in "full light". |
||
222 | // total value = zstar/treeheight*value + 1-zstar/height |
||
223 | // reformulated to: |
||
224 | mImpact = mImpact * dom_height/m_Height ; |
||
225 | m_statAboveZ++; |
||
226 | } |
||
47 | Werner | 227 | if (fabs(mImpact < 0.000001)) |
228 | mImpact = 0.f; |
||
229 | qDebug() << "Tree #"<< id() << "value" << sum << "eigenvalue" << eigenvalue << "Impact" << mImpact; |
||
42 | Werner | 230 | return mImpact; |
40 | Werner | 231 | } |
232 | |||
58 | Werner | 233 | double Tree::readStampMul() |
234 | { |
||
235 | if (!m_stamp) |
||
236 | return 0.; |
||
237 | const Stamp *reader = m_stamp->reader(); |
||
238 | if (!reader) |
||
239 | return 0.; |
||
240 | QPoint pos_reader = m_grid->indexAt(m_Position); |
||
241 | |||
242 | int offset_reader = reader->offset(); |
||
243 | int offset_writer = m_stamp->offset(); |
||
244 | int d_offset = offset_writer - offset_reader; |
||
245 | |||
246 | QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer); |
||
247 | pos_reader-=QPoint(offset_reader, offset_reader); |
||
248 | QPoint p; |
||
249 | |||
250 | float dom_height = (*m_dominanceGrid)[m_Position]; |
||
70 | Werner | 251 | float local_dom; |
58 | Werner | 252 | |
253 | int x,y; |
||
254 | double sum=0.; |
||
255 | double value, own_value; |
||
256 | for (x=0;x<reader->size();++x) { |
||
257 | for (y=0;y<reader->size(); ++y) { |
||
258 | p = pos_reader + QPoint(x,y); |
||
259 | if (m_grid->isIndexValid(p)) { |
||
70 | Werner | 260 | local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) ); |
261 | own_value = 1. - m_stamp->offsetValue(x,y,d_offset)*lafactor /local_dom; // old: dom_height; |
||
59 | Werner | 262 | own_value = qMax(own_value, 0.02); |
58 | Werner | 263 | value = m_grid->valueAtIndex(p) / own_value; // remove impact of focal tree |
264 | if (value>0.) |
||
265 | sum += value * (*reader)(x,y); |
||
266 | //value = (1. - m_stamp->offsetValue(x,y,d_offset)/dom_height); |
||
267 | //value = (1 - m_grid->valueAtIndex(p)) / (m_stamp->offsetValue(x,y,d_offset)/dom_height * lafactor); |
||
268 | if (isDebugging()) { |
||
269 | qDebug() << "Tree" << id() << "Coord" << p << "gridvalue" << m_grid->valueAtIndex(p) << "value" << value << "reader" << (*reader)(x,y) << "writer" << m_stamp->offsetValue(x,y,d_offset); |
||
270 | } |
||
271 | |||
272 | } |
||
273 | } |
||
274 | } |
||
275 | mImpact = sum; |
||
276 | // read dominance field... |
||
277 | if (dom_height < m_Height) { |
||
278 | // if tree is higher than Z*, the dominance height |
||
279 | // a part of the crown is in "full light". |
||
280 | m_statAboveZ++; |
||
281 | mImpact = 1. - (1. - mImpact)*dom_height/m_Height; |
||
282 | } |
||
283 | if (mImpact > 1) |
||
284 | mImpact = 1.f; |
||
61 | Werner | 285 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact; |
58 | Werner | 286 | return mImpact; |
287 | } |
||
288 | |||
289 | |||
40 | Werner | 290 | void Tree::resetStatistics() |
291 | { |
||
292 | m_statPrint=0; |
||
48 | Werner | 293 | m_statAboveZ=0; |
40 | Werner | 294 | m_nextId=1; |
295 | } |