source: trunk/VUT/GtpVisibilityPreprocessor/src/VssTree.cpp @ 446

Revision 446, 45.6 KB checked in by bittner, 19 years ago (diff)

non-functional merged version
Line 
1
2// ================================================================
3// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
4// ****************************************************************
5/**
6   The KD tree based LSDS
7*/
8// Initial coding by
9/**
10   @author Jiri Bittner
11*/
12
13// Standard headers
14#include <stack>
15#include <queue>
16#include <algorithm>
17#include <fstream>
18#include <string>
19
20#include "VssTree.h"
21
22#include "Environment.h"
23#include "VssRay.h"
24#include "Intersectable.h"
25#include "Ray.h"
26
27#define DEBUG_SPLIT_COST 0
28
29// Static variables
30int
31VssTreeLeaf::mailID = 0;
32
33inline void
34AddObject2Pvs(Intersectable *object,
35                          const int side,
36                          int &pvsBack,
37                          int &pvsFront)
38{
39 
40  if (!object)
41        return;
42 
43  if (side <= 0) {
44        if (!object->Mailed() && !object->Mailed(2)) {
45          pvsBack++;
46          if (object->Mailed(1))
47                object->Mail(2);
48          else
49                object->Mail();
50        }
51  }
52 
53  if (side >= 0) {
54        if (!object->Mailed(1) && !object->Mailed(2)) {
55          pvsFront++;
56          if (object->Mailed())
57                object->Mail(2);
58          else
59                object->Mail(1);
60        }
61  }
62}
63
64// Constructor
65VssTree::VssTree()
66{
67  environment->GetIntValue("VssTree.maxDepth", termMaxDepth);
68  environment->GetIntValue("VssTree.minPvs", termMinPvs);
69  environment->GetIntValue("VssTree.minRays", termMinRays);
70  environment->GetFloatValue("VssTree.maxRayContribution", termMaxRayContribution);
71  environment->GetFloatValue("VssTree.maxCostRatio", termMaxCostRatio);
72
73  environment->GetFloatValue("VssTree.minSize", termMinSize);
74  termMinSize = sqr(termMinSize);
75       
76  environment->GetFloatValue("VssTree.refDirBoxMaxSize", refDirBoxMaxSize);
77  refDirBoxMaxSize = sqr(refDirBoxMaxSize);
78 
79  environment->GetFloatValue("VssTree.epsilon", epsilon);
80  environment->GetFloatValue("VssTree.ct_div_ci", ct_div_ci);
81       
82  environment->GetFloatValue("VssTree.maxTotalMemory", maxTotalMemory);
83  environment->GetFloatValue("VssTree.maxStaticMemory", maxStaticMemory);
84 
85 
86
87
88  float refDirAngle;
89  environment->GetFloatValue("VssTree.refDirAngle", refDirAngle);
90 
91  environment->GetIntValue("VssTree.accessTimeThreshold", accessTimeThreshold);
92  //= 1000;
93  environment->GetIntValue("VssTree.minCollapseDepth", minCollapseDepth);
94  //  int minCollapseDepth = 4;
95
96  //  pRefDirThresh = cos(0.5*M_PI - M_PI*refDirAngle/180.0);
97  //  cosRefDir = cos(M_PI*refDirAngle/180.0);
98  //  sinRefDir = sin(M_PI*refDirAngle/180.0);
99 
100 
101  // split type
102  char sname[128];
103  environment->GetStringValue("VssTree.splitType", sname);
104  string name(sname);
105       
106  if (name.compare("regular") == 0)
107    splitType = ESplitRegular;
108  else
109    if (name.compare("heuristic") == 0)
110      splitType = ESplitHeuristic;
111        else
112          if (name.compare("hybrid") == 0)
113                splitType = ESplitHybrid;
114          else {
115                cerr<<"Invalid VssTree split type "<<name<<endl;
116                exit(1);
117          }
118       
119  environment->GetBoolValue("VssTree.randomize", randomize);
120  environment->GetBoolValue("VssTree.splitUseOnlyDrivingAxis", mSplitUseOnlyDrivingAxis);
121  environment->GetBoolValue("VssTree.useRss", mUseRss);
122
123  environment->GetBoolValue("VssTree.interleaveDirSplits", mInterleaveDirSplits);
124  environment->GetIntValue("VssTree.dirSplitDepth", mDirSplitDepth);
125
126  root = NULL;
127 
128  splitCandidates = new vector<SortableEntry>;
129}
130
131
132VssTree::~VssTree()
133{
134  if (root)
135    delete root;
136}
137
138
139
140
141void
142VssStatistics::Print(ostream &app) const
143{
144  app << "===== VssTree statistics ===============\n";
145
146  app << "#N_RAYS ( Number of rays )\n"
147      << rays <<endl;
148
149  app << "#N_INITPVS ( Initial PVS size )\n"
150      << initialPvsSize <<endl;
151 
152  app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
153
154  app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
155
156  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz \n";
157  for (int i=0; i<7; i++)
158    app << splits[i] <<" ";
159  app <<endl;
160
161  app << "#N_RAYREFS ( Number of rayRefs )\n" <<
162    rayRefs << "\n";
163
164  app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
165    rayRefs/(double)rays << "\n";
166
167  app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
168    rayRefs/(double)Leaves() << "\n";
169
170  app << "#N_MAXRAYREFS  ( Max number of rayRefs / leaf )\n" <<
171    maxRayRefs << "\n";
172
173
174  //  app << setprecision(4);
175
176  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
177    maxDepthNodes*100/(double)Leaves()<<endl;
178
179  app << "#N_PMINPVSLEAVES  ( Percentage of leaves with minCost )\n"<<
180    minPvsNodes*100/(double)Leaves()<<endl;
181
182  app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minCost )\n"<<
183    minRaysNodes*100/(double)Leaves()<<endl;
184       
185  app << "#N_PMINSIZELEAVES  ( Percentage of leaves with minSize )\n"<<
186    minSizeNodes*100/(double)Leaves()<<endl;
187
188  app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
189    maxRayContribNodes*100/(double)Leaves()<<endl;
190
191  app << "#N_PMAXCOSTRATIOLEAVES  ( Percentage of leaves with max cost ratio )\n"<<
192    maxCostRatioNodes*100/(double)Leaves()<<endl;
193
194  app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
195    addedRayRefs<<endl;
196
197  app << "#N_REMOVED_RAYREFS  (Number of dynamically removed ray references )\n"<<
198    removedRayRefs<<endl;
199
200  //  app << setprecision(4);
201
202  app << "#N_CTIME  ( Construction time [s] )\n"
203      << Time() << " \n";
204
205  app << "===== END OF VssTree statistics ==========\n";
206
207}
208
209
210void
211VssTreeLeaf::UpdatePvsSize()
212{
213  if (!mValidPvs) {
214        Intersectable::NewMail();
215        int pvsSize = 0;
216        for(VssTreeNode::RayInfoContainer::iterator ri = rays.begin();
217                ri != rays.end();
218                ri++)
219          if ((*ri).mRay->IsActive()) {
220                Intersectable *object;
221#if BIDIRECTIONAL_RAY
222                object = (*ri).mRay->mOriginObject;
223                if (object && !object->Mailed()) {
224                  pvsSize++;
225                  object->Mail();
226                }
227#endif
228                object = (*ri).mRay->mTerminationObject;
229                if (object && !object->Mailed()) {
230                  pvsSize++;
231                  object->Mail();
232                }
233          }
234        mPvsSize = pvsSize;
235        mValidPvs = true;
236
237  }
238}
239
240bool
241VssTree::ClipRay(
242                                 VssTreeNode::RayInfo &rayInfo,
243                                 const AxisAlignedBox3 &box
244                                 )
245{
246  float tmin, tmax;
247  static Ray ray;
248  ray.Init(rayInfo.mRay->GetOrigin(), rayInfo.mRay->GetDir(), Ray::LINE_SEGMENT);
249  box.ComputeMinMaxT(ray, &tmin, &tmax);
250  if (tmin >= tmax)
251        return false;
252       
253  if (tmin > 0.0f)
254        rayInfo.SetMinT(tmin);
255
256  if (tmax < 1.0f)
257        rayInfo.SetMaxT(tmax);
258       
259  //    vssRay->SetupEndPoints(
260  //                                                                                             origin,
261  //                                                                                             termination
262  //                                                                                             );
263  return true;
264}
265
266
267void
268VssTree::Construct(
269                                   VssRayContainer &rays,
270                                   AxisAlignedBox3 *forcedBoundingBox
271                                   )
272{
273  stat.Start();
274 
275  maxMemory = maxStaticMemory;
276
277  if (root)
278    delete root;
279       
280  root = new VssTreeLeaf(NULL, rays.size());
281  // first construct a leaf that will get subdivide
282  VssTreeLeaf *leaf = (VssTreeLeaf *) root;
283
284  stat.nodes = 1;
285 
286  bbox.Initialize();
287  dirBBox.Initialize();
288
289  if (mUseRss)
290        forcedBoundingBox = NULL;
291       
292  for(VssRayContainer::const_iterator ri = rays.begin();
293      ri != rays.end();
294      ri++) {
295               
296        VssTreeNode::RayInfo info(*ri);
297        if (forcedBoundingBox)
298          if (!ClipRay(info, *forcedBoundingBox))
299                continue;
300        leaf->AddRay(info);
301               
302    bbox.Include((*ri)->GetOrigin());
303    bbox.Include((*ri)->GetTermination());
304   
305               
306        dirBBox.Include(Vector3(
307                                                        (*ri)->GetDirParametrization(0),
308                                                        (*ri)->GetDirParametrization(1),
309                                                        0
310                                                        )
311                                        );
312  }
313       
314       
315  if ( forcedBoundingBox )
316        bbox = *forcedBoundingBox;
317       
318  cout<<"Bbox = "<<bbox<<endl;
319  cout<<"Dirr Bbox = "<<dirBBox<<endl;
320
321  stat.rays = leaf->rays.size();
322  leaf->UpdatePvsSize();
323  leaf->ComputeEntropyImportance();
324
325  stat.initialPvsSize = leaf->GetPvsSize();
326  // Subdivide();
327  root = Subdivide(TraversalData(leaf, bbox, 0));
328
329  if (splitCandidates) {
330    // force realease of this vector
331    delete splitCandidates;
332    splitCandidates = new vector<SortableEntry>;
333  }
334 
335  stat.Stop();
336
337  stat.Print(cout);
338  cout<<"#Total memory="<<GetMemUsage()<<endl;
339
340}
341
342int
343VssTree::UpdateSubdivision()
344{
345  priority_queue<TraversalData> tStack;
346  //  stack<TraversalData> tStack;
347 
348  tStack.push(TraversalData(root, bbox, 0));
349       
350  AxisAlignedBox3 backBox;
351  AxisAlignedBox3 frontBox;
352
353  maxMemory = maxTotalMemory;
354  int subdivided = 0;
355  int lastMem = 0;
356  while (!tStack.empty()) {
357               
358        float mem = GetMemUsage();
359               
360        if ( lastMem/10 != ((int)mem)/10) {
361          cout<<mem<<" MB"<<endl;
362        }
363        lastMem = (int)mem;
364               
365        if (  mem > maxMemory ) {
366      // count statistics on unprocessed leafs
367      while (!tStack.empty()) {
368                //                              EvaluateLeafStats(tStack.top());
369                tStack.pop();
370      }
371      break;
372    }
373   
374    TraversalData data = tStack.top();
375    tStack.pop();
376
377        if (data.node->IsLeaf()) {
378          VssTreeNode *node = SubdivideNode((VssTreeLeaf *) data.node,
379                                                                                data.bbox,
380                                                                                backBox,
381                                                                                frontBox
382                                                                                );
383          if (!node->IsLeaf()) {
384                subdivided++;
385                VssTreeInterior *interior = (VssTreeInterior *) node;
386                // push the children on the stack
387                tStack.push(TraversalData(interior->back, backBox, data.depth+1));
388                tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
389          } else {
390                //      EvaluateLeafStats(data);
391          }
392        } else {
393          VssTreeInterior *interior = (VssTreeInterior *) data.node;
394          tStack.push(TraversalData(interior->back, GetBBox(interior->back), data.depth+1));
395          tStack.push(TraversalData(interior->front, GetBBox(interior->front), data.depth+1));
396        }
397  }
398  return subdivided;
399}
400
401
402VssTreeNode *
403VssTree::Subdivide(const TraversalData &tdata)
404{
405  VssTreeNode *result = NULL;
406
407  priority_queue<TraversalData> tStack;
408  //  stack<TraversalData> tStack;
409 
410  tStack.push(tdata);
411
412  AxisAlignedBox3 backBox;
413  AxisAlignedBox3 frontBox;
414
415 
416  int lastMem = 0;
417  while (!tStack.empty()) {
418
419        float mem = GetMemUsage();
420               
421        if ( lastMem/10 != ((int)mem)/10) {
422          cout<<mem<<" MB"<<endl;
423        }
424        lastMem = (int)mem;
425               
426        if (  mem > maxMemory ) {
427      // count statistics on unprocessed leafs
428      while (!tStack.empty()) {
429                EvaluateLeafStats(tStack.top());
430                tStack.pop();
431      }
432      break;
433    }
434   
435    TraversalData data = tStack.top();
436    tStack.pop();
437   
438    VssTreeNode *node = SubdivideNode((VssTreeLeaf *) data.node,
439                                                                          data.bbox,
440                                                                          backBox,
441                                                                          frontBox
442                                                                          );
443    if (result == NULL)
444      result = node;
445   
446    if (!node->IsLeaf()) {
447                       
448      VssTreeInterior *interior = (VssTreeInterior *) node;
449      // push the children on the stack
450      tStack.push(TraversalData(interior->back, backBox, data.depth+1));
451      tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
452                       
453    } else {
454      EvaluateLeafStats(data);
455    }
456  }
457
458  return result;
459}
460
461
462// returns selected plane for subdivision
463int
464VssTree::SelectPlane(
465                                         VssTreeLeaf *leaf,
466                                         const AxisAlignedBox3 &box,
467                                         float &position,
468                                         int &raysBack,
469                                         int &raysFront,
470                                         int &pvsBack,
471                                         int &pvsFront
472                                         )
473{
474
475  int minDirDepth = 6;
476  int axis;
477  float costRatio;
478       
479  costRatio = BestCostRatio(leaf,
480                                                        axis,
481                                                        position,
482                                                        raysBack,
483                                                        raysFront,
484                                                        pvsBack,
485                                                        pvsFront
486                                                        );
487#if DEBUG_SPLIT_COST
488  cout<<axis<<" r="<<costRatio<<endl;
489#endif
490  if (costRatio > termMaxCostRatio) {
491        //              cout<<"Too big cost ratio "<<costRatio<<endl;
492        stat.maxCostRatioNodes++;
493        return -1;
494  }
495       
496#if 0   
497  cout<<
498        "pvs="<<leaf->mPvsSize<<
499        " rays="<<leaf->rays.size()<<
500        " rc="<<leaf->GetAvgRayContribution()<<
501        " axis="<<axis<<endl;
502#endif
503       
504  return axis;
505}
506
507
508float
509VssTree::GetCostRatio(
510                                          VssTreeLeaf *leaf,
511                                          const int axis,
512                                          const float position,
513                                          const int raysBack,
514                                          const int raysFront,
515                                          const int pvsBack,
516                                          const int pvsFront
517                                          )
518{
519  bool costImportance = true;
520
521  float ratio;
522  AxisAlignedBox3 box;
523  float minBox, maxBox;
524
525  if (axis < 3) {
526        box = GetBBox(leaf);
527        minBox = box.Min(axis);
528        maxBox = box.Max(axis);
529  }     else {
530        box = GetDirBBox(leaf);
531        minBox = box.Min(axis-3);
532        maxBox = box.Max(axis-3);
533  }
534       
535  float sizeBox = maxBox - minBox;
536
537  int pvsSize = leaf->GetPvsSize();
538
539  if (!costImportance) {
540        //              float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
541        float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
542        float newCost = ct_div_ci + sum/sizeBox;
543        float oldCost = pvsSize;
544        ratio = newCost/oldCost;
545  } else {
546        // importance based cost
547#if 0
548        float newContrib =
549          ((position - minBox)*sqr(pvsBack/(raysBack + Limits::Small)) +
550           (maxBox - position)*sqr(pvsFront/(raysFront + Limits::Small)))/sizeBox;
551               
552        //                      float newContrib =
553        //                              sqr(pvsBack/(raysBack + Limits::Small)) +
554        //                              sqr(pvsFront/(raysFront + Limits::Small));
555        float oldContrib = sqr(leaf->GetAvgRayContribution());
556        ratio = oldContrib/newContrib;
557#else
558#if 1
559        float newCost = raysBack*pvsBack  + raysFront*pvsFront;
560        float oldCost = leaf->rays.size()*pvsSize;
561        ratio = newCost/oldCost;
562#else
563        float newCost = (pvsBack  + pvsFront)*0.5f;
564        float oldCost = pvsSize;
565        ratio = newCost/oldCost;
566#endif
567#endif
568  }
569       
570  return ratio;
571}
572                                                       
573
574float
575VssTree::EvalCostRatio(
576                                           VssTreeLeaf *leaf,
577                                           const int axis,
578                                           const float position,
579                                           int &raysBack,
580                                           int &raysFront,
581                                           int &pvsBack,
582                                           int &pvsFront
583                                           )
584{
585  raysBack = 0;
586  raysFront = 0;
587  pvsFront = 0;
588  pvsBack = 0;
589
590       
591  Intersectable::NewMail(3);
592       
593  if (axis <= VssTreeNode::SPLIT_Z) {
594    // this is the main ray classification loop!
595    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
596                ri != leaf->rays.end();
597                ri++)
598      if ((*ri).mRay->IsActive()) {
599                               
600                // determine the side of this ray with respect to the plane
601                int side = (*ri).ComputeRayIntersection(axis, position, (*ri).mRay->mT);
602                //                              (*ri).mRay->mSide = side;
603                               
604                if (side <= 0)
605                  raysBack++;
606                               
607                if (side >= 0)
608                  raysFront++;
609                               
610                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
611      }
612               
613  } else {
614               
615        // directional split
616    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
617                ri != leaf->rays.end();
618                ri++)
619      if ((*ri).mRay->IsActive()) {
620                               
621                // determine the side of this ray with respect to the plane
622                int side;
623                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
624                  side = 1;
625                else
626                  side = -1;
627
628                if (side <= 0)
629                  raysBack++;
630                               
631                if (side >= 0)
632                  raysFront++;
633
634                //                              (*ri).mRay->mSide = side;
635                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
636                               
637      }
638  }
639
640  float ratio = GetCostRatio(
641                                                         leaf,
642                                                         axis,
643                                                         position,
644                                                         raysBack,
645                                                         raysFront,
646                                                         pvsBack,
647                                                         pvsFront);
648
649  //    cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
650  //  float oldCost = leaf->rays.size();
651
652  //    cout<<"ratio="<<ratio<<endl;
653
654  return ratio;
655}
656
657float
658VssTree::BestCostRatio(
659                                           VssTreeLeaf *leaf,
660                                           int &axis,
661                                           float &position,
662                                           int &raysBack,
663                                           int &raysFront,
664                                           int &pvsBack,
665                                           int &pvsFront
666                                           )
667{
668  int nRaysBack[6], nRaysFront[6];
669  int nPvsBack[6], nPvsFront[6];
670  float nPosition[6];
671  float nCostRatio[6];
672  int bestAxis = -1;
673       
674  AxisAlignedBox3 sBox = GetBBox(leaf);
675  AxisAlignedBox3 dBox = GetDirBBox(leaf);
676  // int sAxis = box.Size().DrivingAxis();
677  int sAxis = sBox.Size().DrivingAxis();
678  int dAxis = dBox.Size().DrivingAxis() + 3;
679
680
681  float dirSplitBoxSize = 0.01f;
682  bool allowDirSplit = Magnitude(sBox.Size())*dirSplitBoxSize < Magnitude(bbox.Size());
683               
684 
685  for (axis = 0; axis < 5; axis++)
686        if (mInterleaveDirSplits ||
687                (axis < 3 && leaf->depth < mDirSplitDepth) ||
688                (axis >= 3 && leaf->depth >= mDirSplitDepth)
689                ) {
690          if (!mSplitUseOnlyDrivingAxis || axis == sAxis || axis == dAxis) {
691               
692               
693                if (splitType == ESplitRegular) {
694                  if (axis < 3)
695                        nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
696                  else
697                        nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
698                 
699                  nCostRatio[axis] = EvalCostRatio(leaf,
700                                                                                   axis,
701                                                                                   nPosition[axis],
702                                                                                   nRaysBack[axis],
703                                                                                   nRaysFront[axis],
704                                                                                   nPvsBack[axis],
705                                                                                   nPvsFront[axis]
706                                                                                   );
707                } else
708                  if (splitType == ESplitHeuristic) {
709                        nCostRatio[axis] = EvalCostRatioHeuristic(
710                                                                                                          leaf,
711                                                                                                          axis,
712                                                                                                          nPosition[axis],
713                                                                                                          nRaysBack[axis],
714                                                                                                          nRaysFront[axis],
715                                                                                                          nPvsBack[axis],
716                                                                                                          nPvsFront[axis]);
717                  } else
718                        if (splitType == ESplitHybrid) {
719                          if (leaf->depth > 7)
720                                nCostRatio[axis] = EvalCostRatioHeuristic(
721                                                                                                                  leaf,
722                                                                                                                  axis,
723                                                                                                                  nPosition[axis],
724                                                                                                                  nRaysBack[axis],
725                                                                                                                  nRaysFront[axis],
726                                                                                                                  nPvsBack[axis],
727                                                                                                                  nPvsFront[axis]);
728                          else {
729                                if (axis < 3)
730                                  nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
731                                else
732                                  nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
733                               
734                                nCostRatio[axis] = EvalCostRatio(leaf,
735                                                                                                 axis,
736                                                                                                 nPosition[axis],
737                                                                                                 nRaysBack[axis],
738                                                                                                 nRaysFront[axis],
739                                                                                                 nPvsBack[axis],
740                                                                                                 nPvsFront[axis]
741                                                                                                 );
742                          }
743                        } else {
744                          cerr<<"VssTree: Unknown split heuristics\n";
745                          exit(1);
746                        }
747               
748               
749                if ( bestAxis == -1)
750                  bestAxis = axis;
751                else
752                  if ( nCostRatio[axis] < nCostRatio[bestAxis] )
753                        bestAxis = axis;
754          }
755        }
756 
757  axis = bestAxis;
758  position = nPosition[bestAxis];
759
760  raysBack = nRaysBack[bestAxis];
761  raysFront = nRaysFront[bestAxis];
762
763  pvsBack = nPvsBack[bestAxis];
764  pvsFront = nPvsFront[bestAxis];
765       
766  return nCostRatio[bestAxis];
767}
768
769       
770float
771VssTree::EvalCostRatioHeuristic(
772                                                                VssTreeLeaf *leaf,
773                                                                const int axis,
774                                                                float &bestPosition,
775                                                                int &raysBack,
776                                                                int &raysFront,
777                                                                int &pvsBack,
778                                                                int &pvsFront
779                                                                )
780{
781  AxisAlignedBox3 box;
782  float minBox, maxBox;
783       
784  if (axis < 3) {
785        box = GetBBox(leaf);
786        minBox = box.Min(axis);
787        maxBox = box.Max(axis);
788  } else {
789        box = GetDirBBox(leaf);
790        minBox = box.Min(axis-3);
791        maxBox = box.Max(axis-3);
792  }
793       
794  SortSplitCandidates(leaf, axis);
795 
796  // go through the lists, count the number of objects left and right
797  // and evaluate the following cost funcion:
798  // C = ct_div_ci  + (ql*rl + qr*rr)/queries
799       
800  int rl=0, rr = leaf->rays.size();
801  int pl=0, pr = leaf->GetPvsSize();
802  float sizeBox = maxBox - minBox;
803       
804  float minBand = minBox + 0.1*(maxBox - minBox);
805  float maxBand = minBox + 0.9*(maxBox - minBox);
806       
807  float minRatio = 1e20;
808       
809  Intersectable::NewMail();
810  // set all object as belonging to the fron pvs
811  for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
812          ri != leaf->rays.end();
813          ri++)
814        if ((*ri).mRay->IsActive()) {
815          Intersectable *object = (*ri).mRay->mTerminationObject;
816          if (object)
817                if (!object->Mailed()) {
818                  object->Mail();
819                  object->mCounter = 1;
820                } else
821                  object->mCounter++;
822        }
823       
824  Intersectable::NewMail();
825       
826  for(vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
827          ci < splitCandidates->end();
828          ci++) {
829        VssRay *ray;
830        switch ((*ci).type) {
831        case SortableEntry::ERayMin: {
832          rl++;
833          ray = (VssRay *) (*ci).data;
834          Intersectable *object = ray->mTerminationObject;
835          if (object && !object->Mailed()) {
836                object->Mail();
837                pl++;
838          }
839          break;
840        }
841        case SortableEntry::ERayMax: {
842          rr--;
843          ray = (VssRay *) (*ci).data;
844          Intersectable *object = ray->mTerminationObject;
845          if (object) {
846                if (--object->mCounter == 0)
847                  pr--;
848          }
849          break;
850        }
851        }
852
853        float position = (*ci).value;
854               
855        if (position > minBand && position < maxBand) {
856                       
857          float ratio = GetCostRatio(
858                                                                 leaf,
859                                                                 axis,
860                                                                 position,
861                                                                 rl,
862                                                                 rr,
863                                                                 pl,
864                                                                 pr);
865                       
866                       
867          //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
868          //      cout<<"cost= "<<sum<<endl;
869                       
870          if (ratio < minRatio) {
871                minRatio = ratio;
872                bestPosition = position;
873                               
874                raysBack = rl;
875                raysFront = rr;
876                               
877                pvsBack = pl;
878                pvsFront = pr;
879                               
880      }
881    }
882  }
883
884 
885  //  cout<<"===================="<<endl;
886  //  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
887  //      <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
888  return minRatio;
889}
890
891void
892VssTree::SortSplitCandidates(
893                                                         VssTreeLeaf *node,
894                                                         const int axis
895                                                         )
896{
897 
898  splitCandidates->clear();
899 
900  int requestedSize = 2*(node->rays.size());
901  // creates a sorted split candidates array
902  if (splitCandidates->capacity() > 500000 &&
903      requestedSize < (int)(splitCandidates->capacity()/10) ) {
904   
905    delete splitCandidates;
906    splitCandidates = new vector<SortableEntry>;
907  }
908 
909  splitCandidates->reserve(requestedSize);
910
911  // insert all queries
912  for(VssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
913      ri < node->rays.end();
914      ri++) {
915        if ((*ri).mRay->IsActive()) {
916          if (axis < 3) {
917                bool positive = (*ri).mRay->HasPosDir(axis);
918                splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin :
919                                                                                                 SortableEntry::ERayMax,
920                                                                                                 (*ri).ExtrapOrigin(axis),
921                                                                                                 (void *)(*ri).mRay)
922                                                                   );
923                               
924                splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax :
925                                                                                                 SortableEntry::ERayMin,
926                                                                                                 (*ri).ExtrapTermination(axis),
927                                                                                                 (void *)(*ri).mRay)
928                                                                   );
929          } else {
930                float pos = (*ri).mRay->GetDirParametrization(axis-3);
931                splitCandidates->push_back(SortableEntry(SortableEntry::ERayMin,
932                                                                                                 pos - Limits::Small,
933                                                                                                 (void *)(*ri).mRay)
934                                                                   );
935                               
936                splitCandidates->push_back(SortableEntry(SortableEntry::ERayMax,
937                                                                                                 pos + Limits::Small,
938                                                                                                 (void *)(*ri).mRay)
939                                                                   );
940          }
941        }
942  }
943
944  stable_sort(splitCandidates->begin(), splitCandidates->end());
945}
946
947
948void
949VssTree::EvaluateLeafStats(const TraversalData &data)
950{
951
952  // the node became a leaf -> evaluate stats for leafs
953  VssTreeLeaf *leaf = (VssTreeLeaf *)data.node;
954
955  if (data.depth >= termMaxDepth)
956    stat.maxDepthNodes++;
957 
958  //  if ( (int)(leaf->rays.size()) < termMinCost)
959  //    stat.minCostNodes++;
960  if ( leaf->GetPvsSize() < termMinPvs)
961        stat.minPvsNodes++;
962
963  if ( leaf->GetPvsSize() < termMinRays)
964        stat.minRaysNodes++;
965
966  if (0 && leaf->GetAvgRayContribution() > termMaxRayContribution )
967        stat.maxRayContribNodes++;
968       
969  if (SqrMagnitude(data.bbox.Size()) <= termMinSize) {
970        stat.minSizeNodes++;
971  }
972
973  if ( (int)(leaf->rays.size()) > stat.maxRayRefs)
974    stat.maxRayRefs = leaf->rays.size();
975
976}
977
978bool
979VssTree::TerminationCriteriaSatisfied(VssTreeLeaf *leaf)
980{
981  return ( (leaf->GetPvsSize() < termMinPvs) ||
982                   (leaf->rays.size() < termMinRays) ||
983                   //                    (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
984                   (leaf->depth >= termMaxDepth) ||
985                   (SqrMagnitude(GetBBox(leaf).Size()) <= termMinSize)
986                   //              ||
987                   //              (mUseRss && leaf->mPassingRays == leaf->rays.size())
988                   );
989}
990
991
992VssTreeNode *
993VssTree::SubdivideNode(
994                                           VssTreeLeaf *leaf,
995                                           const AxisAlignedBox3 &box,
996                                           AxisAlignedBox3 &backBBox,
997                                           AxisAlignedBox3 &frontBBox
998                                           )
999{
1000 
1001  if (TerminationCriteriaSatisfied(leaf)) {
1002#if 0
1003        if (leaf->depth >= termMaxDepth) {
1004          cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
1005          cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
1006        }
1007#endif
1008               
1009        return leaf;
1010  }
1011       
1012  float position;
1013       
1014  // first count ray sides
1015  int raysBack;
1016  int raysFront;
1017  int pvsBack;
1018  int pvsFront;
1019       
1020  // select subdivision axis
1021  int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
1022  //    cout<<axis<<" ";
1023       
1024  //    cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
1025  //            pvsFront<<endl;
1026
1027  if (axis == -1) {
1028    return leaf;
1029  }
1030 
1031  stat.nodes+=2;
1032  stat.splits[axis]++;
1033
1034  // add the new nodes to the tree
1035  VssTreeInterior *node = new VssTreeInterior(leaf->parent);
1036
1037  node->axis = axis;
1038  node->position = position;
1039  node->bbox = box;
1040  node->dirBBox = GetDirBBox(leaf);
1041 
1042  backBBox = box;
1043  frontBBox = box;
1044
1045  VssTreeLeaf *back = new VssTreeLeaf(node, raysBack);
1046  VssTreeLeaf *front = new VssTreeLeaf(node, raysFront);
1047
1048  // replace a link from node's parent
1049  if (  leaf->parent )
1050    leaf->parent->ReplaceChildLink(leaf, node);
1051  // and setup child links
1052  node->SetupChildLinks(back, front);
1053       
1054  if (axis <= VssTreeNode::SPLIT_Z) {
1055        backBBox.SetMax(axis, position);
1056    frontBBox.SetMin(axis, position);
1057               
1058        for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1059                ri != leaf->rays.end();
1060                ri++) {
1061      if ((*ri).mRay->IsActive()) {
1062                               
1063                // first unref ray from the former leaf
1064                (*ri).mRay->Unref();
1065
1066                // Debug << "computed t: " << (*ri).mRay->mT << endl;
1067                // determine the side of this ray with respect to the plane
1068                int side = node->ComputeRayIntersection(*ri, (*ri).mRay->mT);
1069               
1070                if (side == 0) {
1071                  if ((*ri).mRay->HasPosDir(axis)) {
1072                        back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
1073                                                                                          (*ri).mMinT,
1074                                                                                          (*ri).mRay->mT)
1075                                                 );
1076                        front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
1077                                                                                           (*ri).mRay->mT,
1078                                                                                           (*ri).mMaxT));
1079                  } else {
1080                        back->AddRay(VssTreeNode::RayInfo((*ri).mRay,
1081                                                                                          (*ri).mRay->mT,
1082                                                                                          (*ri).mMaxT));
1083                        front->AddRay(VssTreeNode::RayInfo((*ri).mRay,
1084                                                                                           (*ri).mMinT,
1085                                                                                           (*ri).mRay->mT));
1086                  }
1087                } else
1088                  if (side == 1)
1089                        front->AddRay(*ri);
1090                  else
1091                        back->AddRay(*ri);
1092      } else
1093                (*ri).mRay->Unref();
1094    }
1095  } else {
1096    // rays front/back
1097   
1098   
1099    for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1100                ri != leaf->rays.end();
1101                ri++) {
1102      if ((*ri).mRay->IsActive()) {
1103                // first unref ray from the former leaf
1104                (*ri).mRay->Unref();
1105
1106                int side;
1107                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
1108                  side = 1;
1109                else
1110                  side = -1;
1111                               
1112                if (side == 1)
1113                  front->AddRay(*ri);
1114                else
1115                  back->AddRay(*ri);
1116                               
1117      } else
1118                (*ri).mRay->Unref();
1119    }
1120  }
1121
1122  front->SetPvsSize(pvsFront);
1123  back->SetPvsSize(pvsBack);
1124  // compute entropy as well
1125  front->ComputeEntropyImportance();
1126  back->ComputeEntropyImportance();
1127 
1128  // update stats
1129  stat.rayRefs -= leaf->rays.size();
1130  stat.rayRefs += raysBack + raysFront;
1131
1132 
1133  delete leaf;
1134  return node;
1135}
1136
1137
1138
1139
1140
1141
1142int
1143VssTree::ReleaseMemory(const int time)
1144{
1145  stack<VssTreeNode *> tstack;
1146 
1147  // find a node in the tree which subtree will be collapsed
1148  int maxAccessTime = time - accessTimeThreshold;
1149  int released;
1150
1151  tstack.push(root);
1152
1153  while (!tstack.empty()) {
1154    VssTreeNode *node = tstack.top();
1155    tstack.pop();
1156   
1157 
1158    if (!node->IsLeaf()) {
1159      VssTreeInterior *in = (VssTreeInterior *)node;
1160      //      cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
1161      if (in->depth >= minCollapseDepth &&
1162                  in->lastAccessTime <= maxAccessTime) {
1163                released = CollapseSubtree(node, time);
1164                break;
1165      }
1166     
1167      if (in->back->GetAccessTime() < 
1168                  in->front->GetAccessTime()) {
1169                tstack.push(in->front);
1170                tstack.push(in->back);
1171      } else {
1172                tstack.push(in->back);
1173                tstack.push(in->front);
1174      }
1175    }
1176  }
1177
1178  while (tstack.empty()) {
1179    // could find node to collaps...
1180    //    cout<<"Could not find a node to release "<<endl;
1181    break;
1182  }
1183 
1184  return released;
1185}
1186
1187
1188
1189
1190VssTreeNode *
1191VssTree::SubdivideLeaf(
1192                                           VssTreeLeaf *leaf
1193                                           )
1194{
1195  VssTreeNode *node = leaf;
1196       
1197  AxisAlignedBox3 leafBBox = GetBBox(leaf);
1198
1199  static int pass = 0;
1200  pass ++;
1201       
1202  // check if we should perform a dynamic subdivision of the leaf
1203  if (!TerminationCriteriaSatisfied(leaf)) {
1204   
1205        // memory check and realese...
1206    if (GetMemUsage() > maxTotalMemory) {
1207      ReleaseMemory( pass );
1208    }
1209   
1210    AxisAlignedBox3 backBBox, frontBBox;
1211
1212    // subdivide the node
1213    node =
1214      SubdivideNode(leaf,
1215                                        leafBBox,
1216                                        backBBox,
1217                                        frontBBox
1218                                        );
1219  }
1220       
1221  return node;
1222}
1223
1224
1225
1226void
1227VssTree::UpdateRays(VssRayContainer &remove,
1228                                        VssRayContainer &add
1229                                        )
1230{
1231  VssTreeLeaf::NewMail();
1232
1233  // schedule rays for removal
1234  for(VssRayContainer::const_iterator ri = remove.begin();
1235      ri != remove.end();
1236      ri++) {
1237    (*ri)->ScheduleForRemoval();
1238  }
1239
1240  int inactive=0;
1241
1242  for(VssRayContainer::const_iterator ri = remove.begin();
1243      ri != remove.end();
1244      ri++) {
1245    if ((*ri)->ScheduledForRemoval())
1246          //      RemoveRay(*ri, NULL, false);
1247          // !!! BUG - with true it does not work correctly - aggreated delete
1248      RemoveRay(*ri, NULL, true);
1249    else
1250      inactive++;
1251  }
1252
1253
1254  //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1255 
1256  for(VssRayContainer::const_iterator ri = add.begin();
1257      ri != add.end();
1258      ri++) {
1259        VssTreeNode::RayInfo info(*ri);
1260        if (ClipRay(info, bbox))
1261          AddRay(info);
1262  }
1263}
1264
1265
1266void
1267VssTree::RemoveRay(VssRay *ray,
1268                                   vector<VssTreeLeaf *> *affectedLeaves,
1269                                   const bool removeAllScheduledRays
1270                                   )
1271{
1272       
1273  stack<RayTraversalData> tstack;
1274
1275  tstack.push(RayTraversalData(root, VssTreeNode::RayInfo(ray)));
1276 
1277  RayTraversalData data;
1278
1279  // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1280
1281  while (!tstack.empty()) {
1282    data = tstack.top();
1283    tstack.pop();
1284
1285    if (!data.node->IsLeaf()) {
1286      // split the set of rays in two groups intersecting the
1287      // two subtrees
1288
1289      TraverseInternalNode(data, tstack);
1290     
1291    } else {
1292      // remove the ray from the leaf
1293      // find the ray in the leaf and swap it with the last ray...
1294      VssTreeLeaf *leaf = (VssTreeLeaf *)data.node;
1295     
1296      if (!leaf->Mailed()) {
1297                leaf->Mail();
1298                if (affectedLeaves)
1299                  affectedLeaves->push_back(leaf);
1300       
1301                if (removeAllScheduledRays) {
1302                  int tail = leaf->rays.size()-1;
1303
1304                  for (int i=0; i < (int)(leaf->rays.size()); i++) {
1305                        if (leaf->rays[i].mRay->ScheduledForRemoval()) {
1306                          // find a ray to replace it with
1307                          while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) {
1308                                stat.removedRayRefs++;
1309                                leaf->rays[tail].mRay->Unref();
1310                                leaf->rays.pop_back();
1311                                tail--;
1312                          }
1313
1314                          if (tail < i)
1315                                break;
1316             
1317                          stat.removedRayRefs++;
1318                          leaf->rays[i].mRay->Unref();
1319                          leaf->rays[i] = leaf->rays[tail];
1320                          leaf->rays.pop_back();
1321                          tail--;
1322                        }
1323                  }
1324                }
1325      }
1326     
1327      if (!removeAllScheduledRays)
1328                for (int i=0; i < (int)leaf->rays.size(); i++) {
1329                  if (leaf->rays[i].mRay == ray) {
1330                        stat.removedRayRefs++;
1331                        ray->Unref();
1332                        leaf->rays[i] = leaf->rays[leaf->rays.size()-1];
1333                        leaf->rays.pop_back();
1334                        // check this ray again
1335                        break;
1336                  }
1337                }
1338     
1339    }
1340  }
1341
1342  if (ray->RefCount() != 0) {
1343    cerr<<"Error: Number of remaining refs = "<<ray->RefCount()<<endl;
1344    exit(1);
1345  }
1346 
1347}
1348
1349
1350void
1351VssTree::AddRay(VssTreeNode::RayInfo &info)
1352{
1353
1354  stack<RayTraversalData> tstack;
1355 
1356  tstack.push(RayTraversalData(root, info));
1357 
1358  RayTraversalData data;
1359 
1360  while (!tstack.empty()) {
1361    data = tstack.top();
1362    tstack.pop();
1363
1364    if (!data.node->IsLeaf()) {
1365      TraverseInternalNode(data, tstack);
1366    } else {
1367      // remove the ray from the leaf
1368      // find the ray in the leaf and swap it with the last ray...
1369      VssTreeLeaf *leaf = (VssTreeLeaf *)data.node;
1370      leaf->AddRay(data.rayData);
1371      stat.addedRayRefs++;
1372    }
1373  }
1374}
1375
1376void
1377VssTree::TraverseInternalNode(
1378                                                          RayTraversalData &data,
1379                                                          stack<RayTraversalData> &tstack)
1380{
1381  VssTreeInterior *in =  (VssTreeInterior *) data.node;
1382 
1383  if (in->axis <= VssTreeNode::SPLIT_Z) {
1384
1385    // determine the side of this ray with respect to the plane
1386    int side = in->ComputeRayIntersection(data.rayData,
1387                                                                                  data.rayData.mRay->mT);
1388   
1389 
1390    if (side == 0) {
1391      if (data.rayData.mRay->HasPosDir(in->axis)) {
1392                tstack.push(RayTraversalData(in->back,
1393                                                                         VssTreeNode::RayInfo(data.rayData.mRay,
1394                                                                                                                  data.rayData.mMinT,
1395                                                                                                                  data.rayData.mRay->mT))
1396                                        );
1397                               
1398                tstack.push(RayTraversalData(in->front,
1399                                                                         VssTreeNode::RayInfo(data.rayData.mRay,
1400                                                                                                                  data.rayData.mRay->mT,
1401                                                                                                                  data.rayData.mMaxT
1402                                                                                                                  ))
1403                                        );
1404       
1405      } else {
1406                tstack.push(RayTraversalData(in->back,
1407                                                                         VssTreeNode::RayInfo(data.rayData.mRay,
1408                                                                                                                  data.rayData.mRay->mT,
1409                                                                                                                  data.rayData.mMaxT
1410                                                                                                                  ))
1411                                        );
1412                               
1413                tstack.push(RayTraversalData(in->front,
1414                                                                         VssTreeNode::RayInfo(data.rayData.mRay,
1415                                                                                                                  data.rayData.mMinT,
1416                                                                                                                  data.rayData.mRay->mT))
1417                                        );
1418                               
1419                               
1420      }
1421    } else
1422      if (side == 1)
1423                tstack.push(RayTraversalData(in->front, data.rayData));
1424      else
1425                tstack.push(RayTraversalData(in->back, data.rayData));
1426  }
1427  else {
1428    // directional split
1429        if (data.rayData.mRay->GetDirParametrization(in->axis - 3) > in->position)
1430          tstack.push(RayTraversalData(in->front, data.rayData));
1431        else
1432          tstack.push(RayTraversalData(in->back, data.rayData));
1433  }
1434}
1435
1436
1437int
1438VssTree::CollapseSubtree(VssTreeNode *sroot, const int time)
1439{
1440  // first count all rays in the subtree
1441  // use mail 1 for this purpose
1442  stack<VssTreeNode *> tstack;
1443  int rayCount = 0;
1444  int totalRayCount = 0;
1445  int collapsedNodes = 0;
1446
1447#if DEBUG_COLLAPSE
1448  cout<<"Collapsing subtree"<<endl;
1449  cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
1450  cout<<"depth="<<(int)sroot->depth<<endl;
1451#endif
1452
1453  //    tstat.collapsedSubtrees++;
1454  //    tstat.collapseDepths += (int)sroot->depth;
1455  //    tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1456 
1457  tstack.push(sroot);
1458  VssRay::NewMail();
1459       
1460  while (!tstack.empty()) {
1461    collapsedNodes++;
1462    VssTreeNode *node = tstack.top();
1463    tstack.pop();
1464
1465    if (node->IsLeaf()) {
1466      VssTreeLeaf *leaf = (VssTreeLeaf *) node;
1467      for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1468                  ri != leaf->rays.end();
1469                  ri++) {
1470                               
1471                totalRayCount++;
1472                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed()) {
1473                  (*ri).mRay->Mail();
1474                  rayCount++;
1475                }
1476      }
1477    } else {
1478      tstack.push(((VssTreeInterior *)node)->back);
1479      tstack.push(((VssTreeInterior *)node)->front);
1480    }
1481  }
1482 
1483  VssRay::NewMail();
1484
1485  // create a new node that will hold the rays
1486  VssTreeLeaf *newLeaf = new VssTreeLeaf( sroot->parent, rayCount );
1487  if (  newLeaf->parent )
1488    newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
1489 
1490
1491  tstack.push( sroot );
1492
1493  while (!tstack.empty()) {
1494
1495    VssTreeNode *node = tstack.top();
1496    tstack.pop();
1497
1498    if (node->IsLeaf()) {
1499      VssTreeLeaf *leaf = (VssTreeLeaf *) node;
1500     
1501      for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1502                  ri != leaf->rays.end();
1503                  ri++) {
1504                               
1505                // unref this ray from the old node
1506                               
1507                if ((*ri).mRay->IsActive()) {
1508                  (*ri).mRay->Unref();
1509                  if (!(*ri).mRay->Mailed()) {
1510                        (*ri).mRay->Mail();
1511                        newLeaf->AddRay(*ri);
1512                  }
1513                } else
1514                  (*ri).mRay->Unref();
1515                               
1516      }
1517    } else {
1518      tstack.push(((VssTreeInterior *)node)->back);
1519      tstack.push(((VssTreeInterior *)node)->front);
1520    }
1521  }
1522 
1523  // delete the node and all its children
1524  delete sroot;
1525 
1526  //   for(VssTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
1527  //       ri != newLeaf->rays.end();
1528  //       ri++)
1529  //     (*ri).ray->UnMail(2);
1530
1531
1532#if DEBUG_COLLAPSE
1533  cout<<"Total memory before="<<GetMemUsage()<<endl;
1534#endif
1535
1536  stat.nodes -= collapsedNodes - 1;
1537  stat.rayRefs -= totalRayCount - rayCount;
1538 
1539#if DEBUG_COLLAPSE
1540  cout<<"collapsed nodes"<<collapsedNodes<<endl;
1541  cout<<"collapsed rays"<<totalRayCount - rayCount<<endl;
1542  cout<<"Total memory after="<<GetMemUsage()<<endl;
1543  cout<<"================================"<<endl;
1544#endif
1545
1546  //  tstat.collapsedNodes += collapsedNodes;
1547  //  tstat.collapsedRays += totalRayCount - rayCount;
1548   
1549  return totalRayCount - rayCount;
1550}
1551
1552
1553int
1554VssTree::GetPvsSize(const AxisAlignedBox3 &box) const
1555{
1556  stack<VssTreeNode *> tstack;
1557  tstack.push(root);
1558
1559  Intersectable::NewMail();
1560  int pvsSize = 0;
1561       
1562  while (!tstack.empty()) {
1563    VssTreeNode *node = tstack.top();
1564    tstack.pop();
1565   
1566 
1567    if (node->IsLeaf()) {
1568          VssTreeLeaf *leaf = (VssTreeLeaf *)node;
1569          for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1570                  ri != leaf->rays.end();
1571                  ri++)
1572                if ((*ri).mRay->IsActive()) {
1573                  Intersectable *object;
1574#if BIDIRECTIONAL_RAY
1575                  object = (*ri).mRay->mOriginObject;
1576                  if (object && !object->Mailed()) {
1577                        pvsSize++;
1578                        object->Mail();
1579                  }
1580#endif
1581                  object = (*ri).mRay->mTerminationObject;
1582                  if (object && !object->Mailed()) {
1583                        pvsSize++;
1584                        object->Mail();
1585                  }
1586                }
1587        } else {
1588          VssTreeInterior *in = (VssTreeInterior *)node;
1589          if (in->axis < 3) {
1590                if (box.Max(in->axis) >= in->position )
1591                  tstack.push(in->front);
1592                               
1593                if (box.Min(in->axis) <= in->position )
1594                  tstack.push(in->back);
1595          } else {
1596                // both nodes for directional splits
1597                tstack.push(in->front);
1598                tstack.push(in->back);
1599          }
1600        }
1601  }
1602  return pvsSize;
1603}
1604
1605void
1606VssTree::GetRayContributionStatistics(
1607                                                                          float &minRayContribution,
1608                                                                          float &maxRayContribution,
1609                                                                          float &avgRayContribution
1610                                                                          )
1611{
1612  stack<VssTreeNode *> tstack;
1613  tstack.push(root);
1614       
1615  minRayContribution = 1.0f;
1616  maxRayContribution = 0.0f;
1617  float sumRayContribution = 0.0f;
1618  int leaves = 0;
1619
1620  while (!tstack.empty()) {
1621    VssTreeNode *node = tstack.top();
1622    tstack.pop();
1623               
1624    if (node->IsLeaf()) {
1625          leaves++;
1626          VssTreeLeaf *leaf = (VssTreeLeaf *)node;
1627          float c = leaf->GetImportance();
1628          if (c > maxRayContribution)
1629                maxRayContribution = c;
1630          if (c < minRayContribution)
1631                minRayContribution = c;
1632          sumRayContribution += c;
1633                       
1634        } else {
1635          VssTreeInterior *in = (VssTreeInterior *)node;
1636          // both nodes for directional splits
1637          tstack.push(in->front);
1638          tstack.push(in->back);
1639        }
1640  }
1641       
1642  cout<<"sum="<<sumRayContribution<<endl;
1643  cout<<"leaves="<<leaves<<endl;
1644  avgRayContribution = sumRayContribution/(float)leaves;
1645}
1646
1647
1648int
1649VssTree::GenerateRays(const float ratioPerLeaf,
1650                                          SimpleRayContainer &rays)
1651{
1652  stack<VssTreeNode *> tstack;
1653  tstack.push(root);
1654       
1655  while (!tstack.empty()) {
1656    VssTreeNode *node = tstack.top();
1657    tstack.pop();
1658               
1659    if (node->IsLeaf()) {
1660          VssTreeLeaf *leaf = (VssTreeLeaf *)node;
1661          float c = leaf->GetImportance();
1662          int num = (c*ratioPerLeaf + 0.5);
1663          //                    cout<<num<<" ";
1664
1665          for (int i=0; i < num; i++) {
1666                Vector3 origin = GetBBox(leaf).GetRandomPoint();
1667                Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1668                Vector3 direction = VssRay::GetDirection(dirVector.x, dirVector.y);
1669                //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1670                rays.push_back(SimpleRay(origin, direction));
1671          }
1672                       
1673        } else {
1674          VssTreeInterior *in = (VssTreeInterior *)node;
1675          // both nodes for directional splits
1676          tstack.push(in->front);
1677          tstack.push(in->back);
1678        }
1679  }
1680
1681  return rays.size();
1682}
1683
1684void
1685VssTree::CollectLeaves(vector<VssTreeLeaf *> &leaves)
1686{
1687  stack<VssTreeNode *> tstack;
1688  tstack.push(root);
1689       
1690  while (!tstack.empty()) {
1691    VssTreeNode *node = tstack.top();
1692    tstack.pop();
1693               
1694    if (node->IsLeaf()) {
1695          VssTreeLeaf *leaf = (VssTreeLeaf *)node;
1696          leaves.push_back(leaf);
1697        } else {
1698          VssTreeInterior *in = (VssTreeInterior *)node;
1699          // both nodes for directional splits
1700          tstack.push(in->front);
1701          tstack.push(in->back);
1702        }
1703  }
1704}
1705
1706bool
1707VssTree::ValidLeaf(VssTreeLeaf *leaf) const
1708{
1709  return leaf->rays.size() > termMinRays/4;
1710}
1711
1712
1713void
1714VssTree::GenerateLeafRays(VssTreeLeaf *leaf,
1715                                                  const int numberOfRays,
1716                                                  SimpleRayContainer &rays)
1717{
1718  int nrays = leaf->rays.size();
1719  for (int i=0; i < numberOfRays; i++) {
1720        // pickup 3 random rays
1721        int r1 = Random(nrays-1);
1722        int r2 = Random(nrays-1);
1723        int r3 = Random(nrays-1);
1724               
1725        Vector3 o1 = leaf->rays[r1].Extrap(RandomValue(leaf->rays[r1].GetMinT(),
1726                                                                                                   leaf->rays[r1].GetMaxT()));
1727
1728        Vector3 o2 = leaf->rays[r2].Extrap(RandomValue(leaf->rays[r2].GetMinT(),
1729                                                                                                   leaf->rays[r2].GetMaxT()));
1730               
1731        Vector3 o3 = leaf->rays[r3].Extrap(RandomValue(leaf->rays[r3].GetMinT(),
1732                                                                                                   leaf->rays[r3].GetMaxT()));
1733
1734        const float overlap = 0.1;
1735               
1736        Vector3 origin, direction;
1737        bool useExtendedConvexCombination = true;
1738        if (useExtendedConvexCombination) {
1739          float w1, w2, w3;
1740          GenerateExtendedConvexCombinationWeights(w1, w2, w3, overlap);
1741          origin = w1*o1 + w2*o2 + w3*o3;
1742          direction =
1743                w1*leaf->rays[r1].mRay->GetDir() +
1744                w2*leaf->rays[r2].mRay->GetDir() +
1745                w3*leaf->rays[r3].mRay->GetDir();
1746        } else {
1747          origin = GetBBox(leaf).GetRandomPoint();
1748          Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1749          direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1750        }
1751        //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1752        rays.push_back(SimpleRay(origin, direction));
1753  }
1754}
1755
1756int
1757VssTree::GenerateRays(const int numberOfRays,
1758                                          const int numberOfLeaves,
1759                                          SimpleRayContainer &rays)
1760{
1761
1762  vector<VssTreeLeaf *> leaves;
1763       
1764  CollectLeaves(leaves);
1765       
1766  sort(leaves.begin(),
1767           leaves.end(),
1768           GreaterContribution);
1769                         
1770
1771  float sumContrib = 0.0;
1772  int i;
1773  int k = 0;
1774  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
1775        if (ValidLeaf(leaves[i])) {
1776          float c = leaves[i]->GetImportance();
1777          sumContrib += c;
1778          //            cout<<"ray contrib "<<i<<" : "<<c<<endl;
1779          k++;
1780        }
1781                               
1782  float avgContrib = sumContrib/numberOfLeaves;
1783  float ratioPerLeaf = numberOfRays/(avgContrib*numberOfLeaves);
1784  k = 0;
1785  for (i=0; i < leaves.size() && k < numberOfLeaves; i++)
1786        if (ValidLeaf(leaves[i])) {
1787          k++;
1788          VssTreeLeaf *leaf = leaves[i];
1789          float c = leaf->GetImportance();
1790          int num = (c*ratioPerLeaf + 0.5);
1791          GenerateLeafRays(leaf, num, rays);
1792        }
1793
1794  return rays.size();
1795}
1796
1797
1798float
1799VssTree::GetAvgPvsSize()
1800{
1801  stack<VssTreeNode *> tstack;
1802  tstack.push(root);
1803
1804  int sumPvs = 0;
1805  int leaves = 0;
1806  while (!tstack.empty()) {
1807    VssTreeNode *node = tstack.top();
1808    tstack.pop();
1809               
1810    if (node->IsLeaf()) {
1811          VssTreeLeaf *leaf = (VssTreeLeaf *)node;
1812          // update pvs size
1813          leaf->UpdatePvsSize();
1814          sumPvs += leaf->GetPvsSize();
1815          leaves++;
1816        } else {
1817          VssTreeInterior *in = (VssTreeInterior *)node;
1818          // both nodes for directional splits
1819          tstack.push(in->front);
1820          tstack.push(in->back);
1821        }
1822  }
1823
1824
1825  return sumPvs/(float)leaves;
1826}
1827
1828       
1829float
1830VssTreeLeaf::GetImportance()  const
1831{
1832
1833  if (1) {
1834        //      return GetAvgRayContribution();
1835        return GetPvsSize();
1836  } else {
1837        // return GetAvgRayContribution()*mEntropyImportance;
1838        //return GetAvgRayContribution();
1839        return mEntropyImportance;
1840  }
1841}
1842
1843
1844float
1845VssTreeLeaf::ComputePvsEntropy()
1846{
1847  int samples = 0;
1848  Intersectable::NewMail();
1849  // set all object as belonging to the fron pvs
1850  for(VssTreeNode::RayInfoContainer::iterator ri = rays.begin();
1851          ri != rays.end();
1852          ri++)
1853        if ((*ri).mRay->IsActive()) {
1854          Intersectable *object = (*ri).mRay->mTerminationObject;
1855          if (object) {
1856                if (!object->Mailed()) {
1857                  object->Mail();
1858                  object->mCounter = 1;
1859                } else
1860                  object->mCounter++;
1861                samples++;
1862          }
1863        }
1864 
1865  float entropy = 0.0f;
1866 
1867  if (samples > 1) {
1868        Intersectable::NewMail();
1869        for(RayInfoContainer::const_iterator ri = rays.begin();
1870                ri != rays.end();
1871                ri++)
1872          if ((*ri).mRay->IsActive()) {
1873                Intersectable *object = (*ri).mRay->mTerminationObject;
1874                if (object) {
1875                  if (!object->Mailed()) {
1876                        object->Mail();
1877                        float p = object->mCounter/(float)samples;
1878                        entropy -= p*log(p);
1879                  }
1880                }
1881          }
1882        entropy = entropy/log((float)samples);
1883  }
1884  else
1885        entropy = 1.0f;
1886 
1887  return entropy;
1888}
1889
1890float
1891VssTreeLeaf::ComputeRayLengthEntropy()
1892{
1893  // get sum of all ray lengths
1894  // consider only passing rays or originating rays
1895  float sum = 0.0f;
1896  int samples = 0;
1897
1898  for(RayInfoContainer::const_iterator ri = rays.begin();
1899      ri != rays.end();
1900      ri++) {
1901        int rayClass = (*ri).GetRayClass();
1902        if (
1903                rayClass == RayInfo::PASSING_RAY
1904                //              ||
1905                //              rayClass == RayInfo::SOURCE_RAY
1906                //              rayClass == RayInfo::TERMINATION_RAY
1907                ) {
1908          float c = 1.0f - (*ri).GetMaxT();
1909          sum += (*ri).mRay->GetSize()*c;
1910          samples++;
1911        }
1912  }
1913 
1914  float entropy = 0.0f;
1915 
1916  if (samples > 1) {
1917        for(RayInfoContainer::const_iterator ri = rays.begin();
1918                ri != rays.end();
1919                ri++) {
1920          int rayClass = (*ri).GetRayClass();
1921          if (
1922                  rayClass == RayInfo::PASSING_RAY
1923                  //              ||
1924                  //              rayClass == RayInfo::SOURCE_RAY
1925                  // rayClass == RayInfo::TERMINATION_RAY
1926                  ) {
1927                float c = 1.0f - (*ri).GetMaxT();
1928                float p = (*ri).mRay->GetSize()*c/sum;
1929                entropy -= p*log(p);
1930          }
1931        }
1932        entropy = entropy/log((float)samples);
1933  } else
1934        entropy = 1.0f;
1935 
1936  return entropy;
1937}
1938 
1939float
1940VssTreeLeaf::ComputeRayTerminationEntropy()
1941{
1942
1943  return 0.0f;
1944}
1945
1946
1947void
1948VssTreeLeaf::ComputeEntropyImportance()
1949{
1950  //  mEntropy = 1.0f - ComputeRayLengthEntropy();
1951  mEntropyImportance = 1.0f - ComputePvsEntropy();
1952 
1953  //  cout<<"ei="<<mEntropyImportance<<" ";
1954}
Note: See TracBrowser for help on using the repository browser.