source: trunk/VUT/GtpVisibilityPreprocessor/src/MeshKdTree.cpp @ 177

Revision 177, 10.4 KB checked in by bittner, 19 years ago (diff)
Line 
1#include <stack>
2#include <algorithm>
3#include <queue>
4#include "Environment.h"
5#include "Mesh.h"
6#include "MeshKdTree.h"
7
8float MeshKdTree::mSplitBorder;
9int MeshKdTree::mTermMaxDepth;
10int MeshKdTree::mTermMinCost;
11float MeshKdTree::mMaxCostRatio;
12float MeshKdTree::mCt_div_ci;
13int MeshKdTree::mSplitMethod;
14
15MeshKdTree::MeshKdTree(Mesh *mesh):mMesh(mesh)
16{
17  mRoot = new MeshKdLeaf;
18  mSplitCandidates = NULL;
19}
20
21void
22MeshKdTree::ParseEnvironment()
23{
24  environment->GetIntValue("MeshKdTree.Termination.maxDepth", mTermMaxDepth);
25  environment->GetIntValue("MeshKdTree.Termination.minCost", mTermMinCost);
26  environment->GetFloatValue("MeshKdTree.Termination.maxCostRatio", mMaxCostRatio);
27  environment->GetFloatValue("MeshKdTree.Termination.ct_div_ci", mCt_div_ci);
28  environment->GetFloatValue("MeshKdTree.splitBorder", mSplitBorder);
29 
30  char splitType[64];
31  environment->GetStringValue("MeshKdTree.splitMethod", splitType);
32 
33  mSplitMethod = SPLIT_SPATIAL_MEDIAN;
34  if (strcmp(splitType, "spatialMedian") == 0)
35    mSplitMethod = SPLIT_SPATIAL_MEDIAN;
36  else
37    if (strcmp(splitType, "objectMedian") == 0)
38      mSplitMethod = SPLIT_OBJECT_MEDIAN;
39    else
40      if (strcmp(splitType, "SAH") == 0)
41        mSplitMethod = SPLIT_SAH;
42      else {
43        cerr<<"Wrong kd split type "<<splitType<<endl;
44        exit(1);
45      }
46}
47
48
49
50int
51MeshKdTree::SelectPlane(MeshKdLeaf *leaf,
52                        const AxisAlignedBox3 &box,
53                        float &position
54                        )
55{
56  int axis = -1;
57 
58  switch (mSplitMethod)
59    {
60    case SPLIT_SPATIAL_MEDIAN: {
61      axis = box.Size().DrivingAxis();
62      position = (box.Min()[axis] + box.Max()[axis])*0.5f;
63      break;
64    }
65    case SPLIT_SAH: {
66      int objectsBack, objectsFront;
67      float costRatio;
68      bool mOnlyDrivingAxis = false;
69      if (mOnlyDrivingAxis) {
70        axis = box.Size().DrivingAxis();
71        costRatio = BestCostRatio(leaf,
72                                  box,
73                                  axis,
74                                  position,
75                                  objectsBack,
76                                  objectsFront);
77      } else {
78        costRatio = MAX_FLOAT;
79        for (int i=0; i < 3; i++) {
80          float p;
81          float r = BestCostRatio(leaf,
82                                  box,
83                                  i,
84                                  p,
85                                  objectsBack,
86                                  objectsFront);
87          if (r < costRatio) {
88            costRatio = r;
89            axis = i;
90            position = p;
91          }
92        }
93      }
94     
95      if (costRatio > mMaxCostRatio) {
96        //      cout<<"Too big cost ratio "<<costRatio<<endl;
97        axis = -1;
98      }
99      break;
100    }
101   
102    }
103  return axis;
104}
105
106MeshKdNode *
107MeshKdTree::SubdivideNode(
108                          MeshKdLeaf *leaf,
109                          MeshKdInterior *parent,
110                          const AxisAlignedBox3 &box,
111                          const int depth,
112                          AxisAlignedBox3 &backBBox,
113                          AxisAlignedBox3 &frontBBox
114                          )
115{
116 
117  if (TerminationCriteriaMet(leaf, depth))
118    return leaf;
119 
120  float position;
121  // select subdivision axis
122  int axis = SelectPlane( leaf, box, position );
123 
124  if (axis == -1) {
125    return leaf;
126  }
127 
128  // add the new nodes to the tree
129  MeshKdInterior *node = new MeshKdInterior;
130 
131  node->mAxis = axis;
132  node->mPosition = position;
133 
134  backBBox = box;
135  frontBBox = box;
136 
137  // first count ray sides
138 
139  backBBox.SetMax(axis, position);
140  frontBBox.SetMin(axis, position);
141
142  vector<int>::const_iterator fi;
143  vector<int> objectsFront, objectsBack;
144 
145  for ( fi = leaf->mFaces.begin();
146        fi != leaf->mFaces.end();
147        fi++) {
148    // determine the side of this ray with respect to the plane
149    AxisAlignedBox3 box = mMesh->GetFaceBox(*fi);
150
151    if (box.Max(axis) > position )
152      objectsFront.push_back(*fi);
153   
154    if (box.Min(axis) < position )
155      objectsBack.push_back(*fi);
156  }
157 
158  MeshKdLeaf *back = new MeshKdLeaf(objectsBack);
159  MeshKdLeaf *front = new MeshKdLeaf(objectsFront);
160
161  // replace a link from node's parent
162  if (  parent )
163    parent->ReplaceChildLink(leaf, node);
164 
165  // and setup child links
166  node->SetupChildLinks(back, front);
167 
168  delete leaf;
169  return node;
170}
171
172
173
174
175void
176MeshKdTree::SortSplitCandidates(
177                            MeshKdLeaf *node,
178                            const int axis
179                            )
180{
181  mSplitCandidates->clear();
182 
183  int requestedSize = 2*node->mFaces.size();
184  // creates a sorted split candidates array
185  if (mSplitCandidates->capacity() > 500000 &&
186      requestedSize < (int)(mSplitCandidates->capacity()/10) ) {
187    delete mSplitCandidates;
188    mSplitCandidates = new vector<SortableEntry>;
189  }
190 
191  mSplitCandidates->reserve(requestedSize);
192 
193  vector<int>::const_iterator fi;
194  // insert all queries
195  for(fi = node->mFaces.begin();
196      fi < node->mFaces.end();
197      fi++) {
198   
199    AxisAlignedBox3 box = mMesh->GetFaceBox(*fi);
200   
201    mSplitCandidates->push_back(SortableEntry(SortableEntry::FACE_MIN,
202                                             box.Min(axis),
203                                             *fi)
204                               );
205   
206   
207    mSplitCandidates->push_back(SortableEntry(SortableEntry::FACE_MAX,
208                                             box.Max(axis),
209                                             *fi)
210                               );
211  }
212 
213  stable_sort(mSplitCandidates->begin(), mSplitCandidates->end());
214}
215
216
217float
218MeshKdTree::BestCostRatio(
219                          MeshKdLeaf *node,
220                          const AxisAlignedBox3 &box,
221                          const int axis,
222                          float &position,
223                          int &objectsBack,
224                          int &objectsFront
225                          )
226{
227
228  SortSplitCandidates(node, axis);
229 
230  // go through the lists, count the number of objects left and right
231  // and evaluate the following cost funcion:
232  // C = ct_div_ci  + (ol + or)/queries
233
234  int objectsLeft = 0, objectsRight = node->mFaces.size();
235 
236  float minBox = box.Min(axis);
237  float maxBox = box.Max(axis);
238  float boxArea = box.SurfaceArea();
239 
240  float minBand = minBox + mSplitBorder*(maxBox - minBox);
241  float maxBand = minBox + (1.0f - mSplitBorder)*(maxBox - minBox);
242 
243  float minSum = 1e20;
244  vector<SortableEntry>::const_iterator ci;
245
246  for(ci = mSplitCandidates->begin();
247      ci != mSplitCandidates->end();
248      ci++) {
249    switch ((*ci).type) {
250    case SortableEntry::FACE_MIN:
251      objectsLeft++;
252      break;
253    case SortableEntry::FACE_MAX:
254      objectsRight--;
255      break;
256    }
257   
258    if ((*ci).value > minBand && (*ci).value < maxBand) {
259      AxisAlignedBox3 lbox = box;
260      AxisAlignedBox3 rbox = box;
261      lbox.SetMax(axis, (*ci).value);
262      rbox.SetMin(axis, (*ci).value);
263
264      float sum = objectsLeft*lbox.SurfaceArea() + objectsRight*rbox.SurfaceArea();
265     
266      //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
267      //      cout<<"cost= "<<sum<<endl;
268     
269      if (sum < minSum) {
270        minSum = sum;
271        position = (*ci).value;
272       
273        objectsBack = objectsLeft;
274        objectsFront = objectsRight;
275      }
276    }
277  }
278 
279  float oldCost = node->mFaces.size();
280  float newCost = mCt_div_ci + minSum/boxArea;
281  float ratio = newCost/oldCost;
282
283#if 0
284  cout<<"===================="<<endl;
285  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
286      <<"\t o=("<<objectsBack<<","<<objectsFront<<")"<<endl;
287#endif
288  return ratio;
289}
290
291int
292MeshKdTree::CastRay(
293                    Ray &ray,
294                    MeshInstance *instance
295                    )
296{
297  int hits = 0;
298 
299  stack<RayTraversalData> tStack;
300 
301  float maxt = 1e6;
302  float mint = 0;
303 
304  AxisAlignedBox3 box = GetBox();
305  if (!box.GetMinMaxT(ray, &mint, &maxt))
306    return 0;
307
308  if (mint < 0)
309    mint = 0;
310 
311  maxt += Limits::Threshold;
312 
313  Vector3 entp = ray.Extrap(mint);
314  Vector3 extp = ray.Extrap(maxt);
315 
316  MeshKdNode *node = mRoot;
317  MeshKdNode *farChild;
318  float position;
319  int axis;
320 
321  while (1) {
322    if (!node->IsLeaf()) {
323      MeshKdInterior *in = (MeshKdInterior *) node;
324      position = in->mPosition;
325      axis = in->mAxis;
326
327      if (entp[axis] <= position) {
328        if (extp[axis] <= position) {
329          node = in->mBack;
330          // cases N1,N2,N3,P5,Z2,Z3
331          continue;
332        } else {
333          // case N4
334          node = in->mBack;
335          farChild = in->mFront;
336        }
337      }
338      else {
339        if (position <= extp[axis]) {
340          node = in->mFront;
341          // cases P1,P2,P3,N5,Z1
342          continue;
343        } else {
344          node = in->mFront;
345          farChild = in->mBack;
346          // case P4
347        }
348      }
349      // $$ modification 3.5.2004 - hints from Kamil Ghais
350      // case N4 or P4
351      float tdist = (position - ray.GetLoc(axis)) / ray.GetDir(axis);
352      tStack.push(RayTraversalData(farChild, extp, maxt));
353      extp = ray.GetLoc() + ray.GetDir()*tdist;
354      maxt = tdist;
355    } else {
356      // compute intersection with all objects in this leaf
357      MeshKdLeaf *leaf = (MeshKdLeaf *) node;
358      //      cout<<"leaf mfaces size="<<leaf->mFaces.size()<<endl<<flush;
359      hits += instance->CastRay(ray, leaf->mFaces);
360     
361      if (hits && ray.GetType() == Ray::LOCAL_RAY)
362        if (ray.intersections[0].mT <= maxt)
363          break;
364     
365      // get the next node from the stack
366      if (tStack.empty())
367        break;
368     
369      entp = extp;
370      mint = maxt;
371      RayTraversalData &s  = tStack.top();
372      node = s.mNode;
373      extp = s.mExitPoint;
374      maxt = s.mMaxT;
375      tStack.pop();
376    }
377  }
378 
379 
380  return hits;
381}
382
383bool
384MeshKdTree::TerminationCriteriaMet(const MeshKdLeaf *leaf, const int depth)
385{
386  //  cerr<<"\n OBJECTS="<<leaf->mObjects.size()<<endl;
387  return
388    (leaf->mFaces.size() <= mTermMinCost) ||
389    (depth >= mTermMaxDepth);
390 
391}
392
393bool
394MeshKdTree::Construct()
395{
396  if (!mSplitCandidates)
397    mSplitCandidates = new vector<SortableEntry>;
398 
399  // first construct a leaf that will get subdivide
400  MeshKdLeaf *leaf = (MeshKdLeaf *) mRoot;
401
402  mRoot = Subdivide(TraversalData(leaf, NULL, GetBox(), 0));
403 
404  // remove the allocated array
405  delete mSplitCandidates;
406  mSplitCandidates = NULL;
407 
408  return true;
409}
410
411
412MeshKdNode *
413MeshKdTree::Subdivide(const TraversalData &tdata)
414{
415  MeshKdNode *result = NULL;
416 
417  priority_queue<TraversalData> tStack;
418  //  stack<STraversalData> tStack;
419 
420  tStack.push(tdata);
421  AxisAlignedBox3 backBox, frontBox;
422
423 
424  while (!tStack.empty()) {
425
426    TraversalData data = tStack.top();
427    tStack.pop();
428   
429    MeshKdNode *node = SubdivideNode((MeshKdLeaf *) data.mNode,
430                                     data.mParent,
431                                     data.mBox,
432                                     data.mDepth,
433                                     backBox,
434                                     frontBox
435                                     );
436    if (result == NULL)
437      result = node;
438   
439    if (!node->IsLeaf()) {
440
441      MeshKdInterior *interior = (MeshKdInterior *) node;
442      // push the children on the stack
443      tStack.push(TraversalData(interior->mBack, interior, backBox, data.mDepth+1));
444      tStack.push(TraversalData(interior->mFront, interior, frontBox, data.mDepth+1));
445     
446    }
447  }
448 
449  return result;
450
451}
Note: See TracBrowser for help on using the repository browser.