source: GTP/trunk/Lib/Vis/Preprocessing/src/havran/ktbs.cpp @ 2632

Revision 2632, 50.6 KB checked in by bittner, 17 years ago (diff)
Line 
1// ===================================================================
2// $Id: $
3//
4// ktbs.cpp
5//     The classes for building up CKTB tree by statistical optimization.
6//     This building up is based on sampling.
7//
8// REPLACEMENT_STRING
9//
10// Copyright by Vlastimil Havran, 2007 - email to "vhavran AT seznam.cz"
11// Initial coding by Vlasta Havran, January 2008.
12
13// standard headers
14#include <algorithm>
15#include <iostream>
16#include <string>
17
18// GOLEM headers
19#include "ktbs.h"
20#include "timer.h"
21#include "Environment.h"
22
23//#define _DEBUG
24
25#define AVOIDCHECKLIST
26
27#define _USE_OPTIMIZE_DIVIDE_AX
28
29// Note that the RADIX SORT does not work properly with -mtune=pentium3 for GCC.4.0 !!!!
30// Note that the RADIX SORT does not work properly with -mtune=pentium4 for GCC.3.4 !!!!
31// The reason is unknown, probably the bug in compiler
32
33namespace GtpVisibilityPreprocessor {
34
35// ------------------------------------------------
36// for debugging cost function in the scene
37//#define _DEBUG_COSTFUNCTION
38#ifdef _DEBUG_COSTFUNCTION
39
40const int indexToFindSS = 0;
41
42static CFileIO ffcost;
43static int     costEvalCnt;
44static int     costCntObjects;
45static int     costIndexDebug;
46static string  fnameCost;
47static float   maxCost;
48static float   minCost;
49static float   maxPos;
50static float   minPos;
51static float   minPosAxis;
52static float   maxPosAxis;
53static bool _alreadyDebugged = false;
54static bool _streamOpened = false;
55
56static void
57InitCostStream(int indexDebug, int cntObjects,
58               float minx, float maxx)
59{
60  maxCost = -INFINITY;
61  minCost = INFINITY;
62  costEvalCnt = 0;
63  costIndexDebug = indexDebug;
64  costCntObjects = cntObjects;
65  minPosAxis = minx;
66  maxPosAxis = maxx;
67 
68  // and init the stream - rewrite
69  string name = "debugcostfile"
70  //Environment::GetSingleton()->GetBoolValue("OutputFileName", name);
71  char lns[100];
72  sprintf(lns, "%05d.res", indexDebug);
73  ChangeFilenameExtension(name, string(lns), fnameCost);
74
75  ffcost.SetFilename(fnameCost.c_str());
76  if (ffcost.OpenInMode(CFileIO::EE_Mode_w)) {
77    cerr << "Opening of stream failed" << endl;
78    abort();;
79  }
80  cout << "Saving debug to " << fnameCost << endl;
81  ffcost.WriteChars("\n");
82  sprintf(lns, "#CntObjects = %d\n", costCntObjects);
83  ffcost.WriteChars(lns);
84  sprintf(lns, "#IndexToDebug = %d\n", costIndexDebug);
85  ffcost.WriteChars(lns);
86  ffcost.WriteChars("#File ");
87  ffcost.WriteChars(fnameCost.c_str());
88  ffcost.WriteChars("\n");
89  ffcost.WriteChars("#==============================================\n");
90
91  cout << "#---------------------------------------" << endl;
92  _streamOpened = true;
93}
94
95static void
96CloseCostStream()
97{
98  if (_streamOpened) {   
99    char lns[255];
100    ffcost.WriteChars("#==============================================\n");
101    ffcost.WriteChars("#EndOfFile\n");
102    sprintf(lns, "#minCost = %8.6f minPos = %8.6f\n", minCost, minPos);
103    ffcost.WriteChars(lns);   
104    sprintf(lns, "#maxCost = %8.6f maxPos = %8.6f\n", maxCost, maxPos);
105    ffcost.WriteChars(lns);   
106    sprintf(lns, "#ratioMin/Max =                                %6.5f\n",
107            minCost/maxCost);
108    ffcost.WriteChars(lns);   
109    sprintf(lns, "#CntEval = %d\n", costEvalCnt);
110    ffcost.WriteChars(lns);
111    sprintf(lns, "#CntObjects = %d\n", costCntObjects);
112    ffcost.WriteChars(lns);   
113    sprintf(lns, "#IndexToDebug = %d\n", costIndexDebug);
114    ffcost.WriteChars(lns);
115    sprintf(lns, "#MinPosAxis = %8.6f maxPosAxis = %8.6f\n", minPosAxis, maxPosAxis);
116    ffcost.WriteChars(lns);
117    ffcost.WriteChars("#File ");
118    ffcost.WriteChars(fnameCost.c_str());
119    ffcost.WriteChars("\n");
120    ffcost.Close();
121    //cout << "#=======================================" << endl;
122    _streamOpened = false;
123    _alreadyDebugged = true;
124    // Temporarily
125    exit(0);
126  }
127}
128
129static void
130ReportCostStream(float pos, float cost)
131{
132  // Here we normalize a position to the interval <0-1>
133  float posn = (pos - minPosAxis)/(maxPosAxis - minPosAxis);
134
135  if (_streamOpened) {
136    costEvalCnt++;
137    if (cost < minCost) {
138      minCost = cost;
139      minPos = pos;
140    }
141    if (cost > maxCost) {
142      maxCost = cost;
143      maxPos = pos;
144    }
145    char lns[255];
146    sprintf(lns, "%8.6f %8.6f\n", posn, cost);
147    //cout << lns << endl;
148    ffcost.WriteChars(lns);   
149  }
150}
151#endif // _DEBUG_COSTFUNCTION
152
153// -------------------------------------------------------
154// the next two axes for a given one
155const CKTBAxes::Axes
156CKTBSBuildUp::oaxes[3][2]=
157{{CKTBAxes::EE_Y_axis, CKTBAxes::EE_Z_axis},
158 {CKTBAxes::EE_Z_axis, CKTBAxes::EE_X_axis},
159 {CKTBAxes::EE_X_axis, CKTBAxes::EE_Y_axis}};
160 
161//----------------------------------------------------------------
162//---------- Implementation of CKTB tree with irregular -----------
163//---------- placed splitting planes -----------------------------
164//----------------------------------------------------------------
165
166// default constructor
167CKTBSBuildUp::CKTBSBuildUp()
168{
169  // verbose is set
170  verbose = true;
171
172  // Maximum depth of the tree is set and stack is allocated
173  maxTreeDepth = MAX_HEIGHT;
174  stackDepth = maxTreeDepth + 2;
175  stackID = new GALIGN16 SInputData[stackDepth];
176  assert(stackID);
177  stackIndex = 0;
178  return;
179}
180
181//  destructor
182CKTBSBuildUp::~CKTBSBuildUp()
183{
184  delete [] stackID;
185  stackID = 0;
186}
187
188void
189CKTBSBuildUp::ProvideID(ostream &app)
190{
191  bool tempvar;
192 
193  string termCrit;
194  Environment::GetSingleton()->GetStringValue("BSP.termCrit", termCrit);
195  app << "#BSP.termCrit\t\tTermination criteria to build up BSP tree\n";
196  app << termCrit << "\n";
197
198  maxCountTrials = 0;
199  if ( (!termCrit.compare("auto")) &&
200       (!termCrit.compare("auto2")) ) {
201    // everything except auto settings, when depth is determined
202    // It should be reported for adhoc and possibly others
203    Environment::GetSingleton()->GetIntValue("BSP.maxDepthAllowed",
204                                             maxDepthAllowed);
205    Environment::GetSingleton()->GetIntValue("BSP.maxListLength",
206                                             maxListLength);
207    app << "#MaxDepth (Dmax)\tMaximal allowed depth of the BSP tree\n";
208    app << maxDepthAllowed << "\n";
209    app << "#MaxListLength (Noilf)\tMaximal number of solids in one cell\n";
210    app << maxListLength << "\n";
211    // maximum number of trials to further subdivide
212    maxCountTrials = (int)(1.0 + 0.2 * (float)maxDepthAllowed);
213    if (maxCountTrials < 3)
214      maxCountTrials = 3;
215  }
216  else {
217    if (!termCrit.compare("auto2") ) {
218      Environment::GetSingleton()->GetIntValue("BSP.maxListLength",
219                                               maxListLength);
220      app << "#MaxDepth (Dmax)\tMaximal allowed depth of the BSP tree\n";
221    }
222  }
223 
224  float eCt, eCi, eCd;
225  Environment::GetSingleton()->GetFloatValue("BSP.decisionCost", eCd);
226  Environment::GetSingleton()->GetFloatValue("BSP.intersectionCost", eCi);
227  Environment::GetSingleton()->GetFloatValue("BSP.traversalCost", eCt);
228
229  // Ci should no be changed from 1.0, only Cd and Ct should be changed
230  app << "#SetBSP.(Cd, Ci, Ct) \tDecision, Intersection, ";
231  app << "and Traversal costs\n";
232  app << " ( " << eCd << ", " << eCi << " ," << eCt << " )\n";
233
234  return;
235}
236
237// Init required parameters
238void
239CKTBSBuildUp::InitReqPref(SReqPrefParams *pars)
240{
241  // nothing is obligatory as default
242  pars->reqPosition = Limits::Infinity;
243  pars->reqAxis = CKTBAxes::EE_Leaf;
244  pars->useReqAxis = false;
245 
246  // the values for automatic termination criteria, since it was
247  // not subdivided
248  pars->ratioLast = 1000.0;
249  pars->ratioLastButOne = 1000.0;
250  pars->failedSubDivCount = 0;
251
252  Environment::GetSingleton()->GetIntValue("BSP.axisSelectionAlg",
253                                           _algorithmForAxisSelection);
254  if (_algorithmForAxisSelection != 0) {
255    // The axis as prefered one - testing all 3 axes. Start
256    // from the axis cutting the longest side of the box
257    pars->reqAxis = (CKTBAxes::Axes)(wBbox.Diagonal().DrivingAxis());
258  }
259
260  return;
261}
262
263// creates all the auxiliary structures for building up CKTB tree
264CKTBSBuildUp::SInputData*
265CKTBSBuildUp::Init(ObjectContainer *objlist, const AxisAlignedBox3 &box)
266{
267#ifdef  _RANDOMIZE_POSITION
268  //world->_environment.GetFloat("SetBSP.randomizePosition.Eps", randomEps);
269#endif // _RANDOMIZE_POSITION
270 
271#ifdef _KTB_CONSTR_STATS
272  _stats_interiorCount = 0;
273  _stats_bboxCount = 0;
274  _stats_leafNodeCount = 0;
275  _stats_emptyLeftNodeCount = 0;
276  // Aggregate statistics
277  _maxDepth = 0; 
278  _sumLeafDepth = 0;
279  _sumFullLeafDepth = 0;
280  _sumObjectRefCount = 0;
281  _maxObjectRefInLeaf = 0;
282  // surface areas
283  _sumSurfaceAreaLeaves = 0.f;
284  _sumSurfaceAreaMULcntLeaves = 0.f;
285  _sumSurfaceAreaInteriorNodes = 0.f;
286#endif
287
288  // If to print out the tree during contruction
289  Environment::GetSingleton()->GetBoolValue("BSP.printCuts", _printCuts);
290
291  // if to cut off empty space in postprocessing and how it is performed
292  Environment::GetSingleton()->GetIntValue("BSP.absMaxAllowedDepth",
293                                           absMaxAllowedDepth);
294  Environment::GetSingleton()->GetIntValue("BSP.maxEmptyCutDepth",
295                                           maxEmptyCutDepth);
296  Environment::GetSingleton()->GetIntValue("BSP.algAutoTermination",
297                                           algorithmAutoTermination);
298  Environment::GetSingleton()->GetFloatValue("BSP.biasFreeCuts",
299                                             biasFreeCuts);
300
301  // if to make tagging lists inside the tree
302  Environment::GetSingleton()->GetBoolValue("BSP.minBoxes.use",
303                                            makeMinBoxes);
304  Environment::GetSingleton()->GetBoolValue("BSP.minBoxes.tight",
305                                            makeTightMinBoxes);
306  Environment::GetSingleton()->GetIntValue("BSP.minBoxes.minObjects",
307                                            minObjectsToCreateMinBox);
308  Environment::GetSingleton()->GetIntValue("BSP.minBoxes.minDepthDistance",
309                                           minDepthDistanceBetweenMinBoxes);
310  Environment::GetSingleton()->GetFloatValue("BSP.minBoxes.minSA2ratio",
311                                             minSA2ratioMinBoxes);
312
313  // How many items can be allocated at once
314  int maxItemsAtOnce = 1;
315  if (makeMinBoxes) {
316#ifdef _KTB8Bytes   
317    // We need to allocate for boxes the memory in a row
318    maxItemsAtOnce = 5; // 8x5=40 = 16+24;
319#else
320    maxItemsAtOnce = 4; // 12x4=48 = 24+24;
321#endif // _KTB8Bytes
322    assert(minObjectsToCreateMinBox >= 1);
323    assert(minDepthDistanceBetweenMinBoxes >= 0);
324    assert(minDepthDistanceBetweenMinBoxes <= 50);
325  }
326
327  // the array of preferred parameters used as a stack
328  pars = new SReqPrefParams[MAX_HEIGHT];
329  assert(pars);
330 
331  // Initiate the allocator
332  InitAux(0, MAX_HEIGHT - 1, maxItemsAtOnce);
333 
334  // since six planes are enough to cull empty space from leaves
335  if ( (maxEmptyCutDepth < 0) ||
336       (maxEmptyCutDepth > 6) ) {
337    cerr << "BSP.maxEmptyCutDepth = " << maxEmptyCutDepth
338          << " must be in range <0,6>\n";
339    abort();;
340  }
341
342  wBbox.Convert(box); // the box of the whole scene
343  boxSize = box.Diagonal(); // the size of the box along the axes
344 
345  wholeBoxArea = wBbox.SA2(); // the half of the surface area of the box
346  // the array of bounding boxes and duplications to the objects
347  int objlistcnt = objlist->size();
348
349  SInputData* data = AllocNewData();
350  // set the box
351  data->box = wBbox;
352  data->tightbox = wBbox;
353  data->objlist = new ObjectContainer(objlist->begin(), objlist->end());
354  data->count = objlistcnt;
355
356  // Initialization of the first value to insert the min boxes
357  data->lastMinBoxSA2 = box.SurfaceArea() * 10000.f;
358  data->lastDepthForMinBoxes = -20; // for root you always put the node
359  data->lastMinBoxNode = 0;
360 
361  // Init the parameters from the environment file
362  InitReqPref(&(data->pars));
363
364  // the number of buckets;
365  state.bucketN = 128;
366  state.bucketMin = new int[state.bucketN];
367  assert(state.bucketMin);
368  state.bucketMax = new int[state.bucketN];
369  assert(state.bucketMin);
370 
371  // Init the termination criteria
372  string termCrit;
373  Environment::GetSingleton()->GetStringValue("BSP.termCrit", termCrit);
374
375  maxCountTrials = 0;
376  if (!termCrit.compare("adhoc")) {
377    // termination criteria are the max depth of the hierarchy, number
378    // of primitives
379    data->modeSubDiv = EE_SUBDIVADHOC;
380    Environment::GetSingleton()->GetIntValue("BSP.maxListLength",
381                                             maxListLength);
382    Environment::GetSingleton()->GetIntValue("BSP.maxDepthAllowed",
383                                             maxDepthAllowed);
384  }
385  else {
386    if ( (!termCrit.compare("auto")) ||
387         (!termCrit.compare("auto2")) ) {
388      // automatic termination criteria are used, everything is computed
389      // automatically from number of objects etc.
390      int estHeight = (int)(log((float)initcnt)/log((float)2.0) + 0.9);
391      // cout << "EstHeight=" << estHeight << endl;
392      maxDepthAllowed = (int)((float)estHeight * 1.15f + 1.6f);
393
394      // maximum number of trials to further subdivide
395      maxCountTrials = (int)(1.0 + 0.2 * (float)maxDepthAllowed);
396      if (maxCountTrials < 3)
397        maxCountTrials = 3;
398     
399      // for 'auto2' we specify the length of the object list by hand
400      if (!termCrit.compare("auto2"))
401        Environment::GetSingleton()->GetIntValue("BSP.maxListLength",
402                                                 maxListLength);
403      else
404        // for 'auto' we set the length of the object list that is supposed
405        // to be the best for ray shooting
406        maxListLength = 1;
407
408      if (verbose) {
409        cout << "TERMCRIT:maximum height of a leaf set to "
410               << maxDepthAllowed << "\n";
411        cout << "TERMCRIT:maximum number of objects in leaf set to "
412               << maxListLength << "\n" ;
413      }
414      // set the mode to be recognized
415      data->modeSubDiv = EE_SUBDIVAUTOMATIC;
416    }
417    else {
418      cerr << " unknown termination criteria for BSP tree\n";
419      cerr << "It was specified: " << termCrit << "\n";
420      cerr << "Allowed types = auto, auto2, adhoc" << endl;
421      exit(3);;
422    }
423  }
424
425  // for some evaluation it is necessary to determine the following costs
426  Environment::GetSingleton()->GetFloatValue("BSP.traversalCost", Ct);
427  Environment::GetSingleton()->GetFloatValue("BSP.intersectionCost", Ci);
428  if (verbose)
429    cout << "Ct=" << Ct << " Ci=" << Ci << "\n";
430 
431  if (cutEmptySpace) {
432    // correct the maximum abs depth, when late cutting is allowed
433    if (absMaxAllowedDepth < maxDepthAllowed)
434      absMaxAllowedDepth = maxDepthAllowed;
435    if (absMaxAllowedDepth > (maxDepthAllowed + maxEmptyCutDepth) )
436      absMaxAllowedDepth = maxDepthAllowed + maxEmptyCutDepth;
437    if (absMaxAllowedDepth >= MAX_HEIGHT) {
438      absMaxAllowedDepth = MAX_HEIGHT;
439      maxDepthAllowed = MAX_HEIGHT - maxEmptyCutDepth;
440    }
441   
442    if (verbose) {
443      cout << "Cutting off empty spaces ON\n";
444      cout << "BSP.absMaxAllowedDepth = " << absMaxAllowedDepth << "\n";
445      cout << "BSP.maxEmptyCutDepth = " << maxEmptyCutDepth << endl;
446    }
447  }
448
449  if (verbose)
450    cout << "MaxListLength = " << maxListLength << endl; 
451
452  return data;
453}
454
455// interface function for building up CKTB tree
456SKTBNodeT*
457CKTBSBuildUp::BuildUp(ObjectContainer &objlist,
458                      const AxisAlignedBox3 &box,
459                      bool verboseF)
460{
461  // check the number of objects
462  if (objlist.size() == 0)
463    return 0; // nothing
464
465  // ---------------------------------------------------
466  // Rewriting the triangles to the just wrappers to save
467  // the memory during building!!!!
468  bool useWrappers = false;
469  if (useWrappers) {
470    cout << "WARNING: using only wrappers, not objects to save the memory!"
471         << endl;
472    cout << "Size of(AxisAlignedBox3Intersectable) = " << sizeof(AxisAlignedBox3Intersectable) << endl;
473    cout << "Size of(TriangleIntersectable) = " << sizeof(TriangleIntersectable) << endl;
474    for (ObjectContainer::iterator it = objlist.begin();
475         it != objlist.end(); it++) {
476      // take the triangle
477      Intersectable *i = *it;
478      // store the properties to new variables
479      AxisAlignedBox3 b = i->GetBox();
480      int mId = i->mId;
481      // delete the triangle
482      delete i;
483      // create the wrapper with the same box as triangle
484      AxisAlignedBox3Intersectable *a = new AxisAlignedBox3Intersectable(b);
485      a->mId = mId;
486      // and put it back to the list of objects
487      *it = a;
488    } // for   
489    cout << "Rewriting to wrappers is finished!" << endl;
490  } // if
491  // ---------------------------------------------------
492
493  initcnt = objlist.size();
494 
495  // copy verbose to be used consistenly further on
496  this->verbose = verboseF;
497 
498  // the boundary entries
499  SInputData *d = Init(&objlist, box);
500
501 
502  if (verbose)
503    cout << "Building up KTB tree for " << objlist.size()
504                 << " objects by sampling." << endl << flush;
505 
506  // -----------------------------------------------
507  // Start to build the tree
508  if ( (cutEmptySpace) &&
509       (initcnt <= maxListLength)) {
510    // only cutting off empty space for initial leaf
511    d->modeSubDiv = EE_SUBDIVCUTEMPTY;
512    root = SubDiv(d);
513  }
514  else {
515    // normal subdivision scheme, creating whole tree
516    root = SubDiv(d);
517  }
518
519#ifdef _DEBUG
520  if (GetDepth() != 0) {
521    cerr << "Using depth value does not work correctly, depth = "
522          << GetDepth() << endl;
523  }
524#endif 
525  assert(root); 
526
527  // the pointer to the root node
528  return GetRootNode();
529}
530
531int
532CKTBSBuildUp::UpdateEvaluation(float &eval, const float &newEval)
533{
534  if (newEval < eval) {
535    eval = newEval;
536    return 1; // updated
537  }
538  return 0; // not updated
539}
540
541// recursive function for subdividing the CKTB tree into two
542// halves .. creates one node of CKTB tree
543// list .. list of boundaries, count .. number of objects in current node
544// bb .. current bounding box of node, (pars,reqGlobAxis, modeSubDiv) .. cut pref.
545SKTBNodeT*
546CKTBSBuildUp::SubDiv(SInputData *d)
547{
548#if 0
549  static int index = 0;
550  cout << "SubDiv index = " << index << endl;
551  const int indexToFind = 13916;
552  if (index == indexToFind) {
553    cout << " IndexToFind = " << indexToFind << endl;
554    cout << " You should start debug" << endl;
555  }
556  index++;
557#endif
558
559  // This is the node to return from this function
560  SKTBNodeT *nodeToReturn = 0;
561
562  // if to make the interior node extended by the box here
563  makeMinBoxHere = false;
564
565  // -----------------------------------------------------
566  if (makeMinBoxes) {
567    if ( (d->count >= minObjectsToCreateMinBox) &&
568         (GetDepth() - d->lastDepthForMinBoxes >= minDepthDistanceBetweenMinBoxes) )  {
569      makeMinBoxHere = true;
570      d->lastDepthForMinBoxes = GetDepth();
571      d->lastMinBoxSA2 = d->box.SA2();
572    }
573  }
574
575  if ( (makeMinBoxHere) && (makeTightMinBoxes) ) {
576    // In case we should make the box of the current
577    // interior node and this should be tight, we have
578    // to first compute its extent before splitting
579    SBBox tightbox;
580    GetTightBox(*d, tightbox);
581    // Here we decrease the box to the minimum size
582    d->box = tightbox;
583  }
584 
585  // ------------------------------------------------------------------
586  // if we are forced to make subdivision at given point
587  if (d->pars.reqPosition != Limits::Infinity) {
588    // this code preceeds switch(mode) since 2-level evaluation
589    d->position = d->pars.reqPosition;
590    d->axis = d->pars.reqAxis;
591    // next subdivison will not be forced
592    d->pars.reqPosition = Limits::Infinity;
593    d->pars.reqAxis = CKTBAxes::EE_Leaf;
594    IncDepth();
595    SKTBNodeT* n = MakeOneCut(d);
596    if (!nodeToReturn)
597      nodeToReturn = n;
598    DecDepth();
599    return nodeToReturn;
600  }
601 
602  // -------------------------------------------------------------------------
603  // determine between subdivision modes - different termination criteria
604  switch (d->modeSubDiv) { // subdivision mode
605    case EE_SUBDIVCUTEMPTY: { // cutting off empty space
606      if ((GetDepth() >= absMaxAllowedDepth) ||
607          (GetDepth() >= (startEmptyCutDepth + maxEmptyCutDepth)) ||
608          (d->count == 0) ) {
609        SKTBNodeT* n = MakeLeaf(d);
610        if (!nodeToReturn)
611          nodeToReturn = n;
612        // cout << "LEC\n";
613        return nodeToReturn;
614      }
615      // to require the axis has no meaning in cutting off
616      d->pars.reqAxis = d->axis = CKTBAxes::EE_Leaf;
617      break;
618    }
619    // both adhoc, clustering and automatic termination criteria
620    case EE_SUBDIVADHOC:
621    case EE_SUBDIVAUTOMATIC: {
622      // automatic termination criteria are used in this phase as
623      // the absolute threshold similarly to SUBDIVADHOC mode
624      if ( (GetDepth() >= maxDepthAllowed) ||
625           (d->count <= maxListLength)) {
626        SKTBNodeT* n = MakeLeaf(d);
627        if (!nodeToReturn)
628          nodeToReturn = n;
629        // cout << "LLA\n";
630        return nodeToReturn;
631      }
632      break;
633    }
634    default: {
635      cerr << " unknown subdivision mode in ibsp.cpp\n";
636      abort();;
637    }
638  } // switch(modeSubDiv)
639 
640  // the axis where splitting should be done
641  CKTBAxes::Axes axis = CKTBAxes::EE_Leaf;
642  d->twoSplits = -1;
643  d->bestCost = WorstEvaluation() * 0.5f;
644  d->position = MAXFLOAT;
645  const bool secondPositionMax = true;
646
647#define CallGetSplitPlane GetSplitPlaneOpt
648 
649  // if the axis is required in pars structure for some reasons
650  if (d->pars.useReqAxis) {
651    axis = d->pars.reqAxis;
652    // The next subdivision is not prescribed
653    d->pars.useReqAxis = false;
654    switch(axis) {
655      case CKTBAxes::EE_X_axis : {
656        state.InitXaxis(d->count, d->box); // init
657        CallGetSplitPlane(d, 0); // evaluate
658        break;
659      }
660      case CKTBAxes::EE_Y_axis : {
661        state.InitYaxis(d->count, d->box); // init
662        CallGetSplitPlane(d, 1); // evaluate
663        break;
664      }
665      case CKTBAxes::EE_Z_axis : {
666        state.InitZaxis(d->count, d->box); // init
667        CallGetSplitPlane(d, 2); // evaluate
668        break;
669      }
670      default: {
671        cerr << "No other option is allowed here" << endl;
672        abort();;
673      }
674    }
675    // compute the cost, normalized by surface area
676    state.NormalizeCostBySA2();
677    d->bestCost /= state.areaWholeSA2;
678    if (UpdateEvaluation(d->bestCost, state.bestCost)) {
679      // Compute correct position to be used for splitting
680      d->position = state.bestPosition;
681      d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
682      if (d->twoSplits > 1) {
683        if (secondPositionMax)
684          d->position2 = state.bestPosition2;
685        else
686          d->position2 = state.box.Max(axis);
687      }
688    }
689  }
690  else {
691    int algorithmForAxisSel = _algorithmForAxisSelection;
692    if (d->modeSubDiv == EE_SUBDIVCUTEMPTY)
693      algorithmForAxisSel = 0;
694
695    // only a single axis will be tested
696    bool useSingleAxis = true;
697       
698    // the required axis is by global function .. e.g. cyclic change x, y, z
699    switch (algorithmForAxisSel) {
700      case 0: {
701        // all three axis will be used, starting from x, y, z
702        useSingleAxis = false;
703        break;
704      }
705      case 1: {
706        // cyclic order of axes, x, y, z, x, y....
707        axis = d->pars.reqAxis;
708        // compute the axis for the next subdivision step
709        d->pars.reqAxis = oaxes[axis][0];
710        break;
711      }
712      case 2 : {
713        // The next time will be determined the same way
714        axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
715        break;
716      }
717      case 3 : {
718        float p = RandomValue(0.0f, 1.0f);
719        const float thresholdAxisP = 0.8f;
720        if (p < thresholdAxisP) {         
721          axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
722        }
723        else {
724          // Use the axis prescribed from previous step
725          axis = d->pars.reqAxis;
726        }
727        // for the next time, different axis
728        d->pars.reqAxis = oaxes[axis][0];
729        break;
730      }
731      case 4 : {
732        float p = RandomValue(0.0f, 1.0f);
733        const float thresholdAxisP = 0.3f;
734        if (p < thresholdAxisP) {         
735          axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
736        }
737        else {
738          // Use the axis prescribed from previous step
739          axis = d->pars.reqAxis;
740        }
741        // for the next time, different axis
742        d->pars.reqAxis = oaxes[axis][0];
743        break;
744      }
745      case 5 : {
746        float p = RandomValue(0.0f, 1.0f);
747        const float thresholdAxisP = 0.2f;
748        d->pars.reqAxis = (CKTBAxes::Axes)-1;
749        if (p < thresholdAxisP) {
750          axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
751          // for the next time, different axis
752        }
753        else {
754          // Use the axis for which cost model gets minimum,
755          // compute it for all axes x, y, z
756          useSingleAxis = false;
757        }
758        break;
759      }
760      default: {
761        cerr << "Selection algorithm = " << algorithmForAxisSel
762              << " for axis is not implemented" << endl;
763        abort();;
764      }
765    } // switch
766
767    // How the algorithm is used
768    if (useSingleAxis) {
769      // -------------------------------------------
770      // Compute subdivision only for one axis
771      switch(axis) {
772        case CKTBAxes::EE_X_axis : {
773          state.InitXaxis(d->count, d->box); // init
774          CallGetSplitPlane(d, 0); // evaluate
775          break;
776        }
777        case CKTBAxes::EE_Y_axis : {
778          state.InitYaxis(d->count, d->box); // init
779          CallGetSplitPlane(d, 1); // evaluate
780          break;
781        }
782        case CKTBAxes::EE_Z_axis : {
783          state.InitZaxis(d->count, d->box); // init
784          CallGetSplitPlane(d, 2); // evaluate
785          break;
786        }
787        //case CKTBAxes::EE_Leaf : {
788        //goto MAKE_LEAF;
789        //}
790        default: {
791          cerr << "Selection algorithm = " << algorithmForAxisSel
792                << " for axis = " << axis <<" is not implemented" << endl;
793          abort();;
794        }
795      }
796      state.NormalizeCostBySA2();
797      d->bestCost /= state.areaWholeSA2;
798      if (UpdateEvaluation(d->bestCost, state.bestCost)) {
799        // Compute correct position to be used for splitting
800        d->position = state.bestPosition;
801        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
802        if (d->twoSplits > 1) {
803          if (secondPositionMax)
804            d->position2 = state.bestPosition2;
805          else
806            d->position2 = state.box.Max(axis);
807        }
808      }
809    }
810    else {
811      // no axis is required, we can pickup any we like. We test all three axes
812      // and we pickup the minimum cost for all three axes.
813      state.InitXaxis(d->count, d->box); // init for x axis
814      CallGetSplitPlane(d, 0); // evaluate for x axis
815      if (UpdateEvaluation(d->bestCost, state.bestCost)) {     
816        axis = CKTBAxes::EE_X_axis; // the x-axis should be used
817        // Compute correct position to be used for splitting
818        d->bestCost = state.bestCost;
819        d->position = state.bestPosition;
820        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
821        if (d->twoSplits > 1) {
822          if (secondPositionMax)
823            d->position2 = state.bestPosition2;
824          else
825            d->position2 = state.box.Max(axis);
826        }
827      }
828      // -----------------------
829      // now test for y axis
830      state.InitYaxis(d->count, d->box); // init for x axis
831      CallGetSplitPlane(d, 1); // evaluate for x axis
832      if (UpdateEvaluation(d->bestCost, state.bestCost)) {
833        axis = CKTBAxes::EE_Y_axis; // the y-axis should be used
834        // Compute correct position to be used for splitting
835        d->bestCost = state.bestCost;
836        d->position = state.bestPosition;
837        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
838        if (d->twoSplits > 1) {
839          if (secondPositionMax)
840            d->position2 = state.bestPosition2;
841          else
842            d->position2 = state.box.Max(axis);
843        }
844      }
845      // -----------------------
846      // now test for z axis
847      state.InitZaxis(d->count, d->box); // init for x axis
848      CallGetSplitPlane(d, 2); // evaluate for x axis
849      if (UpdateEvaluation(d->bestCost, state.bestCost)) {
850        axis = CKTBAxes::EE_Z_axis; // the z-axis should be used
851        // Compute correct position to be used for splitting
852        d->bestCost = state.bestCost;
853        d->position = state.bestPosition;
854        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
855        if (d->twoSplits > 1) {
856          if (secondPositionMax)
857            d->position2 = state.bestPosition2;
858          else
859            d->position2 = state.box.Max(axis);
860        }
861      }
862      // Now we have to renormalize the cost for purpose of termination the building
863      d->bestCost /= state.areaWholeSA2;
864    } // arbitrary order of axes
865  } // for which axis to compute surface area heuristics
866
867  // ----------------------------------------------------------
868  // when automatic termination criteria should be used
869  if (d->modeSubDiv == EE_SUBDIVAUTOMATIC) {
870    bool stopSubdivision = false;
871    float ratio;
872    // which algorithm to use for automatic termination criteria
873    switch (algorithmAutoTermination) {
874      // --------------------------------------------------
875      case 0: { // This is described in my thesis.
876        // This is the cost of unsubdivided node
877        float leafCost = (float)(d->count);
878        ratio = d->bestCost / leafCost;
879        const float minRatio = 0.75;
880        if (ratio > minRatio)
881          stopSubdivision = true;
882        break;
883      }
884      // --------------------------------------------------
885      case 1: {
886        // This is the cost of unsubdivided node, with uniformly distributed
887        // objects and the spatial median subdivision, without any object straddling
888        // the splitting plane. In some sense ideal setting for uniform distribution.
889        Vector3 s = d->box.Diagonal();
890        // Taking the widest side of the object to be subdivided
891        int diagAxis = s.DrivingAxis();
892        // Half of the width
893        // ----------------------
894        float w2 = s[diagAxis] * 0.5;
895        float t = s[oaxes[diagAxis][0]];
896        float h = s[oaxes[diagAxis][1]];
897        float sah2 = w2*(t+h) + h*t; // half of the surface area of a child
898        // This is the cost for case when the node is subdivided in the half and
899        // the half of the object is on the left and half of the objects on the right
900        // This is not the best cost possible !!!
901        // = float subdividedCost = 2 * (sah2/state.areaWholeSA2 * count/2);
902        float subdividedCost = sah2/state.areaWholeSA2 * (float)(d->count);
903        ratio = d->bestCost / subdividedCost;
904        // This is the max ratio allowed for splitting node subdivision
905
906        // The computed ratio = 1.0 corresponds to 'ideal case', so the threshold
907        // must be higher than 1.0 !
908        const float minRatio = 1.1;
909        if (ratio > minRatio)
910          stopSubdivision = true;
911        break;
912      }
913      default : {
914        cerr << "Uknown algorithm for automatic termination criteria = "
915              << algorithmAutoTermination << endl;
916      }
917       
918    } // switch
919   
920    //cout << "R=" << ratio << " ";
921    if (stopSubdivision) {
922      // cout << "F" << pars.failedSubDivCount << " ";
923      // when small improvement by the subdivision is reached
924      d->pars.failedSubDivCount++;
925      if (d->pars.failedSubDivCount > maxCountTrials) {
926        // make leaf and finish with this KD-tree branch
927        // cout << "L, cnt=" << count << " ,d=" << GetDepth()()
928        //     << " ,c=" << bestQ << " ,rat=" << ratio << "\n";
929        // possibly cut empty space before leaf is created
930        if (cutEmptySpace) {
931          // setting initial depth if empty space cutting
932          startEmptyCutDepth = GetDepth();
933          IncDepth();
934          d->modeSubDiv = EE_SUBDIVCUTEMPTY;
935          // create the leaf first, with late empty cutting
936          SKTBNodeT *n = SubDiv(d);
937          if (!nodeToReturn)
938            nodeToReturn = n;
939          DecDepth();
940          return nodeToReturn; // and finish by constructing a leaf
941        }
942
943        // make leaf and finish
944        SKTBNodeT *n = MakeLeaf(d);
945        if (!nodeToReturn)
946          nodeToReturn = n;
947        return nodeToReturn;
948      }
949    }
950
951    // Do not stop subdivision now
952   
953    // update the progress in the parameters
954    d->pars.ratioLastButOne = d->pars.ratioLast;
955    d->pars.ratioLast = ratio;
956
957    // assure that axis is defined, if finding the splitting plane has failed
958    if ( (axis == CKTBAxes::EE_Leaf) ||
959         (d->position == MAXFLOAT) ) {
960      // compute the spatial median for arbitrary axis
961      axis = (CKTBAxes::Axes)int(RandomValue(0.f, 2.98f));
962      // set position based algorithm for the next cut and children
963      d->position = (d->box.Min(axis) + d->box.Max(axis)) * 0.5f;
964      d->twoSplits = 1;
965    }
966  } // subdivision automatic
967
968  // -----------------------------------------------
969  // We have evaluated surface area heuristics above
970  // if no winner for subdivision exists, make a leaf
971  if ( (axis == CKTBAxes::EE_Leaf) ||
972       (d->position == MAXFLOAT) )
973  {
974    //cerr << "This shall not happen" << endl;
975    SKTBNodeT *n = MakeLeaf(d);
976    if (!nodeToReturn)
977      nodeToReturn = n;
978    // cout << "LL\n";
979    return nodeToReturn;
980  }
981 
982  if (verbose) {
983#ifdef _DEBUG
984
985//#define _KTBPRINTCUTS
986#ifndef _KTBPRINTCUTS
987    if (GetDepth() == 0)
988#else
989    if (_printCuts)
990#endif
991    {
992      cout << "position: " << d->position << ", axis: "
993             << (int)axis << ", depth=" << GetDepth() << ", box:" << endl;
994      Describe(d->box, cout, 0);
995      cout << endl;
996    }
997#endif // _DEBUG
998  }
999
1000#ifdef _DEBUG
1001  if (d->twoSplits == -1) {
1002    cerr << "Some problem in implementation" << endl;
1003    abort();;
1004  }
1005#endif
1006
1007  // copy the parameters to use for splitting
1008  d->axis = axis;
1009
1010  if (d->twoSplits == 1) {
1011    // create interior node with two successors
1012    // case (1)
1013    // DEBUG << "case 1 \n";
1014    // Make easy cut, using one position
1015    IncDepth();
1016    SKTBNodeT* n = MakeOneCut(d);
1017    if (!nodeToReturn)
1018      nodeToReturn = n;
1019    DecDepth();
1020    return nodeToReturn;
1021  }
1022
1023#if 1
1024  cerr << "Unexpected use in this implementation" << endl;
1025  cerr << "ABORTING " << endl;
1026  abort();;
1027#endif
1028 
1029//  case       case
1030//  (2)   or   (3) structure must be created
1031//   O          O                #
1032//  / \        / \               #
1033// L  RI      LI  R              #
1034//   /  \    / \                 #
1035//   RL RR  LL LR                #
1036
1037#ifdef _DEBUG
1038  if ( (d->box.Min()[axis] >= d->position) ||
1039       (d->position >= d->position2) ||
1040       (d->position2 >= d->box.Max()[axis]) ) {
1041    cerr << " the case 2 or 3 is defined incorrectly\n";
1042    abort();;
1043  }
1044#endif
1045
1046  // To be reworked !!!
1047  abort();;
1048  if (d->twoSplits == 2) {
1049    // -----------------------------------------
1050    // case (2)
1051    // DEBUG << "case 2 \n";
1052
1053    // make L part, linked in DFS order
1054    int cntLeft = 1;
1055    int cntRight = 0;
1056    IncDepth();
1057    SKTBNodeT *n = MakeOneCut(d);
1058    if (!nodeToReturn)
1059      nodeToReturn = n;
1060
1061    // make RI part, linked to the left node explictly
1062    SBBox rbb = d->box; // the bounding box
1063    rbb.Reduce(axis, 1, d->position);
1064    d->box = rbb;
1065    d->pars.reqAxis = axis;
1066    d->pars.reqPosition = d->position2;   
1067    SKTBNodeT *n2 = SubDiv(d);
1068    DecDepth();
1069    //Link the node
1070    SetInteriorNodeLinks(n, 0, n2);
1071    return nodeToReturn;
1072  }
1073
1074  // -----------------------------------------
1075  // case (3) ????????????????? I am not sure if this is correct implementation !!! VH
1076  // DEBUG << "case 3 \n";
1077 
1078  // make LI part
1079  SBBox lbb = d->box; // the bounding box
1080  lbb.Reduce(axis, 0, d->position2);
1081  int cntLeft = 0;
1082  int cntRight = 1;
1083  d->box = lbb;
1084  d->position = d->position2;
1085  d->axis = axis;
1086  // the next subdivision step for R part
1087  d->pars.reqAxis = axis;
1088  d->pars.reqPosition = d->position;
1089 
1090  IncDepth();
1091  SKTBNodeT *nl = SubDiv(d);
1092  if (!nodeToReturn)
1093    nodeToReturn = nl;
1094 
1095  // make R part
1096  SKTBNodeT *n = MakeOneCut(d);
1097  DecDepth();
1098  SetInteriorNodeLinks(nl, 0, n);
1099  // return left node, linked in DFS order to the previous one
1100  return nodeToReturn;
1101}
1102
1103
1104
1105// makes cutting on the given position on given axis
1106// at value position, bb is the bounding box of current node,
1107// count is the number of objects in the node
1108// flags LeftF and RightF declares if to continue down after the recursion
1109SKTBNodeT *
1110CKTBSBuildUp::MakeOneCut(SInputData *d)
1111{
1112#ifdef _DEBUG
1113  assert(d->position >= d->box.Min(d->axis));
1114  assert(d->position <= d->box.Max(d->axis));
1115#endif
1116 
1117  SKTBNodeT *nodeToReturn = 0;
1118
1119#ifdef  _KTB_CONSTR_STATS
1120  // surface area of the interior nodes to be subdivided
1121  _sumSurfaceAreaInteriorNodes += d->box.SA2();
1122#endif // _KTB_CONSTR_STATS
1123 
1124  //if (splitClip) {
1125  //  // empty object list to store the pointers to the split objects
1126  //  ObjectContainer objlist;
1127  //  cerr << "Not yet handled" << endl;
1128  //  abort();;
1129  // }
1130
1131#if 0
1132  // check if the lists are correctly sorted .. in debug
1133  static int index = 0;
1134#if 1
1135  cout << "MakeOneCut index = " << index << " .. check axis = "
1136        << d->axis << " count = " << d->count
1137        << " depth = " << GetDepth() << endl;
1138  const int indexToFind = 3743;
1139  if (index == indexToFind) {
1140   
1141    cout << "Debug should start " << endl;
1142    cout << "index = indexToFind" << endl;
1143  }
1144#endif
1145  index++; 
1146#endif
1147   
1148  // We have to allocate a new node for right child
1149  SInputData* rightData = AllocNewData();
1150  // Copy the basic data from parent
1151  rightData->CopyBasicData(d);
1152  // We compute the tight box after the splitting
1153  rightData->tightbox.Initialize();
1154  rightData->objlist = new ObjectContainer;
1155  assert(rightData->objlist);
1156
1157  // we compute the tight box for the objects on the
1158  // left of the splitting plane
1159  d->tightbox.Initialize();
1160 
1161  // Simply go through the list of primitives and those that belong
1162  // to the left keep on the left (d), those that belong to the right
1163  // assign to the right list
1164  int cntL = 0, cntR = 0;
1165  {
1166    ObjectContainer::iterator endit = d->objlist->end();
1167    int axis = d->axis;
1168    AxisAlignedBox3 obox;
1169    for (ObjectContainer::iterator it = d->objlist->begin(); it != endit; it++) {
1170      Intersectable *obj = *it;
1171      obox = obj->GetBox();
1172      if (obox.Max(axis) > d->position) {
1173        rightData->objlist->push_back(obj);
1174        cntR++;
1175        // the number of objects on the right of the splitting plane
1176        rightData->tightbox.Include(obox);
1177      }
1178      if (obox.Min(axis) >= d->position) {
1179        // we simply shorten the list - this object belongs only to the right list
1180        (*it) = *(endit-1);
1181        it--;
1182        endit--;
1183      }
1184      else {
1185        // the number of objects on the right of the splitting plane
1186        cntL++;
1187        d->tightbox.Include(obox);
1188      }
1189    } // for
1190  }
1191
1192  assert(cntL + cntR >= d->count);
1193
1194  // ---------------------------------------------
1195  // initialize the left bounding box 'bb'
1196  d->box.Reduce(d->axis, 0, d->position);
1197  d->objlist->resize(cntL);
1198 
1199  // initialize the right bounding box 'rbb'
1200  rightData->box.Reduce(d->axis, 1, d->position);
1201
1202  // Set correct number of objects
1203  d->count = cntL;
1204  rightData->count = cntR;
1205
1206  d->tightbox = Intersect(d->tightbox, d->box);
1207  rightData->tightbox =
1208    Intersect(rightData->tightbox, rightData->box);
1209 
1210  // Allocate the representation of the interior node
1211  SKTBNodeT *node = 0;
1212
1213  if (makeMinBoxHere) {
1214    // -------------------------------------------------------
1215    // interior node with the box
1216    node = AllocInteriorNodeWithBox(// interior node data
1217                                    d->axis, d->position, cntL, cntR,
1218                                    // box data
1219                                    d->box, d->lastMinBoxNode,
1220                                    d->lastDepthForMinBoxes);
1221
1222    //cout << "depth = " << GetDepth() << " node = " << node
1223    //   << " lastMinBoxNode = " << d->lastMinBoxNode << endl;
1224   
1225    d->lastMinBoxNode = node; // we have to remember the last node
1226    rightData->lastMinBoxNode = node;
1227    makeMinBoxHere = false; // reset for the next time
1228  }
1229  else
1230    // normal interior node (=without box)
1231    node = AllocInteriorNode(d->axis, d->position, cntL, cntR);
1232 
1233#ifdef _DEBUG
1234  if (!node) {
1235    cerr << "Allocation of Interior Node failed.\n";
1236    abort();;
1237  }
1238#endif // _DEBUG
1239  // set the node to be returned
1240  if (!nodeToReturn)
1241    nodeToReturn = nodeToLink;
1242  assert(node);
1243  assert(nodeToReturn);
1244 
1245  // ---------------------------------------------
1246  // initialize the left bounding box 'bb'
1247  d->box.Reduce(d->axis, 0, d->position);
1248
1249  // initialize the right bounding box 'rbb'
1250  rightData->box.Reduce(d->axis, 1, d->position);
1251
1252  // Set correct number of objects
1253  d->count = cntL;
1254  rightData->count = cntR;
1255 
1256  // The children nodes to be created
1257  SKTBNodeT *nodeL=0, *nodeR=0;
1258
1259  // Recurse to the left child
1260  if (d->makeSubdivisionLeft) {
1261    // subdivide branches
1262    ESubdivMode modeLeft = d->modeSubDiv ;
1263    if ( (cutEmptySpace) &&
1264         (modeLeft != EE_SUBDIVCUTEMPTY) &&
1265         (d->count <= maxListLength) ) {
1266      // setting initial depth if empty space cutting
1267      startEmptyCutDepth = GetDepth();
1268      d->modeSubDiv = EE_SUBDIVCUTEMPTY;
1269    }
1270    // DEBUG << "Left cnt= " << (cnt[0] >> 1) << " depth= "
1271    //       << GetDepth() << "\n";
1272    nodeL = SubDiv(d);
1273  }
1274 
1275  // Recurse to the right child
1276  if (rightData->makeSubdivisionRight) {
1277    ESubdivMode modeRight = d->modeSubDiv ;
1278    if ( (cutEmptySpace) &&
1279         (modeRight != EE_SUBDIVCUTEMPTY) &&
1280         (rightData->count <= maxListLength) ) {
1281      // setting initial depth if empty space cutting
1282      startEmptyCutDepth = GetDepth();
1283      rightData->modeSubDiv = EE_SUBDIVCUTEMPTY;
1284    }
1285    // DEBUG << "Right cnt= " << (cnt[1] >> 1) << " depth= "
1286    // << GetDepth() << "\n";
1287    nodeR = SubDiv(rightData);
1288
1289  }
1290
1291  delete rightData->objlist;
1292  rightData->objlist = 0;
1293  delete d->objlist;
1294  d->objlist = 0;
1295 
1296  // Free the data to be reused further on
1297  FreeLastData();
1298
1299  // Set the node pointers to the children for created node
1300  SetInteriorNodeLinks(node, nodeL, nodeR);
1301
1302  return nodeToReturn; // return the node
1303}
1304
1305// returns a box enclosing all the objects in the node
1306void
1307CKTBSBuildUp::GetTightBox(const SInputData &d, SBBox &tbox)
1308{
1309  tbox = d.tightbox;
1310}
1311
1312// make the full leaf from current node
1313SKTBNodeT*
1314CKTBSBuildUp::MakeLeaf(SInputData *d)
1315{
1316#ifdef  _KTBPRINTCUTS
1317  cout << " Creating leaf, depth = " << GetDepth()
1318       << " cntObjs = " << count << endl;
1319  //exit(3);
1320#endif
1321 
1322#ifdef _KTB_CONSTR_STATS
1323  // update the statistics
1324  float sa2bb;
1325  _sumSurfaceAreaLeaves += (sa2bb = d->box.SA2());
1326  _sumSurfaceAreaMULcntLeaves += ((float)d->count * sa2bb);
1327#endif // _KTB_CONSTR_STATS
1328 
1329  // ------------------------------------------------------
1330  // The version with allocating empty leaves. This makes
1331  // sense since in DFS order it works correctly
1332  SKTBNodeT *node = AllocLeaf(d->count);
1333
1334  if (d->count > 0) {
1335    // copy the list of objects
1336    ObjectContainer *objlist = d->objlist;
1337    assert(objlist->size() == d->count);
1338    // Do not delete the object list, it is used as allocated above
1339    SetFullLeaf(node, objlist);
1340    // We delete the object list since this is used at the level of leaves
1341    d->objlist = 0;
1342  }
1343
1344  // We assume that the allocation of the single leaf cannot
1345  // fail, since this one shall be the last in the sequence.
1346  assert(node == nodeToLink);
1347 
1348  // Return the node to be used in the linking
1349  return nodeToLink;
1350}
1351
1352// --------------------------------------------------------------------
1353// class CKTBSBuildUp::CTestAx
1354
1355// The initialization for the first axis to be tested, this time *****Z axis*******
1356void
1357CKTBSBuildUp::SSplitState::InitXaxis(int cnt, const SBBox &boxN)
1358{
1359  // initialize the variables, mainly for surface area estimates
1360  cntAll = cnt;     // the number of all objects in the bounding box
1361  cntLeft = 0;      // the count of bounding boxes on the left
1362  cntRight = cnt;   // the count of bounding boxes on the right
1363  thickness = 0;    // the count of bounding boxes straddling the splitting plane
1364  axis = CKTBAxes::EE_X_axis; // the axis, where the splitting is proposed
1365  box = boxN;        // the box, that is subdivided
1366  //int topAxis = 1;        // axis that is considered in top .. depth
1367  //int frontAxis = 2;      // axis that is considered in front .. height 
1368  // the size of bounding box along the axis to be split
1369  width = box.Max().x - box.Min().x;
1370  // and along two other axes 
1371  topw = box.Max().y - box.Min().y;
1372  frontw = box.Max().z - box.Min().z;
1373  // surface area of the splitting plane .. one face !!
1374  areaSplitPlane = topw * frontw;
1375  // the sum of length of the boxes not to be subdivided
1376  areaSumLength = topw + frontw;
1377  // The surface are of the whole box
1378  areaWholeSA2 = width * areaSumLength + areaSplitPlane;
1379  // initial evaluation .. the worst cost
1380  bestCost = WorstEvaluation();
1381  // the buckets are zeroed
1382  InitBuckets();
1383}
1384
1385// The initialization for the first axis to be tested, this time *****Y axis*******
1386// This can be run independently.
1387void
1388CKTBSBuildUp::SSplitState::InitYaxis(int cnt, const SBBox &boxN)
1389{
1390  // initialize the variables, mainly for surface area estimates
1391  cntAll = cnt;     // the number of all objects in the bounding box
1392  cntLeft = 0;      // the count of bounding boxes on the left
1393  cntRight = cnt;   // the count of bounding boxes on the right
1394  thickness = 0;    // the count of bounding boxes straddling the splitting plane
1395  axis = CKTBAxes::EE_Y_axis; // the axis, where the splitting is proposed
1396  box = boxN;        // the box, that is subdivided
1397  //int frontAxis = oaxes[axis][0]; // axis that is considered in front .. height - x-axis
1398  //int topAxis = oaxes[axis][1]; // axis that is considered in top .. depth - y-axis
1399  // the size along the second axis - height
1400  frontw = box.Max().x - box.Min().x;
1401  // the size of bounding box along the axis to be split - width
1402  width = box.Max().y - box.Min().y;
1403  // The size along third axis - depth
1404  topw = box.Max().z - box.Min().z;
1405  // surface area of the splitting plane .. one face !!
1406  areaSplitPlane = topw * frontw;
1407  // the sum of length of the boxes not to be subdivided
1408  areaSumLength = topw + frontw;
1409  areaWholeSA2 = width * areaSumLength + areaSplitPlane;
1410  // initial evaluation .. the worst cost
1411  bestCost = WorstEvaluation();
1412  // the buckets are zeroed
1413  InitBuckets();
1414}
1415
1416// The initialization for the first axis to be tested, this time *****Z axis*******
1417// This can be run independently.
1418void
1419CKTBSBuildUp::SSplitState::InitZaxis(int cnt, const SBBox &boxN)
1420{
1421  // initialize the variables, mainly for surface area estimates
1422  cntAll = cnt;     // the number of all objects in the bounding box
1423  cntLeft = 0;      // the count of bounding boxes on the left
1424  cntRight = cnt;   // the count of bounding boxes on the right
1425  thickness = 0;    // the count of bounding boxes straddling the splitting plane
1426  axis = CKTBAxes::EE_Z_axis; // the axis, where the splitting is proposed
1427  box = boxN;       // the box, that is subdivided
1428  //int frontAxis = oaxes[axis][1]; // axis that is considered in front .. height - x axis
1429  //int topAxis = oaxes[axis][0]; // axis that is considered in top .. depth - y axis
1430  // The size along third axis - depth
1431  topw = box.Max().x - box.Min().x;
1432  // the size along the second axis - height
1433  frontw = box.Max().y - box.Min().y;
1434  // the size of bounding box along the axis to be split - width
1435  width = box.Max().z - box.Min().z;
1436  // surface area of the splitting plane .. one face !!
1437  areaSplitPlane = topw * frontw;
1438  // the sum of length of the boxes not to be subdivided
1439  areaSumLength = topw + frontw;
1440  areaWholeSA2 = width * areaSumLength + areaSplitPlane;
1441  // initial evaluation .. the worst cost
1442  bestCost = WorstEvaluation();
1443  // the buckets are zeroed
1444  InitBuckets();
1445}
1446
1447// This is the computation of the cost using surface area heuristics
1448// using LINEAR ESTIMATE
1449void
1450CKTBSBuildUp::EvaluateCost(SSplitState &state)
1451{
1452  // the surface area of the bounding box for the left child
1453  assert(state.position > -1e-9);
1454  float areaLeft = state.position * state.areaSumLength + state.areaSplitPlane;
1455
1456  // the surface area of the right bounding box for the right child
1457  assert(state.width - state.position > -1e-9);
1458  float areaRight = (state.width - state.position) * state.areaSumLength +
1459                    state.areaSplitPlane;
1460
1461  // computation of the cost .. smaller is better
1462  float cost = areaLeft * state.cntLeft + areaRight * state.cntRight;
1463
1464#ifdef _DEBUG_COSTFUNCTION
1465  // Put there normalized position and cost also normalized somehow
1466  ReportCostStream(state.position / state.width,
1467                   cost / (state.areaWholeSA2 * state.cntAll));
1468#endif
1469 
1470  // This normalization is not in fact necessary here
1471  //cost //= state.areaWholeSA2;
1472  if (cost < state.bestCost) {
1473    state.bestCost = cost;
1474    // The iterator pointing to the best boundary
1475    state.bestPosition = state.position;
1476    // The number of objects whose boundaries are duplicated
1477    //state.bestThickness = state.thickness;
1478  }
1479  return;
1480}
1481
1482// This is the computation of the cost using surface area heuristics
1483void
1484CKTBSBuildUp::EvaluateCostFreeCut(SSplitState &state)
1485{
1486  // This is assumed to work for free cut (no object intersected)
1487  assert(state.thickness == 0);
1488  assert((state.cntLeft == 0) || (state.cntRight == 0));
1489 
1490  // the surface area of the bounding box for the left child
1491  assert(state.position > -1e-9);
1492  float areaLeft = state.position * state.areaSumLength + state.areaSplitPlane;
1493
1494  // the surface area of the right bounding box for the right child
1495  assert(state.width - state.position > -1e-9);
1496  float areaRight = (state.width - state.position) * state.areaSumLength +
1497                    state.areaSplitPlane;
1498
1499  // computation of the cost .. smaller is better
1500  float cost = biasFreeCuts *(areaLeft * state.cntLeft + areaRight * state.cntRight);
1501
1502#ifdef _DEBUG_COSTFUNCTION
1503  // Put there normalized position and cost also normalized somehow
1504  ReportCostStream(state.position / state.width,
1505                   cost / (state.areaWholeSA2 * state.cntAll));
1506#endif
1507 
1508  // This normalization is not in fact necessary here
1509  //cost //= state.areaWholeSA2;
1510  if (cost < state.bestCost) {
1511    state.bestCost = cost;
1512    // The iterator pointing to the best boundary
1513    state.bestPosition = state.position;
1514    // The number of objects whose boundaries are duplicated
1515    // state.bestThickness = 0;
1516  }
1517  return;
1518}
1519
1520// ---------------------------------------------------------------------
1521// For a given X-axis search for the best splitting plane position.
1522// This is implemented by sampling
1523void
1524CKTBSBuildUp::GetSplitPlaneOpt(SInputData *d, int axisToTest)
1525{
1526  // some necessary space is required to cut off
1527  const float MinPosition = d->tightbox.Min(axisToTest);
1528  const float MaxPosition = d->tightbox.Max(axisToTest);
1529  const float bMinPosition = d->box.Min(axisToTest);
1530  state.bestTwoSplits = 1;
1531
1532  // This makes no sense to compute, if the width of the box
1533  // is zero!
1534  if (bMinPosition == d->box.Max(axisToTest)) {
1535    state.bestPosition = bMinPosition;
1536    return;
1537  }
1538
1539  int bucketN = state.bucketN;
1540
1541  for (int i = 0; i < bucketN; i++) {
1542    state.bucketMin[i] = 0;
1543    state.bucketMax[i] = 0;
1544  }
1545 
1546  float mult = (float)bucketN / (MaxPosition - MinPosition);
1547
1548  ObjectContainer::iterator endit = d->objlist->end();
1549  int axis = axisToTest;
1550  AxisAlignedBox3 obox;
1551  for (ObjectContainer::iterator it = d->objlist->begin(); it != endit; it++) {
1552    Intersectable *obj = *it;
1553    obox = obj->GetBox();
1554    float v = obox.Min(axis);
1555    if (v < MinPosition)
1556      state.bucketMin[0]++;
1557    else
1558      if (v > MaxPosition)
1559        state.bucketMin[bucketN-1]++;
1560      else {   
1561        assert(v >= MinPosition);
1562        assert(v <= MaxPosition);   
1563        int bucketPos = (int)((v - MinPosition) * mult);
1564        assert(bucketPos >= 0);
1565        assert(bucketPos <= bucketN);
1566        if (bucketPos == bucketN)
1567          bucketPos = bucketN-1;
1568        state.bucketMin[bucketPos]++;
1569      }
1570
1571    float v2 = obox.Max(axis);
1572    if (v2 < MinPosition)
1573      state.bucketMax[0]++;
1574    else
1575      if (v2 > MaxPosition)
1576        state.bucketMax[bucketN-1]++;
1577      else {         
1578        assert(v2 >= MinPosition);
1579        assert(v2 <= MaxPosition);   
1580        int bucketPos2 = (int)((v2 - MinPosition) * mult);
1581        assert(bucketPos2 >= 0);
1582        assert(bucketPos2 <= bucketN);
1583        if (bucketPos2 == bucketN)
1584          bucketPos2 = bucketN-1;
1585        state.bucketMax[bucketPos2]++;
1586      }
1587  } // for
1588
1589  // First evaluate the cost, when the splitting plane is on the left
1590  state.position = MinPosition - bMinPosition;
1591  assert(state.position >= 0.f);
1592  assert(state.position <= state.width);
1593  state.cntLeft = 0;
1594  state.cntRight = d->count;
1595  state.thickness = 0;
1596  EvaluateCostFreeCut(state);
1597  int cntLeft = 0;
1598  int cntRight = d->count;
1599  int thickness = 0;
1600  // --------------------------------------
1601  float imult = (MaxPosition - MinPosition) / (float)bucketN;
1602  for (int i = 0; i < bucketN-1; i++) {
1603    cntLeft += state.bucketMin[i];
1604    cntRight -= state.bucketMax[i];
1605    thickness = cntLeft + cntRight - d->count;
1606    assert(thickness >= 0);
1607    state.cntLeft = cntLeft;
1608    state.cntRight = cntRight;
1609    state.thickness = thickness;
1610    state.position = MinPosition - bMinPosition + imult * ((float)(i+1));
1611    assert(state.position >= 0);
1612    assert(state.position <= state.width);
1613
1614    EvaluateCost(state);
1615  } // for
1616
1617  // free cut
1618  state.position = MaxPosition - bMinPosition;
1619  assert(state.position >= 0);
1620  assert(state.position <= state.width);
1621  state.cntLeft = d->count;
1622  state.cntRight = 0;
1623  state.thickness = 0;
1624  EvaluateCostFreeCut(state);
1625
1626  // we move back to the space
1627  state.bestPosition += bMinPosition;
1628 
1629  // In this case the best result is kept in 'state' variable
1630  return;
1631} // CTestAx::GetSplitPlaneOpt (for X, Y, Z)
1632
1633}
1634
Note: See TracBrowser for help on using the repository browser.