65 else |
65 else |
66 delete[] getLeftChild(); |
66 delete[] getLeftChild(); |
67 } |
67 } |
68 |
68 |
69 // kd-tree recursive build algorithm, inspired by PBRT (www.pbrt.org) |
69 // kd-tree recursive build algorithm, inspired by PBRT (www.pbrt.org) |
70 void KdNode::subdivide(BBox bounds, int maxdepth) |
70 void KdTree::recursive_build(KdNode *node, BBox bounds, int maxdepth) |
71 { |
71 { |
72 if (maxdepth <= 0 || getShapes()->size() <= 2) |
72 ShapeList *shapes = node->getShapes(); |
73 { |
73 |
74 setLeaf(); |
74 if (maxdepth <= 0 || shapes->size() <= 2) |
|
75 { |
|
76 node->setLeaf(); |
75 return; |
77 return; |
76 } |
78 } |
77 |
79 |
78 // choose split axis |
80 // choose split axis |
79 /*axis = 0; |
81 /*axis = 0; |
83 axis = 2; |
85 axis = 2; |
84 */ |
86 */ |
85 // create sorted list of shape bounds (= find all posible splits) |
87 // create sorted list of shape bounds (= find all posible splits) |
86 vector<ShapeBound> edges[3]; |
88 vector<ShapeBound> edges[3]; |
87 ShapeList::iterator shape; |
89 ShapeList::iterator shape; |
88 for (shape = getShapes()->begin(); shape != getShapes()->end(); shape++) |
90 for (shape = shapes->begin(); shape != shapes->end(); shape++) |
89 { |
91 { |
90 BBox shapebounds = (*shape)->get_bbox(); |
92 BBox shapebounds = (*shape)->get_bbox(); |
91 for (int ax = 0; ax < 3; ax++) |
93 for (int ax = 0; ax < 3; ax++) |
92 { |
94 { |
93 edges[ax].push_back(ShapeBound(*shape, shapebounds.L[ax], 0)); |
95 edges[ax].push_back(ShapeBound(*shape, shapebounds.L[ax], 0)); |
99 |
101 |
100 // choose best split pos |
102 // choose best split pos |
101 const Float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection |
103 const Float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection |
102 Float SAV = (bounds.w()*bounds.h() + // surface area of node |
104 Float SAV = (bounds.w()*bounds.h() + // surface area of node |
103 bounds.w()*bounds.d() + bounds.h()*bounds.d()); |
105 bounds.w()*bounds.d() + bounds.h()*bounds.d()); |
104 Float cost = SAV * (K + getShapes()->size()); // initial cost = non-split cost |
106 Float cost = SAV * (K + shapes->size()); // initial cost = non-split cost |
105 |
107 |
106 vector<ShapeBound>::iterator edge, splitedge = edges[2].end(); |
108 vector<ShapeBound>::iterator edge, splitedge = edges[2].end(); |
107 int axis = 0; |
109 int axis = 0; |
108 for (int ax = 0; ax < 3; ax++) |
110 for (int ax = 0; ax < 3; ax++) |
109 { |
111 { |
110 int lnum = 0, rnum = getShapes()->size(); |
112 int lnum = 0, rnum = shapes->size(); |
111 BBox lbb = bounds; |
113 BBox lbb = bounds; |
112 BBox rbb = bounds; |
114 BBox rbb = bounds; |
113 for (edge = edges[ax].begin(); edge != edges[ax].end(); edge++) |
115 for (edge = edges[ax].begin(); edge != edges[ax].end(); edge++) |
114 { |
116 { |
115 if (edge->end) |
117 if (edge->end) |
125 if (splitcost < cost) |
127 if (splitcost < cost) |
126 { |
128 { |
127 axis = ax; |
129 axis = ax; |
128 splitedge = edge; |
130 splitedge = edge; |
129 cost = splitcost; |
131 cost = splitcost; |
130 split = edge->pos; |
|
131 } |
132 } |
132 |
133 |
133 if (!edge->end) |
134 if (!edge->end) |
134 lnum++; |
135 lnum++; |
135 } |
136 } |
136 } |
137 } |
137 |
138 |
138 if (splitedge == edges[2].end()) |
139 if (splitedge == edges[2].end()) |
139 { |
140 { |
140 setLeaf(); |
141 node->setLeaf(); |
141 return; |
142 return; |
142 } |
143 } |
|
144 |
|
145 node->setSplit(splitedge->pos); |
143 |
146 |
144 #if 0 |
147 #if 0 |
145 // export kd-tree as .obj for visualization |
148 // export kd-tree as .obj for visualization |
146 // this would be hard to reconstruct later |
149 // this would be hard to reconstruct later |
147 static ofstream F("kdtree.obj"); |
150 static ofstream F("kdtree.obj"); |
148 Vector3 v; |
151 Vector3 v; |
149 static int f=0; |
152 static int f=0; |
150 v.cell[axis] = split; |
153 v.cell[axis] = node->getSplit(); |
151 v.cell[(axis+1)%3] = bounds.L.cell[(axis+1)%3]; |
154 v.cell[(axis+1)%3] = bounds.L.cell[(axis+1)%3]; |
152 v.cell[(axis+2)%3] = bounds.L.cell[(axis+2)%3]; |
155 v.cell[(axis+2)%3] = bounds.L.cell[(axis+2)%3]; |
153 F << "v " << v.x << " " << v.y << " " << v.z << endl; |
156 F << "v " << v.x << " " << v.y << " " << v.z << endl; |
154 v.cell[(axis+1)%3] = bounds.L.cell[(axis+1)%3]; |
157 v.cell[(axis+1)%3] = bounds.L.cell[(axis+1)%3]; |
155 v.cell[(axis+2)%3] = bounds.H.cell[(axis+2)%3]; |
158 v.cell[(axis+2)%3] = bounds.H.cell[(axis+2)%3]; |
163 F << "f " << f+1 << " " << f+2 << " " << f+3 << " " << f+4 << endl; |
166 F << "f " << f+1 << " " << f+2 << " " << f+3 << " " << f+4 << endl; |
164 f += 4; |
167 f += 4; |
165 #endif |
168 #endif |
166 |
169 |
167 // split this node |
170 // split this node |
168 delete getShapes(); |
171 delete shapes; |
169 |
172 |
170 BBox lbb = bounds; |
173 BBox lbb = bounds; |
171 BBox rbb = bounds; |
174 BBox rbb = bounds; |
172 lbb.H.cell[axis] = split; |
175 lbb.H.cell[axis] = node->getSplit(); |
173 rbb.L.cell[axis] = split; |
176 rbb.L.cell[axis] = node->getSplit(); |
174 setChildren(new KdNode[2]); |
177 node->setChildren(new KdNode[2]); |
175 setAxis(axis); |
178 node->setAxis(axis); |
176 |
179 |
177 for (edge = edges[axis].begin(); edge != splitedge; edge++) |
180 for (edge = edges[axis].begin(); edge != splitedge; edge++) |
178 if (!edge->end && edge->shape->intersect_bbox(lbb)) |
181 if (!edge->end && edge->shape->intersect_bbox(lbb)) |
179 getLeftChild()->addShape(edge->shape); |
182 node->getLeftChild()->addShape(edge->shape); |
180 for (edge = splitedge; edge < edges[axis].end(); edge++) |
183 for (edge = splitedge; edge < edges[axis].end(); edge++) |
181 if (edge->end && edge->shape->intersect_bbox(rbb)) |
184 if (edge->end && edge->shape->intersect_bbox(rbb)) |
182 getRightChild()->addShape(edge->shape); |
185 node->getRightChild()->addShape(edge->shape); |
183 |
186 |
184 getLeftChild()->subdivide(lbb, maxdepth-1); |
187 recursive_build(node->getLeftChild(), lbb, maxdepth-1); |
185 getRightChild()->subdivide(rbb, maxdepth-1); |
188 recursive_build(node->getRightChild(), rbb, maxdepth-1); |
186 } |
189 } |
187 |
190 |
188 void KdTree::build() |
191 void KdTree::build() |
189 { |
192 { |
190 dbgmsg(1, "* building kd-tree\n"); |
193 dbgmsg(1, "* building kd-tree\n"); |
191 root = new KdNode(); |
194 root = new KdNode(); |
192 ShapeList::iterator shape; |
195 ShapeList::iterator shape; |
193 for (shape = shapes.begin(); shape != shapes.end(); shape++) |
196 for (shape = shapes.begin(); shape != shapes.end(); shape++) |
194 root->addShape(*shape); |
197 root->addShape(*shape); |
195 root->subdivide(bbox, max_depth); |
198 recursive_build(root, bbox, max_depth); |
196 built = true; |
199 built = true; |
197 } |
200 } |
198 |
201 |
199 /* algorithm by Vlastimil Havran, Heuristic Ray Shooting Algorithms, appendix C */ |
202 /* algorithm by Vlastimil Havran, Heuristic Ray Shooting Algorithms, appendix C */ |
200 Shape *KdTree::nearest_intersection(const Shape *origin_shape, const Ray &ray, |
203 Shape *KdTree::nearest_intersection(const Shape *origin_shape, const Ray &ray, |