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

Revision 863, 10.7 KB checked in by mattausch, 19 years ago (diff)

working on preprocessor integration
added iv stuff

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