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

Revision 2629, 11.3 KB checked in by bittner, 17 years ago (diff)

commit after merge with vlastimil

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(); fi != leaf->mFaces.end(); fi++)
151  {
152          // determine the side of this ray with respect to the plane
153          AxisAlignedBox3 box = mMesh->GetFaceBox(*fi);
154
155          if (box.Max(axis) > position )
156                  objectsFront.push_back(*fi);
157
158          if (box.Min(axis) < position )
159                  objectsBack.push_back(*fi);
160  }
161
162  MeshKdLeaf *back = new MeshKdLeaf(objectsBack);
163  MeshKdLeaf *front = new MeshKdLeaf(objectsFront);
164
165  // replace a link from node's parent
166  if (  parent )
167    parent->ReplaceChildLink(leaf, node);
168 
169  // and setup child links
170  node->SetupChildLinks(back, front);
171 
172  delete leaf;
173  return node;
174}
175
176
177
178
179void
180MeshKdTree::SortSubdivisionCandidates(
181                            MeshKdLeaf *node,
182                            const int axis
183                            )
184{
185  mSubdivisionCandidates->clear();
186 
187  int requestedSize = 2*(int)node->mFaces.size();
188  // creates a sorted split candidates array
189  if (mSubdivisionCandidates->capacity() > 500000 &&
190      requestedSize < (int)(mSubdivisionCandidates->capacity()/10) ) {
191    delete mSubdivisionCandidates;
192    mSubdivisionCandidates = new vector<SortableEntry>;
193  }
194 
195  mSubdivisionCandidates->reserve(requestedSize);
196 
197  vector<int>::const_iterator fi;
198  // insert all queries
199  for(fi = node->mFaces.begin();
200      fi < node->mFaces.end();
201      fi++) {
202   
203    AxisAlignedBox3 box = mMesh->GetFaceBox(*fi);
204   
205    mSubdivisionCandidates->push_back(SortableEntry(SortableEntry::FACE_MIN,
206                                             box.Min(axis),
207                                             *fi)
208                               );
209   
210   
211    mSubdivisionCandidates->push_back(SortableEntry(SortableEntry::FACE_MAX,
212                                             box.Max(axis),
213                                             *fi)
214                               );
215  }
216 
217  stable_sort(mSubdivisionCandidates->begin(), mSubdivisionCandidates->end());
218}
219
220
221float
222MeshKdTree::BestCostRatio(
223                          MeshKdLeaf *node,
224                          const AxisAlignedBox3 &box,
225                          const int axis,
226                          float &position,
227                          int &objectsBack,
228                          int &objectsFront
229                          )
230{
231
232  SortSubdivisionCandidates(node, axis);
233 
234  // go through the lists, count the number of objects left and right
235  // and evaluate the following cost funcion:
236  // C = ct_div_ci  + (ol + or)/queries
237
238  int objectsLeft = 0, objectsRight = (int)node->mFaces.size();
239 
240  float minBox = box.Min(axis);
241  float maxBox = box.Max(axis);
242  float boxArea = box.SurfaceArea();
243 
244  float minBand = minBox + mSplitBorder*(maxBox - minBox);
245  float maxBand = minBox + (1.0f - mSplitBorder)*(maxBox - minBox);
246 
247  float minSum = 1e20f;
248  vector<SortableEntry>::const_iterator ci;
249
250  for(ci = mSubdivisionCandidates->begin();
251      ci != mSubdivisionCandidates->end();
252      ci++) {
253    switch ((*ci).type) {
254    case SortableEntry::FACE_MIN:
255      objectsLeft++;
256      break;
257    case SortableEntry::FACE_MAX:
258      objectsRight--;
259      break;
260    }
261   
262    if ((*ci).value > minBand && (*ci).value < maxBand) {
263      AxisAlignedBox3 lbox = box;
264      AxisAlignedBox3 rbox = box;
265      lbox.SetMax(axis, (*ci).value);
266      rbox.SetMin(axis, (*ci).value);
267
268      float sum = objectsLeft*lbox.SurfaceArea() + objectsRight*rbox.SurfaceArea();
269     
270      //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
271      //      cout<<"cost= "<<sum<<endl;
272     
273      if (sum < minSum) {
274        minSum = sum;
275        position = (*ci).value;
276       
277        objectsBack = objectsLeft;
278        objectsFront = objectsRight;
279      }
280    }
281  }
282 
283  float oldCost = (float)node->mFaces.size();
284  float newCost = mCt_div_ci + minSum/boxArea;
285  float ratio = newCost/oldCost;
286
287#if 0
288  cout<<"===================="<<endl;
289  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
290      <<"\t o=("<<objectsBack<<","<<objectsFront<<")"<<endl;
291#endif
292  return ratio;
293}
294
295int
296MeshKdTree::CastRay(
297                                        Ray &ray,
298                                        MeshInstance *instance
299                                        )
300{
301  int hits = 0;
302 
303  stack<RayTraversalData> tStack;
304 
305  float maxt = 1e6;
306  float mint = 0;
307 
308  AxisAlignedBox3 box = GetBox();
309
310  if (!box.GetMinMaxT(ray, &mint, &maxt))
311    return 0;
312
313  if (mint < 0)
314    mint = 0;
315 
316
317  if (ray.GetType() == Ray::LOCAL_RAY &&
318          ray.intersections.size() && ray.intersections[0].mT < mint) {
319        return 0;
320  }
321 
322  maxt += Limits::Threshold;
323 
324  Vector3 entp = ray.Extrap(mint);
325  Vector3 extp = ray.Extrap(maxt);
326 
327  MeshKdNode *node = mRoot;
328  MeshKdNode *farChild;
329  float position;
330  int axis;
331
332  while (1) {
333    if (!node->IsLeaf()) {
334      MeshKdInterior *in = (MeshKdInterior *) node;
335      position = in->mPosition;
336      axis = in->mAxis;
337                       
338      if (entp[axis] <= position) {
339                                if (extp[axis] <= position) {
340                                        node = in->mBack;
341                                        // cases N1,N2,N3,P5,Z2,Z3
342                                        continue;
343                                } else {
344                                        // case N4
345                                        node = in->mBack;
346                                        farChild = in->mFront;
347                                }
348      }
349      else {
350                                if (position <= extp[axis]) {
351                                        node = in->mFront;
352                                        // cases P1,P2,P3,N5,Z1
353                                        continue;
354                                } else {
355                                        node = in->mFront;
356                                        farChild = in->mBack;
357                                        // case P4
358                                }
359      }
360      // $$ modification 3.5.2004 - hints from Kamil Ghais
361      // case N4 or P4
362      float tdist = (position - ray.GetLoc(axis)) / ray.GetDir(axis);
363      tStack.push(RayTraversalData(farChild, extp, maxt));
364      extp = ray.GetLoc() + ray.GetDir()*tdist;
365      maxt = tdist;
366    } else {
367      // compute intersection with all objects in this leaf
368      MeshKdLeaf *leaf = (MeshKdLeaf *) node;
369      //      cout<<"leaf mfaces size="<<leaf->mFaces.size()<<endl<<flush;
370      hits += instance->CastRay(ray, leaf->mFaces);
371     
372      if (ray.GetType() == Ray::LOCAL_RAY && ray.intersections.size())
373                                if (ray.intersections[0].mT <= maxt) {
374                                        break;
375                                }
376     
377      // get the next node from the stack
378      if (tStack.empty())
379                                break;
380     
381      entp = extp;
382      mint = maxt;
383
384      if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
385                                break;
386                       
387      RayTraversalData &s  = tStack.top();
388      node = s.mNode;
389      extp = s.mExitPoint;
390      maxt = s.mMaxT;
391      tStack.pop();
392    }
393  }
394 
395 
396  return hits;
397}
398
399bool
400MeshKdTree::TerminationCriteriaMet(const MeshKdLeaf *leaf, const int depth)
401{
402  //  cerr<<"\n OBJECTS="<<leaf->mObjects.size()<<endl;
403  return
404    (leaf->mFaces.size() <= mTermMinCost) ||
405    (depth >= mTermMaxDepth);
406 
407}
408
409bool
410MeshKdTree::Construct()
411{
412  if (!mSubdivisionCandidates)
413    mSubdivisionCandidates = new vector<SortableEntry>;
414 
415  // first construct a leaf that will get subdivide
416  MeshKdLeaf *leaf = (MeshKdLeaf *) mRoot;
417
418  mRoot = Subdivide(TraversalData(leaf, NULL, GetBox(), 0));
419 
420  cout <<  "created " << endl;
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.