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

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

updated view cells, view cell manager. changed rendersimulator

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*(int)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 = (int)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 = 1e20f;
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 = (float)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
312        if (ray.GetType() == Ray::LOCAL_RAY &&
313                        ray.intersections.size() && ray.intersections[0].mT < mint) {
314                return 0;
315        }
316
317  maxt += Limits::Threshold;
318 
319  Vector3 entp = ray.Extrap(mint);
320  Vector3 extp = ray.Extrap(maxt);
321 
322  MeshKdNode *node = mRoot;
323  MeshKdNode *farChild;
324  float position;
325  int axis;
326 
327  while (1) {
328    if (!node->IsLeaf()) {
329      MeshKdInterior *in = (MeshKdInterior *) node;
330      position = in->mPosition;
331      axis = in->mAxis;
332                       
333      if (entp[axis] <= position) {
334                                if (extp[axis] <= position) {
335                                        node = in->mBack;
336                                        // cases N1,N2,N3,P5,Z2,Z3
337                                        continue;
338                                } else {
339                                        // case N4
340                                        node = in->mBack;
341                                        farChild = in->mFront;
342                                }
343      }
344      else {
345                                if (position <= extp[axis]) {
346                                        node = in->mFront;
347                                        // cases P1,P2,P3,N5,Z1
348                                        continue;
349                                } else {
350                                        node = in->mFront;
351                                        farChild = in->mBack;
352                                        // case P4
353                                }
354      }
355      // $$ modification 3.5.2004 - hints from Kamil Ghais
356      // case N4 or P4
357      float tdist = (position - ray.GetLoc(axis)) / ray.GetDir(axis);
358      tStack.push(RayTraversalData(farChild, extp, maxt));
359      extp = ray.GetLoc() + ray.GetDir()*tdist;
360      maxt = tdist;
361    } else {
362      // compute intersection with all objects in this leaf
363      MeshKdLeaf *leaf = (MeshKdLeaf *) node;
364      //      cout<<"leaf mfaces size="<<leaf->mFaces.size()<<endl<<flush;
365      hits += instance->CastRay(ray, leaf->mFaces);
366     
367      if (ray.GetType() == Ray::LOCAL_RAY && ray.intersections.size())
368                                if (ray.intersections[0].mT <= maxt) {
369                                        break;
370                                }
371     
372      // get the next node from the stack
373      if (tStack.empty())
374                                break;
375     
376      entp = extp;
377      mint = maxt;
378
379      if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
380                                break;
381                       
382      RayTraversalData &s  = tStack.top();
383      node = s.mNode;
384      extp = s.mExitPoint;
385      maxt = s.mMaxT;
386      tStack.pop();
387    }
388  }
389 
390 
391  return hits;
392}
393
394bool
395MeshKdTree::TerminationCriteriaMet(const MeshKdLeaf *leaf, const int depth)
396{
397  //  cerr<<"\n OBJECTS="<<leaf->mObjects.size()<<endl;
398  return
399    (leaf->mFaces.size() <= mTermMinCost) ||
400    (depth >= mTermMaxDepth);
401 
402}
403
404bool
405MeshKdTree::Construct()
406{
407  if (!mSplitCandidates)
408    mSplitCandidates = new vector<SortableEntry>;
409 
410  // first construct a leaf that will get subdivide
411  MeshKdLeaf *leaf = (MeshKdLeaf *) mRoot;
412
413  mRoot = Subdivide(TraversalData(leaf, NULL, GetBox(), 0));
414 
415  // remove the allocated array
416  delete mSplitCandidates;
417  mSplitCandidates = NULL;
418 
419  return true;
420}
421
422
423MeshKdNode *
424MeshKdTree::Subdivide(const TraversalData &tdata)
425{
426  MeshKdNode *result = NULL;
427 
428  priority_queue<TraversalData> tStack;
429  //  stack<STraversalData> tStack;
430 
431  tStack.push(tdata);
432  AxisAlignedBox3 backBox, frontBox;
433
434 
435  while (!tStack.empty()) {
436
437    TraversalData data = tStack.top();
438    tStack.pop();
439   
440    MeshKdNode *node = SubdivideNode((MeshKdLeaf *) data.mNode,
441                                     data.mParent,
442                                     data.mBox,
443                                     data.mDepth,
444                                     backBox,
445                                     frontBox
446                                     );
447    if (result == NULL)
448      result = node;
449   
450    if (!node->IsLeaf()) {
451
452      MeshKdInterior *interior = (MeshKdInterior *) node;
453      // push the children on the stack
454      tStack.push(TraversalData(interior->mBack, interior, backBox, data.mDepth+1));
455      tStack.push(TraversalData(interior->mFront, interior, frontBox, data.mDepth+1));
456     
457    }
458  }
459 
460  return result;
461
462}
Note: See TracBrowser for help on using the repository browser.