28 #include <stack> |
28 #include <stack> |
29 |
29 |
30 #include "common.h" |
30 #include "common.h" |
31 #include "kdtree.h" |
31 #include "kdtree.h" |
32 |
32 |
33 class SortableShape |
33 class ShapeBound |
34 { |
34 { |
35 public: |
35 public: |
36 Shape *shape; |
36 Shape *shape; |
37 BBox bbox; |
37 Float pos; |
38 short axis; |
38 bool end; |
39 short mark; |
39 ShapeBound(Shape *ashape, const Float apos, const bool aend): |
40 |
40 shape(ashape), pos(apos), end(aend) {}; |
41 SortableShape(Shape *aShape, short &aAxis): shape(aShape), axis(aAxis), mark(0) |
41 friend bool operator<(const ShapeBound& a, const ShapeBound& b) |
42 { bbox = shape->get_bbox(); }; |
42 { |
43 friend bool operator<(const SortableShape& a, const SortableShape& b) |
43 if (a.pos == b.pos) |
44 { return a.bbox.L.cell[a.axis] < b.bbox.L.cell[b.axis]; }; |
44 return a.shape < b.shape; |
45 void setAxis(short aAxis) { axis = aAxis; }; |
45 else |
46 void setMark() { mark = 1; }; |
46 return a.pos < b.pos; |
47 short hasMark() { return mark; }; |
|
48 }; |
|
49 |
|
50 class SortableShapeList: public vector<SortableShape> |
|
51 { |
|
52 public: |
|
53 SortableShapeList(ShapeList *shapes, short axis) |
|
54 { |
|
55 ShapeList::iterator shape; |
|
56 for (shape = shapes->begin(); shape != shapes->end(); shape++) |
|
57 push_back(SortableShape(*shape, axis)); |
|
58 }; |
47 }; |
59 }; |
|
60 |
|
61 class SplitPos |
|
62 { |
|
63 public: |
|
64 Float pos; |
|
65 int lnum, rnum; |
|
66 SplitPos(): pos(0.0), lnum(0), rnum(0) {}; |
|
67 SplitPos(Float &aPos): pos(aPos), lnum(0), rnum(0) {}; |
|
68 friend bool operator<(const SplitPos& a, const SplitPos& b) |
|
69 { return a.pos < b.pos; }; |
|
70 friend bool operator==(const SplitPos& a, const SplitPos& b) |
|
71 { return a.pos == b.pos; }; |
|
72 }; |
|
73 |
|
74 class SplitList: public vector<SplitPos> |
|
75 { |
|
76 }; |
48 }; |
77 |
49 |
78 // stack element for kd-tree traversal |
50 // stack element for kd-tree traversal |
79 class StackElem |
51 class StackElem |
80 { |
52 { |
94 delete shapes; |
66 delete shapes; |
95 else |
67 else |
96 delete[] children; |
68 delete[] children; |
97 } |
69 } |
98 |
70 |
99 void KdNode::subdivide(BBox bbox, int maxdepth) |
71 // kd-tree recursive build algorithm, inspired by PBRT (www.pbrt.org) |
|
72 void KdNode::subdivide(BBox bounds, int maxdepth) |
100 { |
73 { |
101 if (maxdepth <= 0 || shapes->size() <= 2) |
74 if (maxdepth <= 0 || shapes->size() <= 2) |
102 { |
75 { |
103 setLeaf(); |
76 setLeaf(); |
104 return; |
77 return; |
105 } |
78 } |
106 |
79 |
107 // choose split axis |
80 // choose split axis |
108 axis = 0; |
81 axis = 0; |
109 if (bbox.h() > bbox.w() && bbox.h() > bbox.d()) |
82 if (bounds.h() > bounds.w() && bounds.h() > bounds.d()) |
110 axis = 1; |
83 axis = 1; |
111 if (bbox.d() > bbox.w() && bbox.d() > bbox.h()) |
84 if (bounds.d() > bounds.w() && bounds.d() > bounds.h()) |
112 axis = 2; |
85 axis = 2; |
113 |
86 |
114 // *** find optimal split position |
87 // create sorted list of shape bounds (= find all posible splits) |
115 SortableShapeList sslist(shapes, axis); |
88 vector<ShapeBound> edges[3]; |
116 sort(sslist.begin(), sslist.end()); |
89 ShapeList::iterator shape; |
117 |
90 for (shape = shapes->begin(); shape != shapes->end(); shape++) |
118 SplitList splitlist; |
91 { |
119 |
92 BBox shapebounds = (*shape)->get_bbox(); |
120 SortableShapeList::iterator sh; |
93 // for (int ax = 0; ax < 3; ax++) |
121 for (sh = sslist.begin(); sh != sslist.end(); sh++) |
94 // { |
122 { |
95 edges[axis].push_back(ShapeBound(*shape, shapebounds.L[axis], 0)); |
123 splitlist.push_back(SplitPos(sh->bbox.L.cell[axis])); |
96 edges[axis].push_back(ShapeBound(*shape, shapebounds.H[axis], 1)); |
124 splitlist.push_back(SplitPos(sh->bbox.H.cell[axis])); |
97 // } |
125 } |
98 } |
126 sort(splitlist.begin(), splitlist.end()); |
99 sort(edges[axis].begin(), edges[axis].end()); |
127 // remove duplicate splits |
|
128 splitlist.resize(unique(splitlist.begin(), splitlist.end()) - splitlist.begin()); |
|
129 |
|
130 // find all posible splits |
|
131 SplitList::iterator spl; |
|
132 int rest; |
|
133 for (spl = splitlist.begin(); spl != splitlist.end(); spl++) |
|
134 { |
|
135 for (sh = sslist.begin(), rest = sslist.size(); sh != sslist.end(); sh++, rest--) |
|
136 { |
|
137 if (sh->hasMark()) |
|
138 { |
|
139 spl->lnum++; |
|
140 continue; |
|
141 } |
|
142 // if shape is completely contained in split plane |
|
143 if (spl->pos == sh->bbox.L.cell[axis] == sh->bbox.H.cell[axis]) |
|
144 { |
|
145 if (spl->pos - bbox.L.cell[axis] < bbox.H.cell[axis] - spl->pos) |
|
146 { |
|
147 // left subcell is smaller -> if not empty, put shape here |
|
148 if (spl->lnum) |
|
149 spl->lnum++; |
|
150 else |
|
151 spl->rnum++; |
|
152 } else { |
|
153 // right subcell is smaller |
|
154 if (spl->rnum) |
|
155 spl->rnum++; |
|
156 else |
|
157 spl->lnum++; |
|
158 } |
|
159 sh->setMark(); |
|
160 } else |
|
161 // if shape is on left side of split plane |
|
162 if (sh->bbox.H.cell[axis] <= spl->pos) |
|
163 { |
|
164 spl->lnum++; |
|
165 sh->setMark(); |
|
166 } else |
|
167 // if shape occupies both sides of split plane |
|
168 if (sh->bbox.L.cell[axis] < spl->pos && sh->bbox.H.cell[axis] > spl->pos) |
|
169 { |
|
170 spl->lnum++; |
|
171 spl->rnum++; |
|
172 } else |
|
173 // if shape is on right side of split plane |
|
174 if (sh->bbox.L.cell[axis] >= spl->pos) |
|
175 { |
|
176 spl->rnum += rest; |
|
177 break; |
|
178 } |
|
179 } |
|
180 } |
|
181 |
100 |
182 // choose best split pos |
101 // choose best split pos |
183 const Float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection |
102 const Float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection |
184 Float SAV = (bbox.w()*bbox.h() + bbox.w()*bbox.d() + bbox.h()*bbox.d()); // surface area of node |
103 Float SAV = (bounds.w()*bounds.h() + // surface area of node |
|
104 bounds.w()*bounds.d() + bounds.h()*bounds.d()); |
185 Float cost = SAV * (K + shapes->size()); // initial cost = non-split cost |
105 Float cost = SAV * (K + shapes->size()); // initial cost = non-split cost |
186 bool leaf = true; |
106 BBox lbb = bounds; |
187 BBox lbb = bbox; |
107 BBox rbb = bounds; |
188 BBox rbb = bbox; |
108 |
189 for (spl = splitlist.begin(); spl != splitlist.end(); spl++) |
109 vector<ShapeBound>::iterator edge, splitedge = edges[axis].end(); |
190 { |
110 int lnum = 0, rnum = shapes->size(); |
|
111 for (edge = edges[axis].begin(); edge != edges[axis].end(); edge++) |
|
112 { |
|
113 if (edge->end) |
|
114 rnum--; |
|
115 |
191 // calculate SAH cost of this split |
116 // calculate SAH cost of this split |
192 lbb.H.cell[axis] = spl->pos; |
117 lbb.H.cell[axis] = edge->pos; |
193 rbb.L.cell[axis] = spl->pos; |
118 rbb.L.cell[axis] = edge->pos; |
194 Float SAL = (lbb.w()*lbb.h() + lbb.w()*lbb.d() + lbb.h()*lbb.d()); |
119 Float SAL = (lbb.w()*lbb.h() + lbb.w()*lbb.d() + lbb.h()*lbb.d()); |
195 Float SAR = (rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d()); |
120 Float SAR = (rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d()); |
196 Float splitcost = K*SAV + SAL*(K+spl->lnum) + SAR*(K+spl->rnum); |
121 Float splitcost = K*SAV + SAL*(K + lnum) + SAR*(K + rnum); |
197 |
122 |
198 if (splitcost < cost) |
123 if (splitcost < cost) |
199 { |
124 { |
200 leaf = false; |
125 splitedge = edge; |
201 cost = splitcost; |
126 cost = splitcost; |
202 split = spl->pos; |
127 split = edge->pos; |
203 } |
128 } |
204 } |
129 |
205 |
130 if (!edge->end) |
206 if (leaf) |
131 lnum++; |
|
132 } |
|
133 |
|
134 if (splitedge == edges[axis].end()) |
207 { |
135 { |
208 setLeaf(); |
136 setLeaf(); |
209 return; |
137 return; |
210 } |
138 } |
211 |
139 |
214 // this would be hard to reconstruct later |
142 // this would be hard to reconstruct later |
215 static ofstream F("kdtree.obj"); |
143 static ofstream F("kdtree.obj"); |
216 Vector3 v; |
144 Vector3 v; |
217 static int f=0; |
145 static int f=0; |
218 v.cell[axis] = split; |
146 v.cell[axis] = split; |
219 v.cell[(axis+1)%3] = bbox.L.cell[(axis+1)%3]; |
147 v.cell[(axis+1)%3] = bounds.L.cell[(axis+1)%3]; |
220 v.cell[(axis+2)%3] = bbox.L.cell[(axis+2)%3]; |
148 v.cell[(axis+2)%3] = bounds.L.cell[(axis+2)%3]; |
221 F << "v " << v.x << " " << v.y << " " << v.z << endl; |
149 F << "v " << v.x << " " << v.y << " " << v.z << endl; |
222 v.cell[(axis+1)%3] = bbox.L.cell[(axis+1)%3]; |
150 v.cell[(axis+1)%3] = bounds.L.cell[(axis+1)%3]; |
223 v.cell[(axis+2)%3] = bbox.H.cell[(axis+2)%3]; |
151 v.cell[(axis+2)%3] = bounds.H.cell[(axis+2)%3]; |
224 F << "v " << v.x << " " << v.y << " " << v.z << endl; |
152 F << "v " << v.x << " " << v.y << " " << v.z << endl; |
225 v.cell[(axis+1)%3] = bbox.H.cell[(axis+1)%3]; |
153 v.cell[(axis+1)%3] = bounds.H.cell[(axis+1)%3]; |
226 v.cell[(axis+2)%3] = bbox.H.cell[(axis+2)%3]; |
154 v.cell[(axis+2)%3] = bounds.H.cell[(axis+2)%3]; |
227 F << "v " << v.x << " " << v.y << " " << v.z << endl; |
155 F << "v " << v.x << " " << v.y << " " << v.z << endl; |
228 v.cell[(axis+1)%3] = bbox.H.cell[(axis+1)%3]; |
156 v.cell[(axis+1)%3] = bounds.H.cell[(axis+1)%3]; |
229 v.cell[(axis+2)%3] = bbox.L.cell[(axis+2)%3]; |
157 v.cell[(axis+2)%3] = bounds.L.cell[(axis+2)%3]; |
230 F << "v " << v.x << " " << v.y << " " << v.z << endl; |
158 F << "v " << v.x << " " << v.y << " " << v.z << endl; |
231 F << "f " << f+1 << " " << f+2 << " " << f+3 << " " << f+4 << endl; |
159 F << "f " << f+1 << " " << f+2 << " " << f+3 << " " << f+4 << endl; |
232 f += 4; |
160 f += 4; |
233 #endif |
161 #endif |
234 |
162 |
235 // split this node |
163 // split this node |
236 delete shapes; |
164 delete shapes; |
237 children = new KdNode[2]; |
165 children = new KdNode[2]; |
238 int state = 0; |
166 for (edge = edges[axis].begin(); edge != splitedge; edge++) |
239 for (sh = sslist.begin(); sh != sslist.end(); sh++) |
167 if (!edge->end) |
240 { |
168 children[0].addShape(edge->shape); |
241 // if shape is on left side of split plane |
169 for (edge = splitedge; edge < edges[axis].end(); edge++) |
242 if (state == 1) |
170 if (edge->end) |
243 { // only right |
171 children[1].addShape(edge->shape); |
244 children[1].addShape(sh->shape); |
|
245 continue; |
|
246 } |
|
247 if (state == 0) |
|
248 { |
|
249 if (sh->bbox.H.cell[axis] < split) |
|
250 { // left |
|
251 children[0].addShape(sh->shape); |
|
252 } else |
|
253 if (sh->bbox.H.cell[axis] > split) |
|
254 { |
|
255 if (sh->bbox.L.cell[axis] < split) |
|
256 { // both |
|
257 children[0].addShape(sh->shape); |
|
258 children[1].addShape(sh->shape); |
|
259 } else |
|
260 { // right |
|
261 children[1].addShape(sh->shape); |
|
262 state = 1; |
|
263 } |
|
264 } else |
|
265 { // H == split |
|
266 if (sh->bbox.L.cell[axis] < split) |
|
267 { // left |
|
268 children[0].addShape(sh->shape); |
|
269 } else |
|
270 { // contained |
|
271 if (split - bbox.L.cell[axis] < bbox.H.cell[axis] - split) |
|
272 { |
|
273 // left subcell is smaller -> if not empty, put shape here |
|
274 if (children[0].shapes->size()) |
|
275 children[0].addShape(sh->shape); |
|
276 else |
|
277 children[1].addShape(sh->shape); |
|
278 } else { |
|
279 // right subcell is smaller |
|
280 if (children[1].shapes->size()) |
|
281 children[1].addShape(sh->shape); |
|
282 else |
|
283 children[0].addShape(sh->shape); |
|
284 } |
|
285 } |
|
286 } |
|
287 } |
|
288 } |
|
289 |
172 |
290 lbb.H.cell[axis] = split; |
173 lbb.H.cell[axis] = split; |
291 rbb.L.cell[axis] = split; |
174 rbb.L.cell[axis] = split; |
292 |
175 |
293 children[0].subdivide(lbb, maxdepth-1); |
176 children[0].subdivide(lbb, maxdepth-1); |