source: GTP/trunk/Lib/Vis/Preprocessing/src/MeshKdTree.cpp @ 2176

Revision 2176, 11.3 KB checked in by mattausch, 18 years ago (diff)

removed using namespace std from .h

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