source: GTP/trunk/Lib/Vis/Preprocessing/src/VssTree.cpp @ 2575

Revision 2575, 47.3 KB checked in by bittner, 17 years ago (diff)

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