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

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

commit after merge with vlastimil

Line 
1// ===================================================================
2// $Id: $
3//
4// ktbai.cpp
5//     The classes for building up CKTB tree by statistical optimization.
6//     The bounding boxes of the all objects are sorted along all
7//     axes (x,y,z) in the array-based structure by radix sort algorithm.
8//
9// REPLACEMENT_STRING
10//
11// Copyright by Vlastimil Havran, 2007 - email to "vhavran AT seznam.cz"
12// Initial coding by Vlasta Havran, 2005.
13
14// standard headers
15#include <algorithm>
16#include <iostream>
17#include <string>
18
19// GOLEM headers
20#include "ktbai.h"
21#include "timer.h"
22#include "Environment.h"
23
24//#define _DEBUG
25
26#define AVOIDCHECKLIST
27
28#define _USE_OPTIMIZE_DIVIDE_AX
29
30// Note that the RADIX SORT does not work properly with -mtune=pentium3 for GCC.4.0 !!!!
31// Note that the RADIX SORT does not work properly with -mtune=pentium4 for GCC.3.4 !!!!
32// The reason is unknown, probably the bug in compiler
33
34namespace GtpVisibilityPreprocessor {
35
36// for testing accuracy of setting position from the ideal position
37//#define _RANDOMIZE_POSITION
38
39
40// This is epsilon for randomization
41static
42float randomEps = 0.0f; // By max 5 percent = 0.05
43
44static void
45RandomizePosition(float &value, float vmin, float vmax)
46{
47  assert(vmax > vmin);
48  assert( ((value >= vmin) && (value <= vmax) ) );
49
50  // copy the value
51  const float origvalue = value;
52  float offset = 0.f;
53
54#define _RND_OFFSET 
55#ifdef _RND_OFFSET
56  const int maxCount = 10;
57  int count = 0;
58  do {
59    if (count > maxCount) {
60      value = origvalue;
61      break;
62    }
63    offset = -(vmax - vmin) * randomEps/2.0f;
64    offset += randomEps * RandomValue(vmin,vmax);   
65    value = origvalue + offset;
66    count++;
67  }
68  while (! ((value > vmin) && (value < vmax)) );
69#else
70  offset = (vmax - vmin) * randomEps/2.0f;
71  if (RandomValue(0.f, 1.f) < 0.5f)
72    offset = -offset; 
73  value = origvalue + offset;
74  if ( (value < vmin) || (value > vmax) )
75    value = origvalue - offset;
76  //cout << "offset = " << 100.0f * offset / (vmax-vmin) << "\n";
77#endif
78
79  float epsShift = 1e-5;
80  if ( (value <= (vmin + (vmax-vmin) * epsShift)) ||
81       (value >= (vmax - (vmax-vmin) * epsShift)) ) {
82    value = origvalue;
83  }
84
85  assert( ((value >= vmin) && (value <= vmax) ) );
86 
87  return;
88}
89
90// ------------------------------------------------
91// for debugging cost function in the scene
92//#define _DEBUG_COSTFUNCTION
93#ifdef _DEBUG_COSTFUNCTION
94
95const int indexToFindSS = 0;
96
97static CFileIO ffcost;
98static int     costEvalCnt;
99static int     costCntObjects;
100static int     costIndexDebug;
101static string  fnameCost;
102static float   maxCost;
103static float   minCost;
104static float   maxPos;
105static float   minPos;
106static float   minPosAxis;
107static float   maxPosAxis;
108static bool _alreadyDebugged = false;
109static bool _streamOpened = false;
110
111static void
112InitCostStream(int indexDebug, int cntObjects,
113               float minx, float maxx)
114{
115  maxCost = -INFINITY;
116  minCost = INFINITY;
117  costEvalCnt = 0;
118  costIndexDebug = indexDebug;
119  costCntObjects = cntObjects;
120  minPosAxis = minx;
121  maxPosAxis = maxx;
122 
123  // and init the stream - rewrite
124  string name = "debugcostfile"
125  //Environment::GetSingleton()->GetBoolValue("OutputFileName", name);
126  char lns[100];
127  sprintf(lns, "%05d.res", indexDebug);
128  ChangeFilenameExtension(name, string(lns), fnameCost);
129
130  ffcost.SetFilename(fnameCost.c_str());
131  if (ffcost.OpenInMode(CFileIO::EE_Mode_w)) {
132    cerr << "Opening of stream failed" << endl;
133    abort();;
134  }
135  cout << "Saving debug to " << fnameCost << endl;
136  ffcost.WriteChars("\n");
137  sprintf(lns, "#CntObjects = %d\n", costCntObjects);
138  ffcost.WriteChars(lns);
139  sprintf(lns, "#IndexToDebug = %d\n", costIndexDebug);
140  ffcost.WriteChars(lns);
141  ffcost.WriteChars("#File ");
142  ffcost.WriteChars(fnameCost.c_str());
143  ffcost.WriteChars("\n");
144  ffcost.WriteChars("#==============================================\n");
145
146  cout << "#---------------------------------------" << endl;
147  _streamOpened = true;
148}
149
150static void
151CloseCostStream()
152{
153  if (_streamOpened) {   
154    char lns[255];
155    ffcost.WriteChars("#==============================================\n");
156    ffcost.WriteChars("#EndOfFile\n");
157    sprintf(lns, "#minCost = %8.6f minPos = %8.6f\n", minCost, minPos);
158    ffcost.WriteChars(lns);   
159    sprintf(lns, "#maxCost = %8.6f maxPos = %8.6f\n", maxCost, maxPos);
160    ffcost.WriteChars(lns);   
161    sprintf(lns, "#ratioMin/Max =                                %6.5f\n",
162            minCost/maxCost);
163    ffcost.WriteChars(lns);   
164    sprintf(lns, "#CntEval = %d\n", costEvalCnt);
165    ffcost.WriteChars(lns);
166    sprintf(lns, "#CntObjects = %d\n", costCntObjects);
167    ffcost.WriteChars(lns);   
168    sprintf(lns, "#IndexToDebug = %d\n", costIndexDebug);
169    ffcost.WriteChars(lns);
170    sprintf(lns, "#MinPosAxis = %8.6f maxPosAxis = %8.6f\n", minPosAxis, maxPosAxis);
171    ffcost.WriteChars(lns);
172    ffcost.WriteChars("#File ");
173    ffcost.WriteChars(fnameCost.c_str());
174    ffcost.WriteChars("\n");
175    ffcost.Close();
176    //cout << "#=======================================" << endl;
177    _streamOpened = false;
178    _alreadyDebugged = true;
179    // Temporarily
180    exit(0);
181  }
182}
183
184static void
185ReportCostStream(float pos, float cost)
186{
187  // Here we normalize a position to the interval <0-1>
188  float posn = (pos - minPosAxis)/(maxPosAxis - minPosAxis);
189
190  if (_streamOpened) {
191    costEvalCnt++;
192    if (cost < minCost) {
193      minCost = cost;
194      minPos = pos;
195    }
196    if (cost > maxCost) {
197      maxCost = cost;
198      maxPos = pos;
199    }
200    char lns[255];
201    sprintf(lns, "%8.6f %8.6f\n", posn, cost);
202    //cout << lns << endl;
203    ffcost.WriteChars(lns);   
204  }
205}
206#endif // _DEBUG_COSTFUNCTION
207
208//----------------------------------------------------------------
209//---------- Implementation of CKTB tree with irregular -----------
210//---------- placed splitting planes -----------------------------
211//----------------------------------------------------------------
212
213// default constructor
214CKTBABuildUp::CKTBABuildUp()
215{
216  // verbose is set
217  verbose = true;
218
219  // Maximum depth of the tree is set and stack is allocated
220  maxTreeDepth = MAX_HEIGHT;
221  stackDepth = maxTreeDepth + 2;
222  stackID = new GALIGN16 SInputData[stackDepth];
223  assert(stackID);
224  stackIndex = 0;
225  return;
226}
227
228//  destructor
229CKTBABuildUp::~CKTBABuildUp()
230{
231  delete [] stackID;
232  stackID = 0;
233}
234
235void
236CKTBABuildUp::ProvideID(ostream &app)
237{
238  bool tempvar;
239
240  Environment::GetSingleton()->GetBoolValue("BSP.splitclip", tempvar);
241  if (tempvar)
242    app << "#BSP.splitClip \t\tChange of bboxes(splitclip) "
243        << "dis/en\n" << tempvar << "\n";
244  Environment::GetSingleton()->GetBoolValue("BSP.emptyCut", tempvar);
245  app << "#BSP.emptyCut\t\tLate cutting off empty space in leaves"
246      << " dis/en\n" << tempvar << "\n";
247
248  if (tempvar) {
249    app << "#BSP.absMaxAllowedDepth\tMaximum abs depth for empty cutting\n";
250    app << absMaxAllowedDepth << "\n";
251    app << "#BSP.maxEmptyCutDepth\tMaximum allowed depth for late cutting"
252        << "(0-6)\n";
253    app << maxEmptyCutDepth << "\n";
254  }
255 
256  string termCrit;
257  Environment::GetSingleton()->GetStringValue("BSP.termCrit", termCrit);
258  app << "#BSP.termCrit\t\tTermination criteria to build up BSP tree\n";
259  app << termCrit << "\n";
260
261  maxCountTrials = 0;
262  if ( (!termCrit.compare("auto")) &&
263       (!termCrit.compare("auto2")) ) {
264    // everything except auto settings, when depth is determined
265    // It should be reported for adhoc and possibly others
266    Environment::GetSingleton()->GetIntValue("BSP.maxDepthAllowed",
267                                             maxDepthAllowed);
268    Environment::GetSingleton()->GetIntValue("BSP.maxListLength",
269                                             maxListLength);
270    app << "#MaxDepth (Dmax)\tMaximal allowed depth of the BSP tree\n";
271    app << maxDepthAllowed << "\n";
272    app << "#MaxListLength (Noilf)\tMaximal number of solids in one cell\n";
273    app << maxListLength << "\n";
274    // maximum number of trials to further subdivide
275    maxCountTrials = (int)(1.0 + 0.2 * (float)maxDepthAllowed);
276    if (maxCountTrials < 3)
277      maxCountTrials = 3;
278  }
279  else {
280    if (!termCrit.compare("auto2") ) {
281      Environment::GetSingleton()->GetIntValue("BSP.maxListLength",
282                                               maxListLength);
283      app << "#MaxDepth (Dmax)\tMaximal allowed depth of the BSP tree\n";
284    }
285  }
286 
287  float eCt, eCi, eCd;
288  Environment::GetSingleton()->GetFloatValue("BSP.decisionCost", eCd);
289  Environment::GetSingleton()->GetFloatValue("BSP.intersectionCost", eCi);
290  Environment::GetSingleton()->GetFloatValue("BSP.traversalCost", eCt);
291
292  // Ci should no be changed from 1.0, only Cd and Ct should be changed
293  app << "#SetBSP.(Cd, Ci, Ct) \tDecision, Intersection, ";
294  app << "and Traversal costs\n";
295  app << " ( " << eCd << ", " << eCi << " ," << eCt << " )\n";
296
297  return;
298}
299
300// -------------------------------------------------------
301// the next two axes for a given one
302const CKTBAxes::Axes
303CKTBABuildUp::oaxes[3][2]=
304{{CKTBAxes::EE_Y_axis, CKTBAxes::EE_Z_axis},
305 {CKTBAxes::EE_Z_axis, CKTBAxes::EE_X_axis},
306 {CKTBAxes::EE_X_axis, CKTBAxes::EE_Y_axis}};
307
308// -------------------------------------------------------
309// compare function passed to quicksort
310int
311CKTBABuildUp::Compare(const SItem *p, const SItem *q)
312{
313  register float pv = p->pos;
314  register float qv = q->pos;
315
316  if (pv < qv)
317    return 1;
318  else
319    if (pv > qv)
320      return -1;
321    else
322      // the coordinates are equal
323      if ( (p->IsRightBoundary()) &&
324           (q->IsLeftBoundary()) )
325        return 1; // right_boundary < left_boundary
326      else
327        if ( (p->IsLeftBoundary()) &&
328             (q->IsRightBoundary()) )
329          return -1; // left_boundary > right_boundary
330
331  return 0; // coordinates are equal, the same value and type
332}
333
334// --------------------------------------------
335// This is sorting using quicksort
336void
337CKTBABuildUp::SortOneAxis(SItemVec &itemvec, int cnt, int * const stack)
338{
339  // the pointer to the stack
340  int stackUk = 0;
341  // setting initial beg and end index for the whole array
342  stack[++stackUk] = 0;
343  stack[++stackUk] = cnt - 1;
344  SItem auxSwap;
345  int e, b, i, j;
346 
347  while ( stackUk !=0 ) {
348    // retrieving indeces from the stack
349    j = e = stack[stackUk--]; // begin and end in the array
350    i = b = stack[stackUk--]; // auxiliary indeces
351
352    // selecting pivot .. in the best case median
353    SItem X = itemvec[(i+j) / 2];
354    do {
355      while (Compare(&(itemvec[i]), &X) > 0)
356        i++; // find the array[i] > array[x]=pivot
357      while (Compare(&(itemvec[j]), &X) < 0)
358        j--; // find the array[j] < array[x]=pivot
359      if (i < j) { // swap the elements
360        auxSwap = itemvec[i];
361        itemvec[i] = itemvec[j];
362        itemvec[j] = auxSwap;
363        i++; j--;
364      }
365      else
366        if (i == j) {
367          i++;
368          j--;
369        }
370    }
371    while (i <= j);
372    // Now all the elements on the left are smaller than pivot and
373    // all the right are larger than pivot
374#if 0
375#ifdef _DEBUG
376    int t;
377    for (t = b; t < j; t++) {
378      if (X.pos < itemvec[t].pos) {
379        cout << "Problem 1 for t = " << t << endl;
380        cout << "X->pos = " << X.pos << "a[t]= " << itemvec[t].pos << endl;
381      }
382    }
383    for (t = i; t < e; t++) {
384      if (X.pos > itemvec[t].pos) {
385        cout << "Problem 2 for t = " << t << endl;
386        cout << "X->pos = " << X.pos << "a[t]= " << itemvec[t].pos << endl;
387      }
388    }
389#endif   
390#endif   
391    if (b < j) { // go to the left
392      stack[++stackUk] = b;
393      stack[++stackUk] = j;
394    }
395    if (i < e) { // go to the right
396      stack[++stackUk] = i;
397      stack[++stackUk] = e;
398    }
399    // Just a check
400    assert(stackUk <= cnt);
401  } // while (stackUk)
402
403  return;
404}
405
406void
407CKTBABuildUp::CopyToAuxArray(const SItemVec &bounds, SItemVecRadix &aux)
408{
409  // assume that both arrays are of the same size
410  assert(bounds.size() <= aux.capacity());
411
412#if 1
413  // source iterator
414  SItemVec::const_iterator src = bounds.begin();
415  for (SItemVecRadix::iterator it = aux.begin();
416       it != aux.end();
417       it++, src++)
418  {
419    *it = *src; // copy one elem
420  }
421#else
422  // Attempt for vectorization
423  int i;
424  const int cnt = aux.size();
425  for (i = 0; i < cnt; i++) {
426    aux[i] = bounds[i];
427  }
428#endif
429 
430  return;
431}
432
433// -----------------------------------------------------------------
434// The first pass of radix sort setting the links
435void
436CKTBABuildUp::RadixPassHoffset11(SItemVecRadix &bounds, int bit,
437                                 SRadix *rb, float offset,
438                                 SItemRadix **start)
439{
440  // preparing the hash table
441  memset(rb, 0, sizeof(SRadix) * RXBUFS30);
442
443  SItemRadix *p;
444  SItemRadix *q;
445  //SRadix *r;
446
447  // linearly process the input array
448  //for (SItemVecRadix::iterator it = bounds.begin(); it != bounds.end(); it++)
449  const int cnt = bounds.size();
450  for (int i = 0; i < cnt; i++)
451  {
452    // Take the element linearly from the array
453    //p = &(*it);
454    p = &(bounds[i]);
455
456    // add offset to get positive value to be sorted
457    p->pos += offset;
458    assert(p->pos >= 0.f);
459   
460    // compute the index of the bucket
461    unsigned int *val = (unsigned int*)&(p->pos);
462    unsigned int bucket = (( (*val) >> bit) & (RXBUFS30 - 1));
463    assert( bucket < RXBUFS30 );
464    // compute the address of correct bucket
465    //r = rb + bucket;
466    if (!rb[bucket].beg) {
467      rb[bucket].end = rb[bucket].beg = p;   // inserting into empty cell, mark begin
468      p->next = 0; // new item is the end of chain
469    }
470    else {
471      // if the end of bounding box
472      if (p->IsRightBoundary()) {
473        // right boundary - store item at the beginning of the group
474        p->next = rb[bucket].beg;
475        rb[bucket].beg = p;
476      }
477      else {
478        // left boundary - store item at the end of the group
479        p->next = 0; // it is the end elem
480        rb[bucket].end->next = p;   // chain last item to new item
481        rb[bucket].end = p;
482      }
483    }
484  } // for all input data
485
486  // append groups together
487  // q used as start item
488  SItemRadix **fset = &q;
489
490  for(int j = 0; j < RXBUFS30; j++) {
491    // only used elements
492    if (rb[j].beg) {
493      *fset = rb[j].beg;             // store pointer to the next group
494      fset = &(rb[j].end->next); // prepare group's end for update
495    }
496  } // for all buckets
497
498  *start = q;
499  return;
500}
501
502// radix sort pass starting with 'bit' over 'axis', using buckets in 'rb'
503void
504CKTBABuildUp::RadixPass11(SItemRadix **start, int cnt, int bit, SRadix *rb)
505{
506  // preparing the hash table
507  memset(rb, 0, sizeof(SRadix) * RXBUFS30);
508
509  // list grouping
510  SItemRadix *p = *start;
511  SItemRadix *q;
512  //SRadix *r;
513
514  // process the items using links
515  for (int i = 0; i < cnt; i++) {
516    assert(p);
517    unsigned int *val = (unsigned int*)&(p->pos);
518    unsigned int bucket = (( (*val) >> bit) & (RXBUFS30 - 1));
519    assert(bucket < RXBUFS30);
520    //r = rb + bucket;
521    if (!rb[bucket].beg)
522      rb[bucket].beg = p;   // inserting into empty cell, mark begin
523    else
524      rb[bucket].end->next = p;   // chain last item to new item
525
526    rb[bucket].end = p;  // update the end of group
527    q = p->next; // save the next position
528    p->next = 0; // new item is the end of chain
529    p = q;
530  } // for all items
531
532  // append groups together
533  //r = rb;
534  // q used as start item
535  SItemRadix **fset = &q;
536
537  for(int j = 0; j < RXBUFS30; j++)
538    // only used elements
539    if (rb[j].beg) {
540      *fset = rb[j].beg;       // store pointer to the next group
541      fset = &(rb[j].end->next); // prepare group's end for update
542    }
543
544  *start = q;
545  return;
546}
547
548// radix sort pass starting with 'bit' over 'axis', using buckets in 'rb'
549void
550CKTBABuildUp::RadixPassOffset10(SItemRadix **start, int cnt, int bit,
551                                SRadix *rb, float offset)
552{
553  // preparing the hash table
554  memset(rb, 0, sizeof(SRadix) * RXBUFS30_2);
555
556  // list grouping
557  SItemRadix *p = *start;
558  SItemRadix *q;
559  //SRadix *r;
560 
561  // process the items using links
562  for (int i = 0; i < cnt; i++) {
563    assert(p);
564    unsigned int *val = (unsigned int*)&(p->pos);
565    unsigned int bucket = (( (*val) >> bit) & (RXBUFS30_2 - 1));
566    assert(bucket < RXBUFS30_2);
567    //r = rb + bucket;
568    if (!(rb[bucket].beg)) {
569      rb[bucket].beg = p;   // inserting into empty cell, mark begin
570    }
571    else {
572      rb[bucket].end->next = p;   // chain last item to new item
573    }
574
575    rb[bucket].end = p;  // update the end of group
576
577    // Compute the correct value back, using the offset for the
578    // first pass
579    p->pos -= offset;
580
581    // go to the next item
582    q = p->next; // save the next position
583    p->next = 0; // new item is the end of chain
584    p = q; // take the next item
585  } // for all items
586
587  // append groups together and set bidirectional links
588  //r = rb;
589  SItemRadix *prev = 0;
590  for (int j = 0; j < RXBUFS30_2; j++) {
591    if (rb[j].beg) {
592      if (prev) {
593        // set forward link
594        prev->next = rb[j].beg;
595      }
596      prev = rb[j].end;
597    }
598  } // for i
599
600  // Find the first item in the list
601  //r = rb;
602  for (int k = 0; k < RXBUFS30_2; k++) {
603    if (rb[k].beg) {
604      // find the beginning of the chain
605      *start = rb[k].beg;
606      break;
607    }
608  } // for
609 
610  return;
611}
612
613void
614CKTBABuildUp::CopyFromAuxArray(SItemRadix *sorted, SItemVec &bounds)
615{
616  // source iterator
617  SItemRadix *src = sorted;
618  //for (SItemVec::iterator it = bounds.begin(); it != bounds.end(); it++)
619
620  int i;
621  const int cnt = bounds.size();
622  for (i = 0; i < cnt; i++)
623  {
624    bounds[i] = *src; // copy one elem, forgetting the pointer 'next'
625    src = src->next;
626  }
627
628  return;
629}
630
631
632// ---------------------------------------------------------------
633// sort the boundaries of objects' bounding boxes along all 3 axes
634void
635CKTBABuildUp::SortAxes(SInputData *data)
636{
637  // Note that cnt is the number of boundaries = odd number
638  assert(data);
639  assert(data->count);
640 
641  CTimer timer;
642  timer.Start();
643
644  // the number of boundaries
645  int cnt = 2*data->count;
646
647  if (!_useRadixSort) {
648    cout << "Sorting by quicksort " << flush;   
649#if 1
650    // this is the worst case allocation !!!
651    int *stackQuickSort = new GALIGN16 int [cnt+1];   
652    if (stackQuickSort == NULL) {
653      cerr << " ktbai.cpp: quicksort allocation failed\n";
654      exit(3);;
655    }
656    cout << " Qsort, sizeof(SItem)=" << sizeof(SItem)
657         << " size(vec[0])= " << sizeof(data->xvec[0]) << endl;
658    assert(cnt == data->xvec->size());
659    assert(cnt == data->yvec->size());
660    assert(cnt == data->zvec->size());
661    // Sort all three axis by Quicksort
662    //cout << "Sort x" << endl;
663    SortOneAxis(*(data->xvec), data->xvec->size(), stackQuickSort);
664    //cout << "Check x" << endl;
665    Check1List(data->xvec, 0, data->xvec->size()/2);
666    //cout << "Sort y" << endl;
667    SortOneAxis(*(data->yvec), data->yvec->size(), stackQuickSort);
668    //cout << "Check y" << endl;
669    Check1List(data->yvec, 0, data->yvec->size()/2);
670    //cout << "Sort z" << endl;
671    SortOneAxis(*(data->zvec), data->zvec->size(), stackQuickSort);
672    //cout << "Check z" << endl;
673    Check1List(data->zvec, 0, data->zvec->size()/2);
674 
675    // delete the data
676    delete [] stackQuickSort;
677#endif
678#if 0
679    sort(data->xvec->begin(), data->xvec->end());
680    sort(data->yvec->begin(), data->yvec->end());
681    sort(data->zvec->begin(), data->zvec->end());
682#endif
683#if 0
684    qsort(&(data->xvec[0]), data->xvec->size(), sizeof(data->xvec[0]),
685          (int(*)(const void*, const void*))(&Compare));
686    Check1List(data->xvec, 0, data->xvec->size()/2);
687   
688    qsort(&(data->yvec[0]), data->yvec->size(), sizeof(data->yvec[0]),
689          (int(*)(const void*, const void*))(&Compare));
690    qsort(&(data->zvec[0]), data->zvec->size(), sizeof(data->zvec[0]),
691          (int(*)(const void*, const void*))(&Compare));
692#endif
693  }
694  else {
695    // -----------------------------------------------
696    // Implementation by RadixSort with 3 passes, faster by 30% or so
697    // than 4-passes RadixSort, see ktbi.cpp
698    cout << "Sorting by radixsort3 " << flush;
699    SRadix *rb = new GALIGN16 SRadix[RXBUFS30];
700    assert(rb);
701    // the auxiliary array for sorting
702    SItemVecRadix *aux = new GALIGN16 SItemVecRadix();
703    assert(aux);
704    aux->resize(cnt);
705    // change axis
706    for(int i = 0; i < 3; i++) {     
707      // Offset to get only positive values to be sorted
708      float offset = -(wBbox.Min(i));
709      // Get a single array of bounds
710      SItemVec *bounds = data->GetItemVec(i);
711      assert(bounds);
712      // sort 3 times using 10 bits in one pass
713      CopyToAuxArray(*bounds, *aux);
714      // take the first elem and start radix sort
715      SItemRadix *start = 0;
716      // first pass: bits 0-10
717      RadixPassHoffset11(*aux, RXBITS30 * 0, rb, offset, &start);
718      // second pass: bits 11-21
719      RadixPass11(&start, cnt, RXBITS30 * 1, rb);
720      // third pass: bits 22-31 (bit 32 is sign)
721      RadixPassOffset10(&start, cnt, RXBITS30 * 2, rb, offset);
722      // Copy the sorted data back to the array, forgetting pointer 'next'
723      CopyFromAuxArray(start, *bounds);
724    }
725    // delete aux array
726    delete aux;
727    delete [] rb;
728    // end of radix sort
729  }
730
731  timer.Stop();
732  // ------------------------------------------
733  cout << " finished in " << timer.UserTime() << " [s] - user time " << endl;
734
735#ifdef _DEBUG
736  // check if the lists are correctly sorted .. in debug
737  // DEBUG << "SortAxis check\n" << flush;
738  Check3List(data);
739#endif
740
741  return;
742}
743
744// test if the lists are correctly sorted
745void
746CKTBABuildUp::Check3List(SInputData *data)
747{
748#ifndef AVOIDCHECKLIST
749#ifdef _DEBUG
750  assert(data);
751  if (data->xvec->size() != (unsigned int)data->count*2) {
752    cerr << "The number of X boundaries does not match" << endl;
753    abort();;
754  }
755  if (data->yvec->size() != (unsigned int)data->count*2) {
756    cerr << "The number of Y boundaries does not match" << endl;
757    abort();;
758  }
759  if (data->zvec->size() != (unsigned int)data->count*2) {
760    cerr << "The number of Z boundaries does not match" << endl;
761    abort();;
762  }
763  Check1List(data, 0, data->count);
764  Check1List(data, 1, data->count);
765  Check1List(data, 2, data->count);
766#endif // _DEBUG
767#endif // AVOIDCHECKLIST
768  return;
769}
770
771// test if one list is correctly sorted
772void
773CKTBABuildUp::Check1List(SInputData *data, int axis, int countExpected)
774{
775#ifndef AVOIDCHECKLIST
776#ifdef _DEBUG
777  assert((axis >= 0) && (axis < 3)); 
778  SItemVec *vec = data->GetItemVec(axis);
779  assert(vec);
780  if (countExpected)
781    Check1List(vec, axis, countExpected);
782  else {
783    if (vec->size() > 0) {
784      cerr << "The array of objects should be of zero size"<<endl;
785      abort();;
786    }
787  }
788#endif // _DEBUG
789#endif // AVOIDCHECKLIST
790}
791
792// test if one list is correctly sorted
793void
794CKTBABuildUp::Check1List(SItemVec *vec, int axis, int countExpected)
795{
796#ifndef AVOIDCHECKLIST
797#ifdef _DEBUG
798  int cntl,cntr; // the number of left/right boundaries
799 
800  // DEBUG << "Check1List starting\n";
801  cntl = 0; // the number of left boundaries
802  cntr = 0; // the number of right boundaries
803  int cntDups = 0; // the number of duplications for current position
804
805  // the position to check
806  SItemVec::iterator prev = vec->begin();
807 
808  // DEBUG << p->value[i] << "\n";
809  if ( (*prev).IsLeftBoundary()) {
810    cntl++;
811    cntDups++;
812  }
813  else {
814    cntr++;
815    cntDups--;
816  }
817  // the next item
818  SItemVec::iterator it = vec->begin();
819  it++;
820  // start from the second item
821  int index = 1;
822  bool firstDupsProb = false;
823  for ( ; it != vec->end(); it++, prev++, index++)
824  {
825    const SItem &p = *it;
826
827    if (p.obj->ToBeRemoved()) {
828      cout << "Object to be removed, obj = "
829           << p.obj << " -- it should not happen" << endl;
830    }
831   
832    // DEBUG << p->value[i] << "\n";
833    if (p.IsLeftBoundary()) {
834      cntl++;
835      cntDups++;
836    }
837    else {
838      cntr++;
839      cntDups--;
840    }
841
842    if (cntDups < 0) {
843      cerr << "The array is sorted in wrong way, cntDups = "
844            << cntDups << " index= " << index << " indexMax="
845            << countExpected*2 << endl;
846      if (!firstDupsProb) {
847        int imin = index-3;
848        if (imin < 0) imin = 0;
849        int imax = index+3;
850        if (imax > countExpected*2) imax = countExpected*2;
851        for (int i = imin; i < imax; i++) {
852          cout << "i= " << i << "  pos = " << (*vec)[i].pos;
853          if ((*vec)[i].IsLeftBoundary())
854            cout << "  L  ";
855          else
856            cout << "  R  ";
857          cout << (*vec)[i].obj << endl;
858        }
859      }     
860      firstDupsProb = true;
861    }
862   
863    if ( (*prev).pos > p.pos) {
864      cerr << "The array is for axis=" <<(int)axis<<" incorrectly sorted\n";
865      cerr << "prevPos = " << (*prev).pos << " curr->pos="
866            << p.pos << "\n";
867      abort();;
868    }
869    else {
870      if ( ((*prev).pos == p.pos) &&
871           ( (*prev).IsLeftBoundary()) &&
872           ( (*it).IsRightBoundary()) )  {
873        if (countExpected > 1) {
874          cerr << "The right and left boundary are incorrectly sorted, axis = "
875                << axis << "\n";
876          int i = 0;
877          for (SItemVec::iterator itt = vec->begin(); itt != vec->end(); itt++, i++)
878            {
879              cout << i << " = ";
880              if ((*itt).IsLeftBoundary())
881                cout << "L"; else cout << "R";   
882              cout << " obj = " << (*itt).obj << " pos = " << (*itt).pos << endl;
883            }
884        } // if more than one object
885      }
886    }
887  } // for all items
888
889  if (cntl != cntr) {
890    cerr << "The #left boundaries <> #right boundaries\n";
891    cerr << " #left boundaries= " << cntl << " #right boundaries= "
892          << cntr  << "  axis= " << (int)axis << "\n";
893    abort();;
894  }
895   
896  if (countExpected > 0) {
897    if (cntl != countExpected) {
898      cerr << "The number of found left boundaries = " << cntl
899            << " does not equal to expected count = " << countExpected << endl;
900    }
901    if (cntr != countExpected) {
902      cerr << "The number of found right boundaries = " << cntr
903            << " does not equal to expected count = " << countExpected << endl;
904    }
905  } // if checking of count is allowed
906 
907  cerr << flush;
908#endif // _DEBUG 
909#endif // AVOIDCHECKLIST
910}
911
912// Init required parameters
913void
914CKTBABuildUp::InitReqPref(SReqPrefParams *pars)
915{
916  // nothing is obligatory as default
917  pars->reqPosition = Limits::Infinity;
918  pars->reqAxis = CKTBAxes::EE_Leaf;
919  pars->useReqAxis = false;
920 
921  // the values for automatic termination criteria, since it was
922  // not subdivided
923  pars->ratioLast = 1000.0;
924  pars->ratioLastButOne = 1000.0;
925  pars->failedSubDivCount = 0;
926
927  Environment::GetSingleton()->GetIntValue("BSP.axisSelectionAlg",
928                                           _algorithmForAxisSelection);
929  if (_algorithmForAxisSelection != 0) {
930    // The axis as prefered one - testing all 3 axes. Start
931    // from the axis cutting the longest side of the box
932    pars->reqAxis = (CKTBAxes::Axes)(wBbox.Diagonal().DrivingAxis());
933  }
934
935  return;
936}
937
938// creates all the auxiliary structures for building up CKTB tree
939CKTBABuildUp::SInputData*
940CKTBABuildUp::Init(ObjectContainer *objlist, const AxisAlignedBox3 &box)
941{
942#ifdef  _RANDOMIZE_POSITION
943  //world->_environment.GetFloat("SetBSP.randomizePosition.Eps", randomEps);
944#endif // _RANDOMIZE_POSITION
945 
946#ifdef _KTB_CONSTR_STATS
947  _stats_interiorCount = 0;
948  _stats_bboxCount = 0;
949  _stats_leafNodeCount = 0;
950  _stats_emptyLeftNodeCount = 0;
951  // Aggregate statistics
952  _maxDepth = 0; 
953  _sumLeafDepth = 0;
954  _sumFullLeafDepth = 0;
955  _sumObjectRefCount = 0;
956  _maxObjectRefInLeaf = 0;
957  // surface areas
958  _sumSurfaceAreaLeaves = 0.f;
959  _sumSurfaceAreaMULcntLeaves = 0.f;
960  _sumSurfaceAreaInteriorNodes = 0.f;
961#endif
962  // If to print out the tree during contruction
963  Environment::GetSingleton()->GetBoolValue("BSP.printCuts", _printCuts);
964
965  // if to cut off empty space in postprocessing and how it is performed
966  Environment::GetSingleton()->GetBoolValue("BSP.emptyCut", cutEmptySpace);
967  Environment::GetSingleton()->GetBoolValue("BSP.useRadixSort", _useRadixSort);
968  Environment::GetSingleton()->GetIntValue("BSP.absMaxAllowedDepth",
969                                           absMaxAllowedDepth);
970  Environment::GetSingleton()->GetIntValue("BSP.maxEmptyCutDepth",
971                                           maxEmptyCutDepth);
972  Environment::GetSingleton()->GetIntValue("BSP.algAutoTermination",
973                                           algorithmAutoTermination);
974  Environment::GetSingleton()->GetFloatValue("BSP.biasFreeCuts",
975                                             biasFreeCuts);
976#ifdef _BIDIRLISTS
977  // if to clip the bounding boxes
978  Environment::GetSingleton()->GetBoolValue("BSP.splitclip", splitClip);
979#else
980  // Without bidirectional list we cannot make split clipping
981  splitClip = false;
982#endif // _BIDIRLISTS
983
984  // if to make tagging lists inside the tree
985  Environment::GetSingleton()->GetBoolValue("BSP.minBoxes.use",
986                                            makeMinBoxes);
987  Environment::GetSingleton()->GetBoolValue("BSP.minBoxes.tight",
988                                            makeTightMinBoxes);
989  Environment::GetSingleton()->GetIntValue("BSP.minBoxes.minObjects",
990                                            minObjectsToCreateMinBox);
991  Environment::GetSingleton()->GetIntValue("BSP.minBoxes.minDepthDistance",
992                                           minDepthDistanceBetweenMinBoxes);
993  Environment::GetSingleton()->GetFloatValue("BSP.minBoxes.minSA2ratio",
994                                             minSA2ratioMinBoxes);
995
996  // How many items can be allocated at once
997  int maxItemsAtOnce = 1;
998  if (makeMinBoxes) {
999#ifdef _KTB8Bytes   
1000    // We need to allocate for boxes the memory in a row
1001    maxItemsAtOnce = 5; // 8x5=40 = 16+24;
1002#else
1003    maxItemsAtOnce = 4; // 12x4=48 = 24+24;
1004#endif // _KTB8Bytes
1005    assert(minObjectsToCreateMinBox >= 1);
1006    assert(minDepthDistanceBetweenMinBoxes >= 0);
1007    assert(minDepthDistanceBetweenMinBoxes <= 50);
1008  }
1009
1010  // Initiate the allocator
1011  InitAux(0, MAX_HEIGHT - 1, maxItemsAtOnce);
1012 
1013  // since six planes are enough to cull empty space from leaves
1014  if ( (maxEmptyCutDepth < 0) ||
1015       (maxEmptyCutDepth > 6) ) {
1016    cerr << "BSP.maxEmptyCutDepth = " << maxEmptyCutDepth
1017          << " must be in range <0,6>\n";
1018    abort();;
1019  }
1020
1021  wBbox.Convert(box); // the box of the whole scene
1022  boxSize = box.Diagonal(); // the size of the box along the axes
1023 
1024  wholeBoxArea = wBbox.SA2(); // the half of the surface area of the box
1025 
1026  // the array of bounding boxes and duplications to the objects
1027  int objlistcnt = objlist->size();
1028
1029  // Here create the auxiliary array
1030  //solidArray.reserve(objlistcnt);
1031  // $$JB - crashed below
1032  solidArray.resize(objlistcnt);
1033
1034  // Here create the first entry of boundaries for all objects
1035  SInputData* data = AllocNewData();
1036  // set the box
1037  data->box = wBbox;
1038  // which algorithm to use for break-ax
1039#if 1
1040  data->algorithmBreakAx = 0; // position based
1041#else 
1042  data->algorithmBreakAx = 1; // pointer based
1043#endif
1044
1045#ifdef _RANDOMIZE_POSITION
1046  data->algorithmBreakAx = 0; // position based
1047#endif
1048
1049  resetFlagsForBreakAx = true;
1050 
1051  assert(objlistcnt > 0);
1052 
1053  // allocate the array of boundaries for all 3 axes
1054  data->Alloc(data->count*2);
1055  // set the number of objects
1056  data->count = objlistcnt;
1057  int i = 0;
1058  // cout << data->xvec->size()  << "  " << data->count*2 << endl;
1059  // create the list of bounding boxes
1060  for (ObjectContainer::iterator scr = objlist->begin();
1061       scr != objlist->end(); scr++, i++) {
1062    // set the object
1063    solidArray[i].obj = *scr;
1064    solidArray[i].flags = 0;
1065    SSolid *solid = &(solidArray[i]);
1066    // create bounding box of the object
1067    SBBox abox;
1068    abox.Convert((*scr)->GetBox());   
1069    // copy the boundaries to the arrays
1070    // max boundaries as first because of stable quicksort
1071    SItem itemX;
1072    itemX.pos = abox.Min(0);
1073    itemX.obj = solid;
1074    itemX.axis = 0;
1075    itemX.SetLeftBoundary();
1076    data->xvec->push_back(itemX);
1077
1078    SItem itemY;
1079    itemY.pos = abox.Min(1);
1080    itemY.obj = solid;
1081    itemY.axis = 1;
1082    itemY.SetLeftBoundary();
1083    data->yvec->push_back(itemY);
1084
1085    SItem itemZ;
1086    itemZ.pos = abox.Min(2);
1087    itemZ.obj = solid;
1088    itemZ.axis = 2;
1089    itemZ.SetLeftBoundary();
1090    data->zvec->push_back(itemZ);
1091
1092    // max boundaries
1093    itemX.pos = abox.Max(0);
1094    itemX.SetRightBoundary();
1095    data->xvec->push_back(itemX);
1096
1097    itemY.pos = abox.Max(1);
1098    itemY.SetRightBoundary();
1099    data->yvec->push_back(itemY);
1100
1101    itemZ.pos = abox.Max(2);
1102    itemZ.SetRightBoundary();
1103    data->zvec->push_back(itemZ);
1104  } // for all objects
1105
1106  // cout << data->xvec->size()  << "  " << data->count*2 << endl;
1107 
1108  assert(data->xvec->size() == (unsigned int)data->count*2);
1109  data->xvec->resize(data->count*2);
1110  assert(data->yvec->size() == (unsigned int)data->count*2);
1111  data->yvec->resize(data->count*2);
1112  assert(data->zvec->size() == (unsigned int)data->count*2);
1113  data->zvec->resize(data->count*2);
1114 
1115  // If to use radix sort or quicksort
1116  //_useRadixSort = true;
1117  // Note that quicksort does not work for some reason !! 27/12/2005
1118  //_useRadixSort = false;
1119
1120  // sort all the bounding boxes along the axes
1121  if ((initcnt = objlist->size()) != 0)
1122    SortAxes(data);
1123
1124  // The statistics init
1125  cntDuplicate = 0;
1126
1127  // Init the parameters from the environment file
1128  InitReqPref(&(data->pars));
1129
1130  // Init the termination criteria
1131  string termCrit;
1132  Environment::GetSingleton()->GetStringValue("BSP.termCrit", termCrit);
1133
1134  maxCountTrials = 0;
1135  if (!termCrit.compare("adhoc")) {
1136    // termination criteria are the max depth of the hierarchy, number
1137    // of primitives
1138    data->modeSubDiv = EE_SUBDIVADHOC;
1139    Environment::GetSingleton()->GetIntValue("BSP.maxListLength",
1140                                             maxListLength);
1141    Environment::GetSingleton()->GetIntValue("BSP.maxDepthAllowed",
1142                                             maxDepthAllowed);
1143  }
1144  else {
1145    if ( (!termCrit.compare("auto")) ||
1146         (!termCrit.compare("auto2")) ) {
1147      // automatic termination criteria are used, everything is computed
1148      // automatically from number of objects etc.
1149      int estHeight = (int)(log((float)initcnt)/log((float)2.0) + 0.9);
1150      // cout << "EstHeight=" << estHeight << endl;
1151      maxDepthAllowed = (int)((float)estHeight * 1.2f + 2.0f);
1152
1153      // maximum number of trials to further subdivide
1154      maxCountTrials = (int)(1.0 + 0.2 * (float)maxDepthAllowed);
1155      if (maxCountTrials < 3)
1156        maxCountTrials = 3;
1157     
1158      // for 'auto2' we specify the length of the object list by hand
1159      if (!termCrit.compare("auto2"))
1160        Environment::GetSingleton()->GetIntValue("BSP.maxListLength",
1161                                                 maxListLength);
1162      else
1163        // for 'auto' we set the length of the object list that is supposed
1164        // to be the best for ray shooting
1165        maxListLength = 1;
1166
1167      if (verbose) {
1168        cout << "TERMCRIT:maximum height of a leaf set to "
1169               << maxDepthAllowed << "\n";
1170        cout << "TERMCRIT:maximum number of objects in leaf set to "
1171               << maxListLength << "\n" ;
1172      }
1173      // set the mode to be recognized
1174      data->modeSubDiv = EE_SUBDIVAUTOMATIC;
1175    }
1176    else {
1177      cerr << " unknown termination criteria for BSP tree\n";
1178      cerr << "It was specified: " << termCrit << "\n";
1179      cerr << "Allowed types = auto, auto2, adhoc" << endl;
1180      exit(3);;
1181    }
1182  }
1183
1184  // for some evaluation it is necessary to determine the following costs
1185  Environment::GetSingleton()->GetFloatValue("BSP.traversalCost", Ct);
1186  Environment::GetSingleton()->GetFloatValue("BSP.intersectionCost", Ci);
1187  if (verbose)
1188    cout << "Ct=" << Ct << " Ci=" << Ci << "\n";
1189 
1190  if (cutEmptySpace) {
1191    // correct the maximum abs depth, when late cutting is allowed
1192    if (absMaxAllowedDepth < maxDepthAllowed)
1193      absMaxAllowedDepth = maxDepthAllowed;
1194    if (absMaxAllowedDepth > (maxDepthAllowed + maxEmptyCutDepth) )
1195      absMaxAllowedDepth = maxDepthAllowed + maxEmptyCutDepth;
1196    if (absMaxAllowedDepth >= MAX_HEIGHT) {
1197      absMaxAllowedDepth = MAX_HEIGHT;
1198      maxDepthAllowed = MAX_HEIGHT - maxEmptyCutDepth;
1199    }
1200   
1201    if (verbose) {
1202      cout << "Cutting off empty spaces ON\n";
1203      cout << "BSP.absMaxAllowedDepth = " << absMaxAllowedDepth << "\n";
1204      cout << "BSP.maxEmptyCutDepth = " << maxEmptyCutDepth << endl;
1205    }
1206  }
1207  else {
1208    if (maxDepthAllowed >= MAX_HEIGHT)
1209      maxDepthAllowed = MAX_HEIGHT - 1;
1210  }
1211
1212  if (verbose)
1213    cout << "MaxListLength = " << maxListLength << endl; 
1214
1215  return data;
1216}
1217
1218// interface function for building up CKTB tree
1219SKTBNodeT*
1220CKTBABuildUp::BuildUp(ObjectContainer &objlist,
1221                      const AxisAlignedBox3 &box,
1222                      bool verboseF)
1223{
1224  // check the number of objects
1225  if (objlist.size() == 0)
1226    return 0; // nothing
1227
1228  // ---------------------------------------------------
1229  // Rewriting the triangles to the just wrappers to save
1230  // the memory during building!!!!
1231  bool useWrappers = false;
1232  if (useWrappers) {
1233    cout << "WARNING: using only wrappers, not objects to save the memory!"
1234         << endl;
1235    cout << "Size of(AxisAlignedBox3Intersectable) = " << sizeof(AxisAlignedBox3Intersectable) << endl;
1236    cout << "Size of(TriangleIntersectable) = " << sizeof(TriangleIntersectable) << endl;
1237    for (ObjectContainer::iterator it = objlist.begin();
1238         it != objlist.end(); it++) {
1239      // take the triangle
1240      Intersectable *i = *it;
1241      // store the properties to new variables
1242      AxisAlignedBox3 b = i->GetBox();
1243      int mId = i->mId;
1244      // delete the triangle
1245      delete i;
1246      // create the wrapper with the same box as triangle
1247      AxisAlignedBox3Intersectable *a = new AxisAlignedBox3Intersectable(b);
1248      a->mId = mId;
1249      // and put it back to the list of objects
1250      *it = a;
1251    } // for   
1252    cout << "Rewriting to wrappers is finished!" << endl;
1253  } // if
1254  // ---------------------------------------------------
1255 
1256  // cerr<<"hh44"<<endl;
1257  // initialize the whole box of the kd-tree and sort
1258  // the boundary entries
1259  SInputData *d = Init(&objlist, box);
1260
1261  //cerr<<"hh45"<<endl;
1262
1263  // copy verbose to be used consistenly further on
1264  this->verbose = verboseF;
1265 
1266  if (verbose)
1267    cout << "Building up KTB tree for " << objlist.size()
1268                 << " objects." << endl << flush;
1269
1270  // Initialization of the first value to insert the min boxes
1271  d->lastMinBoxSA2 = box.SurfaceArea() * 10000.f;
1272  d->lastDepthForMinBoxes = -20; // for root you always put the node
1273  d->lastMinBoxNode = 0;
1274 
1275  // -----------------------------------------------
1276  // Start to build the tree
1277  if ( (cutEmptySpace) &&
1278       (initcnt <= maxListLength)) {
1279    // only cutting off empty space for initial leaf
1280    d->modeSubDiv = EE_SUBDIVCUTEMPTY;
1281    root = SubDiv(d);
1282  }
1283  else {
1284    // normal subdivision scheme, creating whole tree
1285    root = SubDiv(d);
1286  }
1287
1288#ifdef _DEBUG
1289  if (GetDepth() != 0) {
1290    cerr << "Using depth value does not work correctly, depth = "
1291          << GetDepth() << endl;
1292  }
1293#endif 
1294  assert(root); 
1295
1296  // ----------------------------------------------------
1297  // Deallocate all auxiliary data
1298  DeleteAuxiliaryData();
1299
1300  // delete the temporary data array
1301  solidArray.resize(0);
1302 
1303  // the pointer to the root node
1304  return GetRootNode();
1305}
1306
1307int
1308CKTBABuildUp::UpdateEvaluation(float &eval, const float &newEval)
1309{
1310  if (newEval < eval) {
1311    eval = newEval;
1312    return 1; // updated
1313  }
1314  return 0; // not updated
1315}
1316
1317// recursive function for subdividing the CKTB tree into two
1318// halves .. creates one node of CKTB tree
1319// list .. list of boundaries, count .. number of objects in current node
1320// bb .. current bounding box of node, (pars,reqGlobAxis, modeSubDiv) .. cut pref.
1321SKTBNodeT*
1322CKTBABuildUp::SubDiv(SInputData *d)
1323{
1324#ifdef _DEBUG
1325  Check3List(d);
1326#endif
1327
1328#if 0
1329  static int index = 0;
1330  cout << "SubDiv index = " << index << endl;
1331  const int indexToFind = 13916;
1332  if (index == indexToFind) {
1333    cout << " IndexToFind = " << indexToFind << endl;
1334    cout << " You should start debug" << endl;
1335  }
1336  index++;
1337#endif
1338
1339  // This is the node to return from this function
1340  SKTBNodeT *nodeToReturn = 0;
1341
1342  // if to make the interior node extended by the box here
1343  makeMinBoxHere = false;
1344
1345  // -----------------------------------------------------
1346  if (makeMinBoxes) {
1347#if 1
1348    if ( (d->count >= minObjectsToCreateMinBox) &&
1349         (GetDepth() - d->lastDepthForMinBoxes >= minDepthDistanceBetweenMinBoxes) )  {
1350        makeMinBoxHere = true;
1351        d->lastDepthForMinBoxes = GetDepth();
1352        d->lastMinBoxSA2 = d->box.SA2();
1353    }
1354#endif
1355#if 0
1356    if (d->count >= minObjectsToCreateMinBox) {
1357      if ( (d->lastMinBoxSA2 >= d->box.SA2() * minSA2ratioMinBoxes) &&
1358           (GetDepth() - d->lastDepthForMinBoxes >= minDepthDistanceBetweenMinBoxes) ) {
1359        makeMinBoxHere = true;
1360        d->lastDepthForMinBoxes = GetDepth();
1361        d->lastMinBoxSA2 = d->box.SA2();
1362      }
1363    }
1364#endif
1365#if 0
1366    if ( (d->count >= minObjectsToCreateMinBox) &&
1367         (GetDepth() % minDepthDistanceBetweenMinBoxes == 0) ) {
1368      makeMinBoxHere = true;
1369      d->lastDepthForMinBoxes = GetDepth();
1370      d->lastMinBoxSA2 = d->box.SA2();
1371    }
1372#endif
1373  }
1374
1375  if ( (makeMinBoxHere) && (makeTightMinBoxes) ) {
1376    // In case we should make the box of the current
1377    // interior node and this should be tight, we have
1378    // to first compute its extent before splitting
1379    SBBox tightbox;
1380    GetTightBox(*d, tightbox);
1381    // Here we decrease the box to the minimum size
1382    d->box = tightbox;
1383  }
1384 
1385  // ------------------------------------------------------------------
1386  // if we are forced to make subdivision at given point
1387  if (d->pars.reqPosition != Limits::Infinity) {
1388    // this code preceeds switch(mode) since 2-level evaluation
1389    d->position = d->pars.reqPosition;
1390    d->axis = d->pars.reqAxis;
1391    // next subdivison will not be forced
1392    d->pars.reqPosition = Limits::Infinity;
1393    d->pars.reqAxis = CKTBAxes::EE_Leaf;
1394    IncDepth();
1395    SKTBNodeT* n = MakeOneCut(d);
1396    if (!nodeToReturn)
1397      nodeToReturn = n;
1398    DecDepth();
1399    return nodeToReturn;
1400  }
1401 
1402  // -------------------------------------------------------------------------
1403  // determine between subdivision modes - different termination criteria
1404  switch (d->modeSubDiv) { // subdivision mode
1405    case EE_SUBDIVCUTEMPTY: { // cutting off empty space
1406      if ((GetDepth() >= absMaxAllowedDepth) ||
1407          (GetDepth() >= (startEmptyCutDepth + maxEmptyCutDepth)) ||
1408          (d->count == 0) ) {
1409        SKTBNodeT* n = MakeLeaf(d);
1410        if (!nodeToReturn)
1411          nodeToReturn = n;
1412        // cout << "LEC\n";
1413        return nodeToReturn;
1414      }
1415      // to require the axis has no meaning in cutting off
1416      d->pars.reqAxis = d->axis = CKTBAxes::EE_Leaf;
1417      break;
1418    }
1419    // both adhoc, clustering and automatic termination criteria
1420    case EE_SUBDIVADHOC:
1421    case EE_SUBDIVAUTOMATIC: {
1422      // automatic termination criteria are used in this phase as
1423      // the absolute threshold similarly to SUBDIVADHOC mode
1424      if ( (GetDepth() >= maxDepthAllowed) ||
1425           (d->count <= maxListLength)) {
1426        SKTBNodeT* n = MakeLeaf(d);
1427        if (!nodeToReturn)
1428          nodeToReturn = n;
1429        // cout << "LLA\n";
1430        return nodeToReturn;
1431      }
1432      break;
1433    }
1434    default: {
1435      cerr << " unknown subdivision mode in ibsp.cpp\n";
1436      abort();;
1437    }
1438  } // switch(modeSubDiv)
1439 
1440  // the axis where splitting should be done
1441  CKTBAxes::Axes axis = CKTBAxes::EE_Leaf;
1442  d->twoSplits = -1;
1443  d->bestCost = WorstEvaluation() * 0.5f;
1444  d->position = MAXFLOAT;
1445  d->cntThickness = 0;
1446
1447// which calls should be used, standard or optimized
1448#if 1
1449#define CallGetSplitPlaneX GetSplitPlaneOpt2
1450#define CallGetSplitPlaneY GetSplitPlaneOpt2
1451#define CallGetSplitPlaneZ GetSplitPlaneOpt2
1452#endif
1453#if 0
1454// This does not work for unknown reason... VH 2006/01/08
1455#define CallGetSplitPlaneX GetSplitPlaneOpt3
1456#define CallGetSplitPlaneY GetSplitPlaneOpt3
1457#define CallGetSplitPlaneZ GetSplitPlaneOpt3
1458#endif
1459#if 0
1460#define CallGetSplitPlaneX GetSplitPlaneOptUnroll4
1461#define CallGetSplitPlaneY GetSplitPlaneOptUnroll4
1462#define CallGetSplitPlaneZ GetSplitPlaneOptUnroll4
1463#endif
1464 
1465  // if the axis is required in pars structure for some reasons
1466  if (d->pars.useReqAxis) {
1467    axis = d->pars.reqAxis;
1468    // The next subdivision is not prescribed
1469    d->pars.useReqAxis = false;
1470    switch(axis) {
1471      case CKTBAxes::EE_X_axis : {
1472        state.InitXaxis(d->count, d->box); // init
1473        CallGetSplitPlaneX(d->xvec, 0); // evaluate
1474        break;
1475      }
1476      case CKTBAxes::EE_Y_axis : {
1477        state.InitYaxis(d->count, d->box); // init
1478        CallGetSplitPlaneY(d->yvec, 1); // evaluate
1479        break;
1480      }
1481      case CKTBAxes::EE_Z_axis : {
1482        state.InitZaxis(d->count, d->box); // init
1483        CallGetSplitPlaneZ(d->zvec, 2); // evaluate
1484        break;
1485      }
1486      default: {
1487        cerr << "No other option is allowed here" << endl;
1488        abort();;
1489      }
1490    }
1491    // compute the cost, normalized by surface area
1492    state.NormalizeCostBySA2();
1493    d->bestCost /= state.areaWholeSA2;
1494    if (UpdateEvaluation(d->bestCost, state.bestCost)) {
1495      // Compute correct position to be used for splitting
1496      d->position = (*(state.bestIterator)).pos;
1497      d->bestIterator = state.bestIterator;
1498      d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
1499      d->cntThickness = state.bestThickness;
1500      if (d->twoSplits > 1) {
1501        SItemVec::iterator it = state.bestIterator; it++;     
1502        if (it != d->GetItemVec(axis)->end())
1503          d->position2 = (*it).pos;
1504        else
1505          d->position2 = state.box.Max(axis);
1506      }
1507    }
1508  }
1509  else {
1510    int algorithmForAxisSel = _algorithmForAxisSelection;
1511    if (d->modeSubDiv == EE_SUBDIVCUTEMPTY)
1512      algorithmForAxisSel = 0;
1513
1514    // only a single axis will be tested
1515    bool useSingleAxis = true;
1516       
1517    // the required axis is by global function .. e.g. cyclic change x, y, z
1518    switch (algorithmForAxisSel) {
1519      case 0: {
1520        // all three axis will be used, starting from x, y, z
1521        useSingleAxis = false;
1522        break;
1523      }
1524      case 1: {
1525        // cyclic order of axes, x, y, z, x, y....
1526        axis = d->pars.reqAxis;
1527        // compute the axis for the next subdivision step
1528        d->pars.reqAxis = oaxes[axis][0];
1529        break;
1530      }
1531      case 2 : {
1532        // The next time will be determined the same way
1533        axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
1534        break;
1535      }
1536      case 3 : {
1537        float p = RandomValue(0.0f, 1.0f);
1538        const float thresholdAxisP = 0.8f;
1539        if (p < thresholdAxisP) {         
1540          axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
1541        }
1542        else {
1543          // Use the axis prescribed from previous step
1544          axis = d->pars.reqAxis;
1545        }
1546        // for the next time, different axis
1547        d->pars.reqAxis = oaxes[axis][0];
1548        break;
1549      }
1550      case 4 : {
1551        float p = RandomValue(0.0f, 1.0f);
1552        const float thresholdAxisP = 0.3f;
1553        if (p < thresholdAxisP) {         
1554          axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
1555        }
1556        else {
1557          // Use the axis prescribed from previous step
1558          axis = d->pars.reqAxis;
1559        }
1560        // for the next time, different axis
1561        d->pars.reqAxis = oaxes[axis][0];
1562        break;
1563      }
1564      case 5 : {
1565        float p = RandomValue(0.0f, 1.0f);
1566        const float thresholdAxisP = 0.2f;
1567        d->pars.reqAxis = (CKTBAxes::Axes)-1;
1568        if (p < thresholdAxisP) {
1569          axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
1570          // for the next time, different axis
1571        }
1572        else {
1573          // Use the axis for which cost model gets minimum,
1574          // compute it for all axes x, y, z
1575          useSingleAxis = false;
1576        }
1577        break;
1578      }
1579      default: {
1580        cerr << "Selection algorithm = " << algorithmForAxisSel
1581              << " for axis is not implemented" << endl;
1582        abort();;
1583      }
1584    } // switch
1585
1586    // How the algorithm is used
1587    if (useSingleAxis) {
1588      // -------------------------------------------
1589      // Compute subdivision only for one axis
1590      switch(axis) {
1591        case CKTBAxes::EE_X_axis : {
1592          state.InitXaxis(d->count, d->box); // init
1593          CallGetSplitPlaneX(d->xvec, 0); // evaluate
1594          break;
1595        }
1596        case CKTBAxes::EE_Y_axis : {
1597          state.InitYaxis(d->count, d->box); // init
1598          CallGetSplitPlaneY(d->yvec, 1); // evaluate
1599          break;
1600        }
1601        case CKTBAxes::EE_Z_axis : {
1602          state.InitZaxis(d->count, d->box); // init
1603          CallGetSplitPlaneZ(d->zvec, 2); // evaluate
1604          break;
1605        }
1606        //case CKTBAxes::EE_Leaf : {
1607        //goto MAKE_LEAF;
1608        //}
1609        default: {
1610          cerr << "Selection algorithm = " << algorithmForAxisSel
1611                << " for axis = " << axis <<" is not implemented" << endl;
1612          abort();;
1613        }
1614      }
1615      state.NormalizeCostBySA2();
1616      d->bestCost /= state.areaWholeSA2;
1617      if (UpdateEvaluation(d->bestCost, state.bestCost)) {
1618        // Compute correct position to be used for splitting
1619        d->position = (*(state.bestIterator)).pos;
1620        d->bestIterator = state.bestIterator;
1621        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
1622        d->cntThickness = state.bestThickness;
1623        if (d->twoSplits > 1) {
1624          SItemVec::iterator it = state.bestIterator; it++;     
1625          if (it != d->GetItemVec(axis)->end())
1626            d->position2 = (*it).pos;
1627          else
1628            d->position2 = state.box.Max(axis);
1629        }
1630      }
1631    }
1632    else {
1633      // no axis is required, we can pickup any we like. We test all three axis
1634      // and we pickup the minimum cost for all three axes.
1635      state.InitXaxis(d->count, d->box); // init for x axis
1636      CallGetSplitPlaneX(d->xvec, 0); // evaluate for x axis
1637      if (UpdateEvaluation(d->bestCost, state.bestCost)) {     
1638        axis = CKTBAxes::EE_X_axis; // the x-axis should be used
1639        // Compute correct position to be used for splitting
1640        d->bestCost = state.bestCost;
1641        d->position = (*(state.bestIterator)).pos;
1642        d->bestIterator = state.bestIterator;
1643        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
1644        d->cntThickness = state.bestThickness;
1645        if (d->twoSplits > 1) {
1646          SItemVec::iterator it = state.bestIterator;
1647          it++;
1648          if (it != d->xvec->end())
1649            d->position2 = (*it).pos;
1650          else
1651            d->position2 = state.box.Max().x;
1652        }
1653      }
1654      // -----------------------
1655      // now test for y axis
1656      state.ReinitYaxis(d->count, d->box); // init for x axis
1657      //state.InitYaxis(d->count, d->box); // init for x axis
1658      CallGetSplitPlaneY(d->yvec, 1); // evaluate for x axis
1659      if (UpdateEvaluation(d->bestCost, state.bestCost)) {
1660        axis = CKTBAxes::EE_Y_axis; // the y-axis should be used
1661        // Compute correct position to be used for splitting
1662        d->position = (*(state.bestIterator)).pos;
1663        d->bestIterator = state.bestIterator;
1664        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
1665        d->cntThickness = state.bestThickness;
1666        if (d->twoSplits > 1) {
1667          SItemVec::iterator it = state.bestIterator;
1668          it++;
1669          if (it != d->yvec->end())
1670            d->position2 = (*it).pos;
1671          else
1672            d->position2 = state.box.Max().y;
1673        }
1674      }
1675      // -----------------------
1676      // now test for z axis
1677      state.ReinitZaxis(d->count, d->box); // init for x axis
1678      //state.InitZaxis(d->count, d->box); // init for x axis
1679      CallGetSplitPlaneZ(d->zvec, 2); // evaluate for x axis
1680      if (UpdateEvaluation(d->bestCost, state.bestCost)) {
1681        axis = CKTBAxes::EE_Z_axis; // the z-axis should be used
1682        // Compute correct position to be used for splitting
1683        d->position = (*(state.bestIterator)).pos;
1684        d->bestIterator = state.bestIterator;
1685        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
1686        d->cntThickness = state.bestThickness;
1687        if (d->twoSplits > 1) {
1688          SItemVec::iterator it = state.bestIterator;
1689          it++;
1690          if (it != d->zvec->end())
1691            d->position2 = (*it).pos;
1692          else
1693            d->position2 = state.box.Max().z;
1694        }
1695      }
1696      // Now we have to renormalize the cost for purpose of termination the building
1697      d->bestCost /= state.areaWholeSA2;
1698    } // arbitrary order of axes
1699  } // for which axis to compute surface area heuristics
1700
1701  // ----------------------------------------------------------
1702  // when automatic termination criteria should be used
1703  if (d->modeSubDiv == EE_SUBDIVAUTOMATIC) {
1704    bool stopSubdivision = false;
1705    float ratio;
1706    // which algorithm to use for automatic termination criteria
1707    switch (algorithmAutoTermination) {
1708      // --------------------------------------------------
1709      case 0: { // This is described in my thesis.
1710        // This is the cost of unsubdivided node
1711        float leafCost = (float)(d->count);
1712        ratio = d->bestCost / leafCost;
1713        const float minRatio = 0.75;
1714        if (ratio > minRatio)
1715          stopSubdivision = true;
1716        break;
1717      }
1718      // --------------------------------------------------
1719      case 1: {
1720        // This is the cost of unsubdivided node, with uniformly distributed
1721        // objects and the spatial median subdivision, without any object straddling
1722        // the splitting plane. In some sense ideal setting for uniform distribution.
1723        Vector3 s = d->box.Diagonal();
1724        // Taking the widest side of the object to be subdivided
1725        int diagAxis = s.DrivingAxis();
1726        // Half of the width
1727        // ----------------------
1728        float w2 = s[diagAxis] * 0.5;
1729        float t = s[oaxes[diagAxis][0]];
1730        float h = s[oaxes[diagAxis][1]];
1731        float sah2 = w2*(t+h) + h*t; // half of the surface area of a child
1732        // This is the cost for case when the node is subdivided in the half and
1733        // the half of the object is on the left and half of the objects on the right
1734        // This is not the best cost possible !!!
1735        // = float subdividedCost = 2 * (sah2/state.areaWholeSA2 * count/2);
1736        float subdividedCost = sah2/state.areaWholeSA2 * (float)(d->count);
1737        ratio = d->bestCost / subdividedCost;
1738        // This is the max ratio allowed for splitting node subdivision
1739
1740        // The computed ratio = 1.0 corresponds to 'ideal case', so the threshold
1741        // must be higher than 1.0 !
1742        const float minRatio = 1.1;
1743        if (ratio > minRatio)
1744          stopSubdivision = true;
1745        break;
1746      }
1747      default : {
1748        cerr << "Uknown algorithm for automatic termination criteria = "
1749              << algorithmAutoTermination << endl;
1750      }
1751       
1752    } // switch
1753   
1754    //cout << "R=" << ratio << " ";
1755    if (stopSubdivision) {
1756      // cout << "F" << pars.failedSubDivCount << " ";
1757      // when small improvement by the subdivision is reached
1758      d->pars.failedSubDivCount++;
1759      if (d->pars.failedSubDivCount > maxCountTrials) {
1760        // make leaf and finish with this KD-tree branch
1761        // cout << "L, cnt=" << count << " ,d=" << GetDepth()()
1762        //     << " ,c=" << bestQ << " ,rat=" << ratio << "\n";
1763        // possibly cut empty space before leaf is created
1764        if (cutEmptySpace) {
1765          // setting initial depth if empty space cutting
1766          startEmptyCutDepth = GetDepth();
1767          IncDepth();
1768          d->modeSubDiv = EE_SUBDIVCUTEMPTY;
1769          // create the leaf first, with late empty cutting
1770          SKTBNodeT *n = SubDiv(d);
1771          if (!nodeToReturn)
1772            nodeToReturn = n;
1773          DecDepth();
1774          return nodeToReturn; // and finish by constructing a leaf
1775        }
1776
1777        // make leaf and finish
1778        SKTBNodeT *n = MakeLeaf(d);
1779        if (!nodeToReturn)
1780          nodeToReturn = n;
1781        return nodeToReturn;
1782      }
1783    }
1784
1785    // Do not stop subdivision now
1786   
1787    // update the progress in the parameters
1788    d->pars.ratioLastButOne = d->pars.ratioLast;
1789    d->pars.ratioLast = ratio;
1790
1791    // assure that axis is defined, if finding the
1792    // splitting plane has failed
1793    if ( (axis == CKTBAxes::EE_Leaf) ||
1794         (d->position == MAXFLOAT) ) {
1795      // compute the spatial median for arbitrary axis
1796      axis = (CKTBAxes::Axes)int(RandomValue(0.f, 2.98f));
1797      // set position based algorithm for the next cut and children
1798      d->algorithmBreakAx = 0;
1799      resetFlagsForBreakAx = true;
1800      d->position = (d->box.Min(axis) + d->box.Max(axis)) * 0.5f;
1801      d->twoSplits = 1;
1802    }
1803  } // subdivision automatic
1804
1805  // -----------------------------------------------
1806  // We have evaluated surface area heuristics above
1807  // if no winner for subdivision exists, make a leaf
1808  if ( (axis == CKTBAxes::EE_Leaf) ||
1809       (d->position == MAXFLOAT) )
1810  {
1811    //cerr << "This shall not happen" << endl;
1812    SKTBNodeT *n = MakeLeaf(d);
1813    if (!nodeToReturn)
1814      nodeToReturn = n;
1815    // cout << "LL\n";
1816    return nodeToReturn;
1817  }
1818 
1819  if (verbose) {
1820#ifdef _DEBUG
1821
1822//#define _KTBPRINTCUTS
1823#ifndef _KTBPRINTCUTS
1824    if (GetDepth() == 0)
1825#else
1826    if (_printCuts)
1827#endif
1828    {
1829      cout << "position: " << d->position << ", axis: "
1830             << (int)axis << ", depth=" << GetDepth() << ", box:" << endl;
1831      Describe(d->box, cout, 0);
1832      cout << endl;
1833    }
1834#endif // _DEBUG
1835  }
1836
1837#ifdef _DEBUG
1838  if (d->twoSplits == -1) {
1839    cerr << "Some problem in implementation" << endl;
1840    abort();;
1841  }
1842#endif
1843
1844  // copy the parameters to use for splitting
1845  d->axis = axis;
1846
1847  if (d->twoSplits == 1) {
1848    // create interior node with two successors
1849    // case (1)
1850    // DEBUG << "case 1 \n";
1851    // Make easy cut, using one position
1852    IncDepth();
1853    SKTBNodeT* n = MakeOneCut(d);
1854    if (!nodeToReturn)
1855      nodeToReturn = n;
1856    DecDepth();
1857    return nodeToReturn;
1858  }
1859
1860#if 1
1861  cerr << "Unexpected use in this implementation" << endl;
1862  cerr << "ABORTING " << endl;
1863  abort();;
1864#endif
1865 
1866//  case       case
1867//  (2)   or   (3) structure must be created
1868//   O          O                #
1869//  / \        / \               #
1870// L  RI      LI  R              #
1871//   /  \    / \                 #
1872//   RL RR  LL LR                #
1873
1874#ifdef _DEBUG
1875  if ( (d->box.Min()[axis] >= d->position) ||
1876       (d->position >= d->position2) ||
1877       (d->position2 >= d->box.Max()[axis]) ) {
1878    cerr << " the case 2 or 3 is defined incorrectly\n";
1879    abort();;
1880  }
1881#endif
1882
1883  // To be reworked !!!
1884  abort();;
1885  if (d->twoSplits == 2) {
1886    // -----------------------------------------
1887    // case (2)
1888    // DEBUG << "case 2 \n";
1889
1890    // make L part, linked in DFS order
1891    int cntLeft = 1;
1892    int cntRight = 0;
1893    IncDepth();
1894    SKTBNodeT *n = MakeOneCut(d);
1895    if (!nodeToReturn)
1896      nodeToReturn = n;
1897
1898    // make RI part, linked to the left node explictly
1899    SBBox rbb = d->box; // the bounding box
1900    rbb.Reduce(axis, 1, d->position);
1901    d->box = rbb;
1902    d->pars.reqAxis = axis;
1903    d->pars.reqPosition = d->position2;   
1904    SKTBNodeT *n2 = SubDiv(d);
1905    DecDepth();
1906    //Link the node
1907    SetInteriorNodeLinks(n, 0, n2);
1908    return nodeToReturn;
1909  }
1910
1911  // -----------------------------------------
1912  // case (3) ????????????????? I am not sure if this is correct implementation !!! VH
1913  // DEBUG << "case 3 \n";
1914 
1915  // make LI part
1916  SBBox lbb = d->box; // the bounding box
1917  lbb.Reduce(axis, 0, d->position2);
1918  int cntLeft = 0;
1919  int cntRight = 1;
1920  d->box = lbb;
1921  d->position = d->position2;
1922  d->axis = axis;
1923  // the next subdivision step for R part
1924  d->pars.reqAxis = axis;
1925  d->pars.reqPosition = d->position;
1926 
1927  IncDepth();
1928  SKTBNodeT *nl = SubDiv(d);
1929  if (!nodeToReturn)
1930    nodeToReturn = nl;
1931 
1932  // make R part
1933  SKTBNodeT *n = MakeOneCut(d);
1934  DecDepth();
1935  SetInteriorNodeLinks(nl, 0, n);
1936  // return left node, linked in DFS order to the previous one
1937  return nodeToReturn;
1938}
1939
1940// makes cutting on the given position on given axis
1941// at value position, bb is the bounding box of current node,
1942// count is the number of objects in the node
1943// flags LeftF and RightF declares if to continue down after the recursion
1944SKTBNodeT *
1945CKTBABuildUp::MakeOneCut(SInputData *d)
1946{
1947// for testing accuracy of setting position from the ideal position
1948#ifdef _RANDOMIZE_POSITION
1949  assert(d->algorithmBreakAx == 0); // position based
1950  RandomizePosition(d->position, d->box.Min(d->axis), d->box.Max(d->axis));
1951#endif //_RANDOMIZE_POSITION
1952
1953#ifdef _DEBUG
1954  assert(d->position >= d->box.Min(d->axis));
1955  assert(d->position <= d->box.Max(d->axis));
1956#endif
1957 
1958  SKTBNodeT *nodeToReturn = 0;
1959
1960#ifdef  _KTB_CONSTR_STATS
1961  // surface area of the interior nodes to be subdivided
1962  _sumSurfaceAreaInteriorNodes += d->box.SA2();
1963#endif // _KTB_CONSTR_STATS
1964 
1965  if (splitClip) {
1966    // empty object list to store the pointers to the split objects
1967    ObjectContainer objlist;
1968    cerr << "Not yet handled" << endl;
1969    abort();;
1970  }
1971
1972#if 0
1973  // check if the lists are correctly sorted .. in debug
1974  static int index = 0;
1975#if 1
1976  cout << "MakeOneCut index = " << index << " .. check axis = "
1977        << d->axis << " count = " << d->count
1978        << " depth = " << GetDepth() << endl;
1979  const int indexToFind = 3743;
1980  if (index == indexToFind) {
1981   
1982    cout << "Debug should start " << endl;
1983    cout << "index = indexToFind" << endl;
1984  }
1985#endif
1986  index++; 
1987#endif
1988 
1989#ifdef _DEBUG
1990  // We do not know the number of boundaries
1991  Check3List(d);
1992#endif
1993 
1994//#define _VYPIS
1995#ifdef _VYPIS
1996  DEBUG << "START: Subdivision at pos="
1997        << position << " axis = " << axis << "\n";
1998  PrintOut(axis, 1, list+axis, sec+axis, Limits::Infinity);
1999#endif
2000
2001  // axis .. the axis perpendicular to splitting plane positioned at position
2002  // break active axis = split the list into two pieces
2003  // list[axis] .. start of the first list in axis direction .. input
2004  // second[axis] .. start of the second list in axis direction .. output
2005  // objlist .. the list of objects that are duplicated
2006#ifdef _DEBUG
2007  Check1List(d, d->axis, d->count);
2008  //cout << "Axis = " << axis << endl;
2009#endif
2010
2011  // We have to allocate a new node for right child
2012  SInputData* rightData = AllocNewData(2*d->count);
2013  // Copy the basic data from parent
2014  rightData->CopyBasicData(d);
2015
2016  // The number of objects for left children and right children
2017  int cntL, cntR;
2018
2019//#ifdef _DEBUG
2020#if 1
2021  if (d->cntThickness < 0) {
2022    cerr << "Problem, d->cntThickness = " << d->cntThickness << endl;
2023    abort();;
2024  }
2025#endif
2026
2027  // Correct for cutting off empty space and many left
2028  // boundaries on the splitting plane.
2029  SItemVec *vec = d->GetItemVec(d->axis);
2030
2031  if (d->algorithmBreakAx) {
2032    // Use pointer-based break-ax algorithm
2033    if (d->bestIterator == vec->begin()) {
2034      assert(d->cntThickness == 0);
2035      d->bestIterator--;
2036    }
2037    else {
2038      if (d->count > 1) {
2039        int origThickness = d->cntThickness;
2040        SItemVec::iterator origIterator = d->bestIterator;
2041        SItemVec::iterator it = d->bestIterator;
2042        int caseBreak = 0;
2043        if ((*(it)).IsLeftBoundary()) {
2044          float pos = d->position;
2045          it--;
2046          int ir = 1;
2047          assert(it != vec->begin()-1);
2048          for ( ; true; ) {
2049            if ( (*(it)).pos < pos) {
2050              caseBreak = 0;
2051              d->bestIterator = it;
2052              if (ir > 1)
2053                d->cntThickness -= ir;
2054              break;
2055            }
2056            if (it == vec->begin()) {
2057              // cutting off, #objects on the left = 0
2058              d->cntThickness = 0;
2059              it--;
2060              d->bestIterator = it;
2061              caseBreak = 1;
2062              break;
2063            }
2064            if ((*(it)).IsRightBoundary()) {
2065              d->bestIterator = it;
2066              if (ir > 1)
2067                d->cntThickness -= ir;
2068              caseBreak = 2;
2069              break;
2070            }
2071            it--;
2072            ir++;
2073          } // for
2074#ifdef _DEBUG
2075          if (d->cntThickness < 0) {
2076            cerr << "Problem, d->cntThickness = " << d->cntThickness << endl;
2077            cerr << " origThickness = " << origThickness << " Depth = "
2078                  << GetDepth() << endl;
2079            for (SItemVec::iterator iq = it; iq != origIterator; iq++) {
2080              cout << "  pos = " << (*iq).pos;
2081              if ((*iq).IsLeftBoundary())
2082                cout << "  L  ";
2083              else
2084                cout << "  R  ";
2085              cout << (*iq).obj << endl;         
2086            } // for
2087            cout << "AAAA" << endl;
2088            abort();;
2089        }
2090#endif
2091        } // if left boundary
2092      } // if d->count > 1
2093    } // if d->vec is not the beginning of the list
2094  } // pointer based break-ax algorithm
2095   
2096#if 0
2097  DEBUG << "MakeOneCut cntL = " << cntL << " cntR = " << cntR << endl;
2098#endif
2099
2100  // Here is the breaking the arrays into two arrays, compute the number
2101  // of objects on the left and on the right of the splitting plane
2102  int dCountTagging = 0;
2103  // Use either
2104  if (d->algorithmBreakAx == 0)
2105    BreakAxPosition(d, d->axis, rightData, cntL, cntR);
2106  else
2107    BreakAx(d, d->axis, rightData, cntL, cntR);
2108   
2109#ifdef _DEBUG
2110  // the number of objects on the left and on the right
2111  int cntL_B = cntL;
2112  int cntR_B = cntR;
2113#if 0
2114  cout << "B - Count*2 = " << d->count*2
2115       << " countSum = " << cntL_B + cntR_B
2116       << " CntL_B = " << cntL_B
2117       << " cntR_B = " << cntR_B << endl;
2118#endif
2119  //Check1List(alist, axis, cntL_B, true, position);
2120  Check1List(d, d->axis, cntL_B);
2121  Check1List(rightData, d->axis, cntR_B);
2122  // cout << "AxisNext = " << oax0 << endl;
2123#endif
2124 
2125  // DEBUG << "THE CHECK at DEPTH=" << GetDepth() << endl;
2126 
2127#ifdef _VYPIS
2128  DEBUG << "After Break at pos="
2129        << position << " axis = " << axis << "\n";
2130  PrintOut(d, axis, d->position);
2131#endif
2132
2133#ifdef _VYPIS
2134  DEBUG << "BEFORE oax0\n";
2135  PrintOut(d, oax0, Limits::Infinity);
2136#endif
2137
2138  // --------------------------------------------------
2139  // break other two axes
2140  CKTBAxes::Axes oax0 = oaxes[d->axis][0];
2141#ifdef _DEBUG
2142  Check1List(d, oax0, d->count);
2143#endif
2144 
2145  SItemVec *vecc = d->GetItemVec(oax0);
2146  if (vecc->size() != 2 * d->count) {
2147    cout << "AAAA 1" << endl;
2148  }
2149 
2150  // Here break the list in the next axis
2151#ifdef _USE_OPTIMIZE_DIVIDE_AX
2152  DivideAx_I_opt(d, oax0, rightData, cntL, cntR);
2153#else
2154  DivideAx_I(d, oax0, rightData, cntL, cntR);
2155#endif
2156
2157  if (vecc->size() != 2 * cntL) {
2158    cout << "AAAA 3" << endl;
2159  }
2160  if (rightData->GetItemVec(oax0)->size() != 2 * cntR) {
2161    cout << "AAAA 4" << endl;
2162    cout << "cntR * 2 = " << 2 * cntR << endl;
2163    cout << "vecSize = " << rightData->GetItemVec(oax0)->size() << endl;
2164  }
2165 
2166#ifdef _DEBUG
2167  // the number of objects on the left and on the right
2168  int cntL_I = cntL;
2169  int cntR_I = cntR;
2170#if 0
2171  cout << "I - Count*2 = " << count*2
2172       << " countSum = " << cntL_I + cntR_I
2173       << " CntL_B = " << cntL_I
2174       << " cntR_B = " << cntR_I << endl;
2175#endif
2176 
2177  if (cntL_I != cntL_B) {
2178    cerr << "FirstList - Problem with DivideAx_I, cntL_B = " << cntL_B
2179          << " and cntL_I = " << cntL_I << endl;
2180  }
2181  if (cntR_I != cntR_B) {
2182    cerr << "SecondList = Problem with DivideAx_I, cntR_B = " << cntR_B
2183          << " and cntR_I = " << cntR_I << endl;
2184  }
2185 
2186  Check1List(d, oax0, cntL_B);
2187  Check1List(rightData, oax0, cntR_B);
2188
2189  //cout << "AxisNextNext = " << oax1 << endl;
2190
2191#endif
2192 
2193#ifdef _VYPIS
2194  DEBUG << "AFTER DivideAx oax0\n";
2195  PrintOut(oax0, 3, alist + oax0, sec+oax0, 0);
2196#endif
2197
2198#ifdef _VYPIS
2199  DEBUG << "BEFORE oax1\n";
2200  PrintOut(oax1, 1, alist + oax1, sec+oax1, Limits::Infinity);
2201#endif
2202
2203  // Her break the list in the next next axis
2204  CKTBAxes::Axes oax1 = oaxes[d->axis][1];
2205
2206  SItemVec *vec2 = d->GetItemVec(oax1);
2207  if (vec2->size() != 2 * d->count) {
2208    cout << "BBBB 1" << endl;
2209  }
2210 
2211#ifdef _USE_OPTIMIZE_DIVIDE_AX
2212  DivideAx_II_opt(d, oax1, rightData, cntL, cntR);
2213#else
2214  DivideAx_II(d, oax1, rightData, cntL, cntR);
2215#endif
2216
2217  if (vec2->size() != 2 * cntL) {
2218    cout << "BBBB 3" << endl;
2219    cout << "cntL * 2 = " << 2 * cntL << endl;
2220    cout << "vec2size = " << vec2->size() << endl;
2221  }
2222  if (rightData->GetItemVec(oax1)->size() != 2 * cntR) {
2223    cout << "BBBB 4" << endl;
2224    cout << "cntR * 2 = " << 2 * cntR << endl;
2225    cout << "vecSize = " << rightData->GetItemVec(oax1)->size() << endl;
2226  }
2227 
2228#ifdef _DEBUG
2229  int cntL_II = cntL;
2230  int cntR_II = cntR;
2231#if 0
2232  cout << "II - Count*2 = " << count*2
2233       << " countSum = " << cntL_II + cntR_II
2234       << " CntL_B = " << cntL_II
2235       << " cntR_B = " << cntR_II << endl;
2236#endif
2237  if (cntL_II != cntL_B) {
2238    cerr << "FirstList - Problem with DivideAx_II, cntL_B = " << cntL_B
2239          << " and cntL_I = " << cntL_I
2240          << " and cntL_II = " << cntL_II << endl;
2241  }
2242  if (cntR_II != cntR_B) {
2243    cerr << "Second List - Problem with DivideAx_II, cntR_B = " << cntR_B
2244          << " and cntR_I = " << cntR_I
2245          << " and cntR_II = " << cntR_II << endl;
2246  }
2247 
2248  Check1List(d, oax1, cntL_B);
2249  Check1List(rightData, oax1, cntR_B);
2250
2251  if ( (cntL_I != cntL_II) ||
2252       (cntR_I != cntR_II) ) {
2253    cout << "Problem with dividing axis II: cntL_I = " << cntL_I
2254         << " cntL_II = " << cntL_II
2255         << " cntR_I = " << cntR_I
2256         << " cntR_II = " << cntR_II << endl;     
2257  }
2258#endif
2259 
2260#ifdef _VYPIS
2261  DEBUG << "AFTER DivideAx oax1\n";
2262  PrintOut(oax1, 3, alist + oax1, sec+oax1, position);
2263#endif
2264 
2265  // Allocate the representation of the interior node
2266  SKTBNodeT *node = 0;
2267
2268  if (makeMinBoxHere) {
2269    // -------------------------------------------------------
2270    // interior node with the box
2271    node = AllocInteriorNodeWithBox(// interior node data
2272                                    d->axis, d->position, cntL, cntR,
2273                                    // box data
2274                                    d->box, d->lastMinBoxNode, d->lastDepthForMinBoxes);
2275
2276    //cout << "depth = " << GetDepth() << " node = " << node
2277    //   << " lastMinBoxNode = " << d->lastMinBoxNode << endl;
2278   
2279    d->lastMinBoxNode = node; // we have to remember the last node
2280    rightData->lastMinBoxNode = node;
2281    makeMinBoxHere = false; // reset for the next time
2282  }
2283  else
2284    // normal interior node (=without box)
2285    node = AllocInteriorNode(d->axis, d->position, cntL, cntR);
2286 
2287#ifdef _DEBUG
2288  if (!node) {
2289    cerr << "Allocation of Interior Node failed.\n";
2290    abort();;
2291  }
2292#endif // _DEBUG
2293  // set the node to be returned
2294  if (!nodeToReturn)
2295    nodeToReturn = nodeToLink;
2296  assert(node);
2297  assert(nodeToReturn);
2298
2299#if 0 
2300  // If we make split clipping, then we want to reduce the boxes for
2301  // the objects straddling the splitting plane now.
2302  if ((splitClip) && (objlist.ListCount() > 0 )) {
2303    // if the objects are split by splitting plane
2304    //DEBUG << "Reduces box ax=" << axis << " position="
2305    // << position << "#split_objs=" << objlist.ListCount() << "\n" << flush;
2306    ReduceBBoxes(d, d->axis, rightData, d->position);
2307  }
2308#endif
2309 
2310  // ---------------------------------------------
2311  // initialize the left bounding box 'bb'
2312  d->box.Reduce(d->axis, 0, d->position);
2313
2314  // initialize the right bounding box 'rbb'
2315  rightData->box.Reduce(d->axis, 1, d->position);
2316
2317  // Set correct number of objects
2318  d->count = cntL;
2319  rightData->count = cntR;
2320 
2321  // The children nodes to be created
2322  SKTBNodeT *nodeL=0, *nodeR=0;
2323
2324  // Recurse to the left child
2325  if (d->makeSubdivisionLeft) {
2326    // subdivide branches
2327    ESubdivMode modeLeft = d->modeSubDiv ;
2328    if ( (cutEmptySpace) &&
2329         (modeLeft != EE_SUBDIVCUTEMPTY) &&
2330         (d->count <= maxListLength) ) {
2331      // setting initial depth if empty space cutting
2332      startEmptyCutDepth = GetDepth();
2333      d->modeSubDiv = EE_SUBDIVCUTEMPTY;
2334    }
2335    // DEBUG << "Left cnt= " << (cnt[0] >> 1) << " depth= "
2336    //       << GetDepth() << "\n";
2337    nodeL = SubDiv(d);
2338  }
2339 
2340  // Recurse to the right child
2341  if (rightData->makeSubdivisionRight) {
2342    ESubdivMode modeRight = d->modeSubDiv ;
2343    if ( (cutEmptySpace) &&
2344         (modeRight != EE_SUBDIVCUTEMPTY) &&
2345         (rightData->count <= maxListLength) ) {
2346      // setting initial depth if empty space cutting
2347      startEmptyCutDepth = GetDepth();
2348      rightData->modeSubDiv = EE_SUBDIVCUTEMPTY;
2349    }
2350    // DEBUG << "Right cnt= " << (cnt[1] >> 1) << " depth= "
2351    // << GetDepth() << "\n";
2352    nodeR = SubDiv(rightData);
2353
2354  }
2355
2356  // Free the data to be reused further on
2357  FreeLastData();
2358
2359  // Set the node pointers to the children for created node
2360  SetInteriorNodeLinks(node, nodeL, nodeR);
2361
2362  return nodeToReturn; // return the node
2363}
2364
2365// returns a box enclosing all the objects in the node
2366void
2367CKTBABuildUp::GetTightBox(const SInputData &i, SBBox &tbox)
2368{
2369  // Index to the last boundary
2370  int cnt = 2*i.count-1; 
2371
2372  tbox.MinX() = (*(i.xvec))[0].pos;
2373  tbox.MinY() = (*(i.yvec))[0].pos;
2374  tbox.MinZ() = (*(i.zvec))[0].pos;
2375  tbox.MaxX() = (*(i.xvec))[cnt].pos;
2376  tbox.MaxY() = (*(i.yvec))[cnt].pos;
2377  tbox.MaxZ() = (*(i.zvec))[cnt].pos;
2378
2379  assert(tbox.MinX() >= i.box.MinX());
2380  assert(tbox.MaxX() <= i.box.MaxX());
2381  assert(tbox.MinY() >= i.box.MinY());
2382  assert(tbox.MaxY() <= i.box.MaxY());
2383  assert(tbox.MinZ() >= i.box.MinZ());
2384  assert(tbox.MaxZ() <= i.box.MaxZ());
2385}
2386
2387
2388// returns a box enclosing all the objects in the node
2389int
2390CKTBABuildUp::GetEBox(const SInputData &i, SBBox &tbox)
2391{
2392  // LEFT BOX
2393  // initialization
2394
2395  // Index to the last boundary
2396  int cnt = 2*i.count-1; 
2397  int changed = 0; // number of planes, that are changed by bounding box
2398
2399  float xMin = (*(i.xvec))[0].pos;
2400  changed = (xMin > i.box.Min(0)) ? changed+1 : changed;   
2401  float xMax = (*(i.xvec))[cnt].pos;
2402  changed = (xMax < i.box.Max(0)) ? changed+1 : changed;   
2403
2404  float yMin = (*(i.yvec))[0].pos;
2405  changed = (yMin > i.box.Min(1)) ? changed+1 : changed;   
2406  float yMax = (*(i.yvec))[cnt].pos;
2407  changed = (yMax < i.box.Max(1)) ? changed+1 : changed;   
2408
2409  float zMin = (*(i.zvec))[0].pos;
2410  changed = (zMin > i.box.Min(2)) ? changed+1 : changed;   
2411  float zMax = (*(i.zvec))[cnt].pos;
2412  changed = (zMax < i.box.Max(2)) ? changed+1 : changed;   
2413
2414  // Set the bounding (reduced) box
2415  tbox.Min() = Vector3(xMin, yMin, zMin);
2416  tbox.Max() = Vector3(xMax, yMax, zMax);
2417
2418  // Return the number of bounding splitting planes with respect to original box
2419  return changed;
2420}
2421
2422
2423// This is only for debugging purposes, to be removed!
2424static
2425CKTBABuildUp::SSolid* objtf = (CKTBABuildUp::SSolid*)0x0;
2426
2427// removes the objects specified in "tagobjlist" from boundary "list"
2428// it supposes the "tagobjlist" is ordered in the left boundary order
2429// in CKTBAxes::EE_X_axis.
2430void
2431CKTBABuildUp::RemoveObjects(SItemVec *vec, int cntObjects)
2432{
2433  assert(vec);
2434  assert(cntObjects);
2435 
2436  // number of removed left boundaries
2437  int cnt = 0;
2438
2439  SItemVec::iterator curr;
2440#ifdef _DEBUG
2441  int cnt2 = 0;
2442  for (curr = vec->begin(); curr != vec->end(); curr++)
2443  {
2444    if (curr->obj->ToBeRemoved()) {
2445      //cout << cnt2 << "  " << curr - vec->begin() << endl;
2446      cnt2++;
2447    }
2448    //if (curr->obj == objtf)
2449    //  cout << "Index found cnt2 = " << cnt2 << endl;
2450  } // for
2451  if (cnt2 != 2*cntObjects) {
2452    cerr << "Some problem with data marking" << endl;
2453    abort();;
2454  }
2455#endif
2456
2457  // Find the first boundary
2458  for (curr = vec->begin(); curr != vec->end(); curr++)
2459  {
2460    //if (curr->obj == objtf)
2461    //cout << "2 Index found cnt2 = " << cnt2 << endl;
2462
2463    if (curr->obj->ToBeRemoved()) {
2464      cnt++; // the number of removed boundaries
2465      assert(curr->IsLeftBoundary());
2466      //cout << cnt << "  objToRemove=" << curr->obj << endl;
2467      break;
2468    }
2469  }
2470  // the first element to be overwritten
2471  SItemVec::iterator it = curr;
2472  curr++;
2473
2474  for ( ; curr != vec->end(); curr++) {
2475    //if (curr->obj == objtf) {
2476    //  cout << "3 Index found cnt2 = " << cnt2 << endl;
2477
2478    if (curr->obj->ToBeRemoved()) {
2479      cnt++; // the number of removed boundaries
2480      //cout << cnt << "  " << curr->obj << endl;
2481      if (cnt == cntObjects*2) {
2482        curr++;
2483        break;
2484      }
2485    }
2486    else {
2487      // copy the element to the right place
2488      *it = *curr;
2489      it++;
2490    }
2491  } // for
2492
2493#ifdef _DEBUG
2494  if (cnt != 2*cntObjects) {
2495    cerr << "Problem with removing objects implementation" << endl;
2496    abort();;
2497  }
2498#endif
2499
2500  // copy the rest of the list to the right place
2501  for ( ; curr != vec->end(); curr++, it++)
2502  {
2503    //if (curr->obj == objtf) {
2504    //  cout << "4 Index found cnt2 = " << cnt2 << endl;
2505
2506    // copy the element to the right place
2507    *it = *curr;
2508  }
2509     
2510  // resize the vector correctly
2511  int oldSize = vec->size();
2512  assert(oldSize >= cnt);
2513  if (oldSize == cnt)
2514    vec->erase(vec->begin(), vec->end());
2515  else
2516    vec->resize(oldSize - cnt);
2517  assert(vec->size() % 2 == 0);
2518  assert(vec->size() == oldSize - cnt);
2519
2520#ifdef _DEBUG   
2521  for (curr = vec->begin(); curr != vec->end(); curr++)
2522  {
2523    if (curr->obj->ToBeRemoved()) {
2524      cerr << "Implementation bug" << endl;
2525      abort();;
2526    } // if left boundary
2527  } // for
2528#endif
2529
2530  //cout << "size = " << vec->size() << endl;
2531
2532  return;
2533}
2534
2535// removes the objects specified in "tagobjlist" from boundary "list"
2536// it supposes the "tagobjlist" is ordered in the left boundary order
2537// in CKTBAxes::EE_X_axis.
2538void
2539CKTBABuildUp::RemoveObjectsReset(SItemVec *vec, int cntObjects)
2540{
2541  assert(vec);
2542  assert(cntObjects);
2543 
2544  // number of removed left boundaries
2545  int cnt = 0;
2546  SItemVec::iterator curr;
2547
2548#ifdef _DEBUG
2549  int cnt2 = 0;
2550  for (curr = vec->begin(); curr != vec->end(); curr++)
2551  {
2552    if (curr->obj->ToBeRemoved()) {
2553      //cout << cnt2 << "  " << curr - vec->begin() << endl;
2554      cnt2++;
2555    }
2556  } // for
2557  if (cnt2 != 2*cntObjects) {
2558    cerr << "Some problem with data marking" << endl;
2559    abort();;
2560  }
2561#endif
2562
2563  // Find the first boundary
2564  for (curr = vec->begin(); curr != vec->end(); curr++)
2565  {
2566    if (curr->obj->ToBeRemoved()) {
2567      cnt++; // the number of removed boundaries
2568      assert(curr->IsLeftBoundary());
2569      curr++;
2570      break;
2571    }
2572  }
2573  SItemVec::iterator it = curr;
2574  it--; // the first element to be overwritten
2575
2576  for ( ; curr != vec->end(); curr++) {
2577    if (curr->obj->ToBeRemoved()) {
2578      cnt++; // the number of removed boundaries
2579      if (curr->IsRightBoundary()) {
2580        curr->obj->ResetFlags();
2581        assert(curr->obj->flags == 0);
2582      }
2583      if (cnt == cntObjects*2) {
2584        curr++;
2585        break;
2586      }
2587    }
2588    else {
2589      // copy the element to the right place
2590      *it = *curr;
2591      it++;
2592    }
2593  } // for
2594 
2595#ifdef _DEBUG
2596  if (cnt != 2*cntObjects) {
2597    cerr << "Problem with removing objects implementation" << endl;
2598    abort();;
2599  }
2600#endif
2601
2602  // copy the rest of the list to the right place
2603  for ( ; curr != vec->end(); curr++, it++)
2604  {
2605    // copy the element to the right place
2606    *it = *curr;
2607  }
2608     
2609  // resize the vector correctly
2610  int oldSize = vec->size();
2611  assert(oldSize >= cnt);
2612  if (oldSize == cnt)
2613    vec->erase(vec->begin(), vec->end());
2614  else
2615    vec->resize(oldSize - cnt);
2616  assert(vec->size() % 2 == 0);
2617  assert(vec->size() == oldSize - cnt);
2618
2619#ifdef _DEBUG   
2620  for (curr = vec->begin(); curr != vec->end(); curr++)
2621  {
2622    if (curr->obj->ToBeRemoved()) {
2623      cerr << "Implementation bug" << endl;
2624      abort();;
2625    } // if left boundary
2626  } // for
2627#endif
2628
2629  //cout << "size = " << vec->size() << endl;
2630 
2631  return;
2632}
2633
2634
2635// make the full leaf from current node
2636SKTBNodeT*
2637CKTBABuildUp::MakeLeaf(SInputData *d)
2638{
2639#ifdef  _KTBPRINTCUTS
2640  cout << " Creating leaf, depth = " << GetDepth()
2641       << " cntObjs = " << count << endl;
2642  //exit(3);
2643#endif
2644 
2645#ifdef _KTB_CONSTR_STATS
2646  // update the statistics
2647  float sa2bb;
2648  _sumSurfaceAreaLeaves += (sa2bb = d->box.SA2());
2649  _sumSurfaceAreaMULcntLeaves += ((float)d->count * sa2bb);
2650#endif // _KTB_CONSTR_STATS
2651 
2652  // ------------------------------------------------------
2653  // The version with allocating empty leaves. This makes
2654  // sense since in DFS order it works correctly
2655  SKTBNodeT *node = AllocLeaf(d->count);
2656
2657  if (d->count > 0) {
2658    // copy the list of objects
2659    ObjectContainer *objlist = NULL;
2660    objlist = new ObjectContainer;
2661
2662    SItemVec *vec = d->GetItemVec(0);
2663    SItemVec::iterator curr;   
2664    for (curr = vec->begin(); curr != vec->end(); curr++) {
2665      if ( (*curr).IsLeftBoundary()) // skip all right boundaries
2666        objlist->push_back( (*curr).obj->obj);
2667    }
2668
2669    assert(objlist->size() != 0);
2670    // Do not delete the object list, it is used as allocated above
2671    SetFullLeaf(node, objlist);
2672  }
2673
2674  // We assume that the allocation of the single leaf cannot
2675  // fail, since this one shall be the last in the sequence.
2676  assert(node == nodeToLink);
2677 
2678  // Return the node to be used in the linking
2679  return nodeToLink;
2680}
2681
2682// split the list along the required axis
2683// first is the begin of the SItem list, axis is the axis, where
2684// to split, val is the splitting value in integer notation.
2685// The output of the function are two lists, first and second
2686// Optionally, objlist store the objects that straddle the splitting plane
2687void
2688CKTBABuildUp::BreakAx(SInputData *d, int axis,
2689                      SInputData *rightData,
2690                      int &cntLout, int &cntRout)
2691{
2692  // the number of object boundaries for left array
2693  int cntLeft = (d->bestIterator) - d->GetItemVec(axis)->begin() + 1;
2694  assert(cntLeft >= 0);
2695  assert(cntLeft <= 2*d->count); 
2696  // the number of object boundaries for right array
2697  int cntRight = d->count*2 - cntLeft;
2698  assert(cntRight >= 0);
2699  assert(cntRight <= 2*d->count);
2700  assert(d->cntThickness >= 0);
2701  cntLeft += d->cntThickness; // duplicated boundaries for left array
2702  cntRight += d->cntThickness; // duplicated boundaries for right array
2703  assert(cntLeft + cntRight == (d->count + d->cntThickness) * 2);
2704
2705#ifdef _DEBUG
2706  if ( (cntLeft % 2 != 0) ||
2707       (cntRight % 2 != 0) ) {
2708    int i = 0;
2709    SItemVec *vec = d->GetItemVec(axis);
2710    cout << "d->count = " << d->count ;
2711    cout << "cntLeft = " << cntLeft << " cntRight" << cntRight << endl;
2712    cout << "bestIterator = " << (int)(d->bestIterator - vec->begin());
2713    cout << " cntThickness= " << d->cntThickness << " position=" << d->position << endl;
2714    for (SItemVec::iterator itt = vec->begin(); itt != vec->end(); itt++, i++)
2715    {
2716      cout << i << " = ";
2717      if ((*itt).IsLeftBoundary())
2718        cout << "L"; else cout << "R";
2719      cout << " obj = " << (*itt).obj << " pos = " << (*itt).pos << endl;
2720    }
2721    abort();;
2722  }
2723#endif
2724 
2725  // A special case, right list is empty
2726  if (cntRight == 0) {
2727    assert(d->cntThickness == 0);
2728    rightData->Alloc(0);
2729    rightData->count = 0;
2730    // The final setting of number of objects
2731    cntLout = cntLeft / 2;
2732    cntRout = cntRight / 2;   
2733    return; // nothing to do
2734  }
2735
2736  // We have some problem that is difficult to detect, when tagging objects
2737  // are allowed, the flags are not always reset to zero in some cases.
2738  // VH 2/1/2006. The flags are always left boundaries (=1).
2739  // Therefore here it is a hack to assure that the flags are definitely zero.
2740  if (resetFlagsForBreakAx) {
2741    SItemVec *vecd = d->GetItemVec(axis);
2742    for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
2743      (*itt).obj->ResetFlags();
2744    } // for
2745  } // if
2746 
2747#ifdef _DEBUG
2748  SItemVec *vecd = d->GetItemVec(axis);
2749  for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
2750    if ((*itt).obj->Flags() != 0) {
2751      cout << "Init ERROR 1 in original list, flags = "
2752           << (*itt).obj->Flags() << " obj = " << (void*)((*itt).obj) << endl;
2753    }
2754  } // for
2755#endif // _DEBUG 
2756
2757  assert(cntRight > 0);
2758 
2759  // Allocate the appropriate size for right array, for all 3 axes
2760  rightData->Alloc(cntRight);
2761 
2762  // First find create all left boundaries for right array
2763  SItem boundary; // create left boundary
2764  SItemVec *vecRight = rightData->GetItemVec(axis);
2765
2766#define _USE_PUSHBACK
2767#ifdef _USE_PUSHBACK
2768  vecRight->resize(0);
2769#else
2770  vecRight->resize(cntRight);
2771#endif 
2772  if (d->cntThickness) {
2773    boundary.pos = d->position;
2774    boundary.axis = axis;
2775    boundary.SetLeftBoundary();
2776#ifdef _USE_PUSHBACK
2777    int i;
2778    for (i = 0; i < d->cntThickness; i++) {
2779      vecRight->push_back(boundary); // Note that the objects pointed to are not set !!!
2780    } // for i;
2781#else
2782    for (int i = 0; i < d->cntThickness; i++) {
2783      (*vecRight)[i] = boundary; // Note that the objects pointed to are not set !!!
2784    } // for i;
2785#endif
2786  }
2787  const SItemVec::iterator itendLeft = d->bestIterator + 1;
2788  SItemVec *vec = d->GetItemVec(axis);
2789  SItemVec::iterator it;
2790  //int cntSetLeft = 0;
2791  for (it = vec->begin(); it != itendLeft; it++) {
2792    // Mark all the items to be in the first list
2793    (*it).obj->SetInFirstList();
2794    //cntSetLeft++;
2795  } // for
2796  // The boundaries that belong definitely to the second list
2797  //int cntSetRight = 0;
2798#ifdef _USE_PUSHBACK
2799  const SItemVec::iterator itend = vec->end();
2800  for (; it != itend; it++) {
2801    // Mark all the items to be in the second list
2802    (*it).obj->SetInSecondList();
2803    // and copy the item to the right array
2804    vecRight->push_back(*it);
2805    //cntSetRight++;
2806  } // for
2807#else
2808  const int imax = vec->end() - (d->bestIterator + 1);
2809  int ir = d->cntThickness;
2810  for (int i = 0; i < imax; i++) {
2811    (*it).obj->SetInSecondList();
2812    // and copy the item to the right array
2813    (*vecRight)[ir] = (*it);
2814    ir++; it++;
2815  }
2816  assert(ir == cntRight);
2817#endif
2818 
2819#ifdef _DEBUG
2820  if (vecRight->size() != cntRight) {
2821    cerr << "Implementation problem vecRight->size() = "
2822          << vecRight->size() << "  cntRight = " << cntRight
2823          << " d->cntThickness = " << d->cntThickness << endl;
2824    abort();;
2825  }
2826#endif 
2827  // Only check, if right boundary is set correctly
2828  assert(vecRight->size() == cntRight);
2829 
2830  if (d->cntThickness) {
2831    // Now go only through the items on the left side of the splitting plane
2832    SItemVec::iterator itR = vecRight->begin(); // to set left boundaries in right array
2833    SItemVec::iterator itL = itendLeft; // to set new right boundaries in left array
2834    boundary.SetRightBoundary();
2835#ifdef _USE_PUSHBACK
2836    for (it = vec->begin(); it != itendLeft; it++) {
2837      if ((*it).obj->InBothLists()) {
2838        // set new right boundaries in left array
2839        (*itL) = boundary;
2840        (*itL).obj = (*it).obj; // set the right object
2841        itL++;
2842        // and set new left boundary in right array, only object is set
2843        (*itR).obj = (*it).obj;
2844        itR++;
2845      }
2846    } // for
2847#else
2848    int imax = itendLeft - vec->begin();
2849    int i;
2850    for (i = 0; i < imax ; i++) {
2851      if ((*vec)[i].obj->InBothLists()) {
2852        // set new right boundaries in left array
2853        (*itL) = boundary;
2854        (*itL).obj = (*vec)[i].obj; // set the right object
2855        itL++;
2856        // and set new left boundary in right array, only object is set
2857        (*itR).obj = (*vec)[i].obj;
2858        itR++;
2859      }
2860    } // for
2861#endif
2862#ifdef _DEBUG
2863    int sizeL = itL - vec->begin();
2864    assert(sizeL == cntLeft);
2865    int sizeR = itR - vecRight->begin();
2866    assert(sizeR == d->cntThickness);
2867#endif
2868  }
2869#ifdef _DEBUG
2870  else {
2871    vec->resize(cntLeft);   
2872    assert(d->cntThickness == 0);
2873    int cntLL = 0;
2874    int cntRR = 0;
2875    for (it = vec->begin(); it != vec->end(); it++) {
2876      if ((*it).obj->InBothLists()) {
2877        cntLL++;
2878        cout << "ERROR in left list " << cntLL << endl;
2879      }
2880    } // for   
2881    for (it = vecRight->begin(); it != vecRight->end(); it++) {
2882      if ((*it).obj->InBothLists()) {
2883        cntRR++;
2884        cout << "ERROR in right list " << cntRR << endl;
2885      }
2886    } // for   
2887  } 
2888#endif // _DEBUG 
2889 
2890  // The boundary on the left, unused items are not used
2891  vec->resize(cntLeft);
2892
2893  // The final setting of size
2894  cntLout = cntLeft / 2;
2895  cntRout = cntRight / 2;
2896
2897  return;
2898}
2899
2900
2901void
2902CKTBABuildUp::BreakAxPosition(SInputData *d, int axis,
2903                              SInputData *rightData,
2904                              int &cntLout, int &cntRout)
2905{
2906  // the number of object boundaries for left array is not known
2907  int cntLeft = 0;
2908  int cntRight = 0;
2909  d->cntThickness = 0;
2910  const float positionToDecide = d->position;
2911
2912  // We have some problem that is difficult to detect, when tagging objects
2913  // are allowed, the flags are not always reset to zero in some cases.
2914  // VH 2/1/2006. The flags are always left boundaries (=1).
2915  // Therefore here it is a hack to assure that the flags are definitely zero.
2916  if (resetFlagsForBreakAx) {
2917    SItemVec *vecd = d->GetItemVec(axis);
2918    for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
2919      (*itt).obj->ResetFlags();
2920    } // for
2921  } // if
2922 
2923#ifdef _DEBUG
2924  SItemVec *vecd = d->GetItemVec(axis);
2925  for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
2926    if ((*itt).obj->Flags() != 0) {
2927      cout << "Init ERROR 2 in original list, flags = "
2928           << (*itt).obj->Flags() << " obj = " << (void*)((*itt).obj) << endl;
2929    }
2930  } // for
2931#endif // _DEBUG 
2932
2933  // Allocate the appropriate size for right array, for all 3 axes
2934  // This is only estimated at this algorithm.
2935  int estimatedSize = (int)(0.6f * (float)d->count);
2936  if (estimatedSize < 10)
2937    estimatedSize = 10;
2938  rightData->Alloc(estimatedSize * 2);
2939 
2940  // First find create all left boundaries for right array
2941  SItemVec *vecRight = rightData->GetItemVec(axis);
2942
2943  // Resize to 0 to start from the beginning
2944  vecRight->resize(0);
2945  SItemVec *vec = d->GetItemVec(axis);
2946  SItemVec::iterator it;
2947  const SItemVec::iterator itEnd = vec->end();
2948  for (it = vec->begin(); it != itEnd; it++) {
2949    const float &pos = (*it).pos; 
2950    if (pos < positionToDecide) {
2951      cntLeft++;
2952      // Mark the item to be in the first list
2953      (*it).obj->SetInFirstList();
2954    }
2955    else {
2956      if (pos > positionToDecide) {
2957        // It belongs to the right list
2958        break;
2959      }
2960      else {
2961        // pos == positionToDecide
2962        if ((*it).IsRightBoundary()) {
2963          // it belongs to the left list
2964          cntLeft++;
2965          (*it).obj->SetInFirstList();
2966        }
2967        else {
2968          // this is the left boundary
2969          // it belongs already to the right list
2970          break;
2971        }
2972      }
2973    }
2974  } // for it
2975
2976  // Remember the position of the iterator, where right
2977  // list has to start
2978  SItemVec::iterator itStart = it;
2979  int cntDups = 0;
2980  SItem boundary;
2981
2982  // create left boundaries for right list
2983  boundary.pos = d->position;
2984  boundary.axis = axis;
2985  boundary.SetLeftBoundary();
2986  for (; it != itEnd; it++) {
2987    assert((*it).pos >= positionToDecide);
2988    // Mark all the items to be in the second list
2989    if ((*it).obj->InFirstList()) {
2990      // and copy the newly created boundary to the right array
2991      boundary.obj = (*it).obj;
2992      vecRight->push_back(boundary);
2993      cntDups++;
2994      cntRight++;
2995    }
2996  } // for
2997
2998  // Now go through the right part of the list again
2999  // and copy the rest of the list to the right list
3000  for (it = itStart; it != itEnd; it++) {
3001    (*it).obj->SetInSecondList();
3002    assert((*it).pos >= positionToDecide);
3003    vecRight->push_back((*it));
3004    cntRight++;
3005  } // for
3006
3007  // in final pass create the new boundaries for the left list
3008  int i;
3009  SItemVec::iterator srcIt = vecRight->begin();
3010  for (it = itStart, i = 0 ; i < cntDups; it++, srcIt++, i++) {
3011    cntLeft++;
3012    // first copy the boundary from the beginning of the right list
3013    *it = *srcIt;
3014    // then set the right boundary in copied item
3015    (*it).SetRightBoundary();
3016    assert( (*srcIt).IsLeftBoundary());
3017  } // for it
3018 
3019  // The boundary on the left, unused items are not used
3020  vec->resize(cntLeft);
3021  assert(cntLeft % 2 == 0);
3022  assert(cntRight % 2 == 0);
3023  assert(cntLeft+cntRight == 2*(d->count+cntDups));
3024 
3025  // The final setting of size
3026  cntLout = cntLeft / 2;
3027  cntRout = cntRight / 2;
3028  // Set correctly also the number of objects
3029  d->cntThickness = cntDups;
3030
3031#ifdef _DEBUG
3032  SItemVec *vecdL = d->GetItemVec(axis);
3033  int cntL = 0; int cntL_LB = 0; int cntL_RB = 0; int cntL_LR = 0;
3034  for (SItemVec::iterator itl = vecdL->begin(); itl != vecdL->end(); itl++) {
3035    if ((*itl).obj->InBothLists())
3036      cntL_LR++;
3037    if ((*itl).IsLeftBoundary()) cntL_LB++;
3038    if ((*itl).IsRightBoundary()) cntL_RB++;
3039    cntL++;
3040  } // for
3041  SItemVec *vecdR = rightData->GetItemVec(axis);
3042  int cntR = 0; int cntR_LB = 0; int cntR_RB = 0; int cntR_LR = 0;
3043  for (SItemVec::iterator itr = vecdR->begin(); itr != vecdR->end(); itr++) {
3044    if ((*itr).obj->InBothLists())
3045      cntR_LR++;
3046    if ((*itr).IsLeftBoundary()) cntR_LB++;
3047    if ((*itr).IsRightBoundary()) cntR_RB++;
3048    cntR++;
3049  } // for
3050  assert(cntL_LR == cntR_LR);
3051  assert(cntL_LB == cntL_RB);
3052  assert(cntR_LB == cntR_RB);
3053#endif
3054 
3055  return;
3056}
3057
3058#if 0
3059// This below was an attempt to do it differently in Vienna, Dec 2007.
3060// VH
3061
3062// ---------------------------------------------------------
3063// split the list along the required axis
3064// first is the begin of the SItem list, axis is the axis, where
3065// to split, val is the splitting value in integer notation.
3066// The output of the function are two lists, first and second
3067// Optionally, objlist store the objects that straddle the splitting plane
3068void
3069CKTBABuildUp::BreakAxPosition(SInputData *d, int axis,
3070                              SInputData *rightData,
3071                              int &cntLout, int &cntRout)
3072{
3073  // the number of object boundaries for left array is not known
3074  int cntLeft = 0;
3075  int cntRight = 0;
3076  d->cntThickness = 0;
3077  const float positionToDecide = d->position;
3078
3079  // We have some problem that is difficult to detect, when tagging objects
3080  // are allowed, the flags are not always reset to zero in some cases.
3081  // VH 2/1/2006. The flags are always left boundaries (=1).
3082  // Therefore here it is a hack to assure that the flags are definitely zero.
3083  if (resetFlagsForBreakAx) {
3084    SItemVec *vecd = d->GetItemVec(axis);
3085    for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
3086      (*itt).obj->ResetFlags();
3087    } // for
3088  } // if
3089 
3090#ifdef _DEBUG
3091  SItemVec *vecd = d->GetItemVec(axis);
3092  for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
3093    if ((*itt).obj->Flags() != 0) {
3094      cout << "Init ERROR 2 in original list, flags = "
3095           << (*itt).obj->Flags() << " obj = " << (void*)((*itt).obj) << endl;
3096    }
3097  } // for
3098#endif // _DEBUG 
3099
3100  // Allocate the appropriate size for right array, for all 3 axes
3101  // This is only estimated at this algorithm.
3102  int estimatedSize = (int)(0.6f * (float)d->count);
3103  if (estimatedSize < 10)
3104    estimatedSize = 10;
3105  assert(estimatedSize > 0);
3106  rightData->Alloc(estimatedSize * 2);
3107 
3108  // First find create all left boundaries for right array
3109  SItemVec *vecRight = rightData->GetItemVec(axis);
3110
3111  // Resize to 0 to start from the beginning
3112  vecRight->resize(0);
3113  SItemVec *vec = d->GetItemVec(axis);
3114  SItemVec::iterator it;
3115  SItemVec::iterator itLeft;
3116  const SItemVec::iterator itEnd = vec->end();
3117  bool foundBoundary = false; 
3118  for (it = vec->begin(); it != itEnd; it++) {
3119    const float &pos = (*it).pos; 
3120    if (pos < positionToDecide) {
3121      cntLeft++;
3122      // Mark the item to be in the first list
3123      (*it).obj->SetInFirstList();
3124    }
3125    else {
3126      if (pos > positionToDecide) {
3127        cntRight++;
3128        (*it).obj->SetInSecondList();
3129        if (!foundBoundary) {
3130          foundBoundary = true;
3131          itLeft = it;
3132        }
3133        // It belongs to the right list
3134      }
3135      else {
3136        // pos == positionToDecide
3137        if ((*it).IsRightBoundary()) {
3138          // it belongs to the left list
3139          cntLeft++;
3140          (*it).obj->SetInFirstList();
3141        }
3142        else {
3143          // this is the left boundary
3144          // it belongs already to the right list
3145          cntRight++;
3146          (*it).obj->SetInSecondList();
3147          if (!foundBoundary) {
3148            foundBoundary = true;
3149            itLeft = it;
3150          }
3151        }
3152      }
3153    }
3154  } // for it
3155
3156  // Remember the position of the iterator, where right
3157  // list has to start
3158  SItem boundary;
3159  int cntDups = 0;
3160  cntLeft = 0;
3161  cntRight = 0;
3162 
3163  // create left boundaries for right list
3164  boundary.pos = d->position;
3165  boundary.axis = axis;
3166  boundary.SetLeftBoundary();
3167  for (it = itLeft; it != itEnd; it++) {
3168    assert((*it).pos >= positionToDecide);
3169    // Mark all the items to be in the second list
3170    if ((*it).obj->InFirstList()) {
3171      // and copy the newly created boundary to the right array
3172      boundary.obj = (*it).obj;
3173      vecRight->push_back(boundary);
3174      cntDups++;
3175      cntRight++;
3176    }
3177  } // for
3178
3179  // Now go through the right part of the list again
3180  // and copy the rest of the list to the right list
3181  for (it = vec->begin(); it != itEnd; it++) {
3182    const float &pos = (*it).pos;
3183    if (pos < positionToDecide) {
3184      // we copy the boundary to the left list
3185      *itLeft = *it;
3186      itLeft++;
3187      cntLeft++;
3188    }   
3189    if (pos > positionToDecide) {
3190      // we copy the boundary to the right list
3191      vecRight->push_back(*it);
3192      cntRight++;
3193    }
3194    else {
3195      if (pos == positionToDecide) {
3196        if ((*it).obj->InBothLists()) {
3197          if ((*it).IsRightBoundary()) {
3198            // right boundary is added to the right as
3199            // we already created left boundaries
3200            vecRight->push_back(*it);
3201            cntRight++;
3202          }
3203          else {
3204            // we copy the left boundary to the left list
3205            *itLeft = *it;
3206            itLeft++;
3207            cntLeft++;
3208          }
3209        }
3210        else {
3211          // This boundary belongs to only a single list
3212          if ((*it).obj->InFirstList()) {
3213            // we copy this boundary  to the left list
3214            *itLeft = *it;
3215            itLeft++;
3216            cntLeft++;
3217          }
3218          else
3219          if ((*it).obj->InSecondList()) {
3220            // we copy this boundary to the right list
3221            vecRight->push_back(*it);
3222            cntRight++;     
3223          }
3224        }
3225      }
3226    }
3227  } // for
3228
3229 
3230  // in final pass create the new boundaries for the left list
3231  int i;
3232  SItemVec::iterator srcIt = vecRight->begin();
3233  for (it = itLeft, i = 0 ; i < cntDups; it++, srcIt++, i++) {
3234    cntLeft++;
3235    // first copy the boundary from the beginning of the right list
3236    assert( (*srcIt).IsLeftBoundary());
3237    *it = *srcIt;
3238    // then set the right boundary in copied item
3239    (*it).SetRightBoundary();
3240    assert( (*it).IsRightBoundary());
3241  } // for it
3242 
3243  // The boundary on the left, unused items are not used
3244  vec->resize(cntLeft);
3245  assert(cntLeft % 2 == 0);
3246  assert(cntRight % 2 == 0);
3247  assert(cntLeft+cntRight == 2*(d->count+cntDups));
3248 
3249  // The final setting of size
3250  cntLout = cntLeft / 2;
3251  cntRout = cntRight / 2;
3252  // Set correctly also the number of objects
3253  d->cntThickness = cntDups;
3254
3255#ifdef _DEBUG
3256  SItemVec *vecdL = d->GetItemVec(axis);
3257  int cntL = 0; int cntL_LB = 0; int cntL_RB = 0; int cntL_LR = 0;
3258  for (SItemVec::iterator itl = vecdL->begin(); itl != vecdL->end(); itl++) {
3259    if ((*itl).obj->InBothLists())
3260      cntL_LR++;
3261    if ((*itl).IsLeftBoundary()) cntL_LB++;
3262    if ((*itl).IsRightBoundary()) cntL_RB++;
3263    cntL++;
3264  } // for
3265  SItemVec *vecdR = rightData->GetItemVec(axis);
3266  int cntR = 0; int cntR_LB = 0; int cntR_RB = 0; int cntR_LR = 0;
3267  for (SItemVec::iterator itr = vecdR->begin(); itr != vecdR->end(); itr++) {
3268    if ((*itr).obj->InBothLists())
3269      cntR_LR++;
3270    if ((*itr).IsLeftBoundary()) cntR_LB++;
3271    if ((*itr).IsRightBoundary()) cntR_RB++;
3272    cntR++;
3273  } // for
3274  assert(vecdL->size() == 2 * cntLout);
3275  assert(vecdR->size() == 2 * cntRout);
3276  assert(vecdL->size() == cntL_LB + cntL_RB);
3277  assert(vecdR->size() == cntR_LB + cntR_RB); 
3278  assert(cntL_LR == cntR_LR);
3279  assert(cntL_LB == cntL_RB);
3280  assert(cntR_LB == cntR_RB);
3281#endif
3282 
3283  return;
3284}
3285#endif
3286
3287// ------------------------------------------------------------------
3288// break the given list of SItem into two parts by reference axis
3289// and reference value, the output is in the first and second SItem
3290void
3291CKTBABuildUp::DivideAx_I(SInputData *d, int axis,
3292                         SInputData *rightData,
3293                         int &cntLout, int &cntRout)
3294{
3295  // First check, if this is not only singular case
3296  if (cntRout == 0) {
3297    assert(d->cntThickness == 0);
3298    return; // nothing to do
3299  }
3300 
3301  // the number of found boundaries on the left and on the right
3302  int cntLeft = 0;
3303  int cntRight = 0;
3304 
3305  // vector for left and right part
3306  SItemVec *vec = d->GetItemVec(axis);
3307#ifdef _DEBUG
3308  if (d->cntThickness == 0) {
3309    SItemVec::iterator it;
3310    for (it = vec->begin(); it != vec->end(); it++) {
3311      if ((*it).obj->InBothLists()) {
3312        cout << "ERROR in left 3 list" << endl;
3313      }
3314      if ((*it).obj->flags > 3) {
3315        cout << "Problem with 4 flags = " << (*it).obj->flags << endl;   
3316      }
3317    } // for   
3318  } 
3319#endif // _DEBUG 
3320 
3321  assert(vec->size() == d->count*2);
3322  SItemVec *vecRight = rightData->GetItemVec(axis);
3323  vecRight->resize(0);
3324 
3325  // Now go through source list
3326  SItemVec::iterator itDest = vec->begin();
3327  SItemVec::iterator itSrc = vec->begin();
3328  // traverse all boundaries, belonging only to the left part
3329  int i = 0;
3330  const int imax = d->count*2;
3331  for (; i != imax; itSrc++, i++) {   
3332    // Check if the item belongs to the second list
3333    if ((*itSrc).obj->InSecondList()) {
3334      break; // we can copy the
3335    }
3336    // Check if the item belongs to the first list
3337    if ((*itSrc).obj->InFirstList()) {
3338      // copy the content to the first list
3339      // Note that *itDest = *itSrc is not necessary, since itDest = itSrc
3340      cntLeft++;
3341    }
3342  } // for
3343
3344  // for newly copied object in the next step
3345  itDest = itSrc;
3346
3347  // the boundaries can belong to both parts
3348  for (; i != imax; itSrc++, i++) {   
3349    // Check if the item belongs to the second list
3350    if ((*itSrc).obj->InSecondList()) {
3351      vecRight->push_back((*itSrc));
3352      cntRight++;
3353    }
3354    // Check if the item belongs to the first list
3355    if ((*itSrc).obj->InFirstList()) {
3356      // copy the content to the first list
3357      *itDest = *itSrc;
3358      itDest++; // and move the iterator
3359      cntLeft++;
3360    }
3361  } // for
3362
3363  assert(cntLeft % 2 == 0);
3364  assert(cntRight % 2 == 0);
3365  assert(cntLeft == 2 * cntLout);
3366  assert(cntRight == 2 * cntRout);
3367
3368  vec->resize(cntLeft); 
3369
3370  return;
3371}
3372
3373// break the given list of SItem into two parts by reference axis
3374// and reference value, the output is in the first and second SItem
3375void
3376CKTBABuildUp::DivideAx_II(SInputData *d, int axis,
3377                          SInputData *rightData,
3378                          int &cntLout, int &cntRout)
3379{
3380  // First check, if this is not only singular case
3381  if (cntRout == 0) {
3382    assert(d->cntThickness == 0);
3383    return; // nothing to do
3384  }
3385 
3386  // the number of found boundaries on the left and on the right
3387  int cntLeft = 0;
3388  int cntRight = 0;
3389 
3390  // vector for left and right part
3391  SItemVec *vec = d->GetItemVec(axis);
3392#ifdef _DEBUG
3393  if (d->cntThickness == 0) {
3394    SItemVec::iterator it;
3395    for (it = vec->begin(); it != vec->end(); it++) {
3396      if ((*it).obj->InBothLists()) {
3397        cout << "ERROR in left list" << endl;
3398      }
3399      if ((*it).obj->flags > 3) {
3400        cout << "Problem with flags = " << (*it).obj->flags << endl;     
3401      }
3402    } // for   
3403  }
3404#endif // _DEBUG 
3405 
3406  assert(vec->size() == d->count*2);
3407  SItemVec *vecRight = rightData->GetItemVec(axis);
3408  vecRight->resize(0);
3409 
3410  // Now go through the source list
3411  SItemVec::iterator itDest = vec->begin();
3412  SItemVec::iterator itSrc = vec->begin();
3413  // traverse all boundaries, belonging only to the left part
3414  int i = 0;
3415  const int imax = d->count*2;
3416  for (; i != imax; itSrc++, i++) {   
3417    // Check if the item belongs to the second list
3418    if ((*itSrc).obj->InSecondList()) {
3419      break; // we can copy the
3420    }
3421    // Check if the item belongs to the first list
3422    if ((*itSrc).obj->InFirstList()) {
3423      // copy the content to the first list
3424      // Note that *itDest = *itSrc is not necessary, since itDest = itSrc
3425      cntLeft++;
3426    }
3427    // Reset flags for the next subdivision step
3428    if ((*itSrc).IsRightBoundary())
3429      (*itSrc).obj->ResetFlags();
3430  } // for
3431
3432  // for newly copied object in the next step
3433  itDest = itSrc;
3434
3435  // !!!!!!!!!!!!!!
3436  // for newly copied object in the next step
3437  itDest = itSrc;
3438
3439  // the boundaries can belong to both parts
3440  for (; i != imax; itSrc++, i++) {   
3441    //if ((*itSrc).obj == objtf) {
3442    //  cout << "DivideAxII - obj found2 , i = " << i << endl;
3443
3444    // Reset flags for the next subdivision step
3445    if ((*itSrc).IsRightBoundary()) {
3446      // This is the right boundary, we have to reset flags
3447      bool inFirstList = (*itSrc).obj->InFirstList();
3448      bool inSecondList = (*itSrc).obj->InSecondList();
3449      (*itSrc).obj->ResetFlags();
3450      if (inSecondList) {
3451        cntRight++;
3452        vecRight->push_back((*itSrc));
3453      }
3454       
3455      // Check if the item belongs to the first list
3456      if (inFirstList) {
3457        cntLeft++;
3458        // copy the content to the first list
3459        *itDest = *itSrc;
3460        itDest++; // and move the iterator
3461      }
3462    }
3463    else { // left boundary
3464      // Check if the item belongs to the second list
3465      if ((*itSrc).obj->InSecondList()) {
3466        cntRight++;
3467        vecRight->push_back((*itSrc));
3468      }
3469      // Check if the item belongs to the first list
3470      if ((*itSrc).obj->InFirstList()) {
3471        cntLeft++;
3472        // copy the content to the first list
3473        *itDest = *itSrc;
3474        itDest++; // and move the iterator
3475      }
3476    } // end of left boundary
3477  } // for
3478 
3479  assert(cntLeft % 2 == 0);
3480  assert(cntRight % 2 == 0);
3481  assert(cntLeft == 2 * cntLout);
3482  assert(cntRight == 2 * cntRout);
3483  vec->resize(cntLeft); 
3484
3485#ifdef _DEBUG
3486  for (SItemVec::iterator itt = vec->begin(); itt != vec->end(); itt++) {
3487    if ((*itt).obj->Flags() != 0) {
3488      cout << "DivideAxII ERROR 1 in final left list, flags = "
3489           << (*itt).obj->Flags() << " obj = " << (void*)((*itt).obj) << endl;
3490      abort();;
3491    }
3492  } // for
3493  for (SItemVec::iterator itt = vecRight->begin(); itt != vecRight->end(); itt++) {
3494    if ((*itt).obj->Flags() != 0) {
3495      cout << "DivideAxII ERROR 1 in final right list, flags = "
3496           << (*itt).obj->Flags() << " obj = " << (void*)((*itt).obj) << endl;
3497      abort();;
3498    }
3499  } // for
3500#endif // _DEBUG 
3501 
3502  return;
3503}
3504
3505
3506// break the given list of SItem into two parts by reference axis
3507// and reference value, the output is in the first and second SItem
3508void
3509CKTBABuildUp::DivideAx_I_opt(SInputData *d, int axis,
3510                             SInputData *rightData,
3511                             int cntLout, int cntRout)
3512{
3513  // First check, if this is not only singular case
3514  if (cntRout == 0) {
3515    assert(d->cntThickness == 0);
3516    return; // nothing to do
3517  }
3518 
3519  // vector for left and right part
3520  SItemVec *vec = d->GetItemVec(axis);
3521
3522#ifdef _DEBUG
3523  if (vec->size() != d->count*2) {
3524    cout << "vec->size() = " << vec->size()
3525         << " d->count*2 = " << d->count*2 << endl;
3526  }
3527
3528  if (d->cntThickness == 0) {
3529    SItemVec::iterator it;
3530    for (it = vec->begin(); it != vec->end(); it++) {
3531      if ((*it).obj->InBothLists()) {
3532        cout << "ERROR in left list" << endl;
3533      }
3534      if ((*it).obj->flags > 3) {
3535        cout << "Problem with flags = " << (*it).obj->flags << endl;     
3536      }
3537    } // for   
3538  } 
3539#endif // _DEBUG 
3540 
3541  assert(vec->size() == d->count*2);
3542  SItemVec *vecRight = rightData->GetItemVec(axis);
3543  vecRight->resize(0);
3544 
3545  // Now go through source list
3546  SItemVec::iterator itDest = vec->begin();
3547  SItemVec::iterator itSrc = vec->begin();
3548#if 1
3549  if (vec->size() != d->count * 2) {
3550    cout << "Something wrong HERE vec->size="
3551         << vec->size() << " d->count*2 = "
3552         << d->count*2 << endl;
3553    abort();
3554  }
3555#endif 
3556
3557  // traverse all boundaries, belonging only to the left part
3558  int i = 0;
3559  const int imax = d->count*2;
3560  for (; i != imax; itSrc++, i++) {   
3561    // Check if the item belongs to the second list
3562    if ((*itSrc).obj->InSecondList()) {
3563      break; // we can copy the
3564    }
3565  } // for
3566
3567  // for newly copied object in the next step
3568  itDest = itSrc;
3569 
3570  // the boundaries can belong to both parts
3571  for (; i != imax; itSrc++, i++) {
3572    SItem &itt = *itSrc;
3573    // Check if the item belongs to the second list   
3574    if (itt.obj->InSecondList()) {
3575      // and first copy the right boundary to a new list
3576      vecRight->push_back((*itSrc));
3577    }
3578    // Check if the item belongs to the first list
3579    if (itt.obj->InFirstList()) {
3580      // copy the content to the first list
3581      *itDest = *itSrc;
3582      itDest++; // and move the iterator
3583    }
3584  } // for
3585 
3586  // Shorten the left array by resizing
3587  vec->resize(cntLout*2);
3588  assert(vecRight->size() == cntRout*2);
3589 
3590  return;
3591}
3592
3593// break the given list of SItem into two parts by reference axis
3594// and reference value, the output is in the first and second SItem
3595void
3596CKTBABuildUp::DivideAx_II_opt(SInputData *d, int axis,
3597                              SInputData *rightData,
3598                              int cntLout, int cntRout)
3599{
3600  // First check, if this is not only singular case
3601  if (cntRout == 0) {
3602    assert(d->cntThickness == 0);
3603    return; // nothing to do
3604  }
3605   
3606  // vector for left and right part
3607  SItemVec *vec = d->GetItemVec(axis);
3608  assert(vec->size() == d->count*2);
3609#ifdef _DEBUG
3610  if (d->cntThickness == 0) {
3611    SItemVec::iterator it;
3612    for (it = vec->begin(); it != vec->end(); it++) {     
3613      if ((*it).obj->InBothLists()) {
3614        cout << "ERROR in left list" << endl;
3615      }
3616      if ((*it).obj->flags > 3) {
3617        cout << "Problem with flags = " << (*it).obj->flags << endl;     
3618      }
3619    } // for   
3620  }
3621  {
3622    SItemVec::iterator ittt;
3623    for (ittt = vec->begin(); ittt != vec->end(); ittt++) {
3624      if ((*ittt).obj->ToBeRemoved()) {
3625        cout << "ERROR - to be removed in the left list" << endl;
3626      }
3627    }
3628  }
3629#endif // _DEBUG 
3630 
3631  assert(vec->size() == d->count*2);
3632  SItemVec *vecRight = rightData->GetItemVec(axis);
3633  vecRight->resize(0);
3634 
3635  // Now go through the source list
3636  SItemVec::iterator itDest = vec->begin();
3637  SItemVec::iterator itSrc = vec->begin();
3638#if 1
3639  if (vec->size() != d->count * 2) {
3640    cout << "Something wrong HERE II vec->size="
3641         << vec->size() << " d->count*2 = "
3642         << d->count*2 << endl;
3643    abort();
3644  }
3645#endif 
3646  // traverse all boundaries, belonging only to the left part
3647  int i = 0;
3648  const int imax = d->count*2;
3649  for (; i != imax; itSrc++, i++) {
3650    //if ((*itSrc).obj == objtf) {
3651    //  cout << "DivideAxII - obj found1 , i = " << i << endl;
3652
3653    // Check if the item belongs to the second list
3654    if ((*itSrc).obj->InSecondList()) {
3655      break; // we can copy the
3656    }
3657    // Reset flags for the next subdivision step
3658    if ((*itSrc).IsRightBoundary())
3659      (*itSrc).obj->ResetFlags();
3660  } // for
3661
3662  // for newly copied object in the next step
3663  itDest = itSrc;
3664
3665  // the boundaries can belong to both parts
3666  for (; i != imax; itSrc++, i++) {   
3667    //if ((*itSrc).obj == objtf) {
3668    //  cout << "DivideAxII - obj found2 , i = " << i << endl;
3669
3670    bool inFirstList = (*itSrc).obj->InFirstList();
3671    bool inSecondList = (*itSrc).obj->InSecondList();
3672
3673    // Reset flags for the next subdivision step
3674    if ((*itSrc).IsRightBoundary()) {
3675      // This is the right boundary, we have to reset flags
3676      if (inSecondList)
3677        vecRight->push_back((*itSrc));
3678      // Check if the item belongs to the first list
3679      if (inFirstList) {
3680        // copy the content to the first list
3681        *itDest = *itSrc;
3682        itDest++; // and move the iterator
3683      }
3684    }
3685    else { // left boundary
3686      // Check if the item belongs to the second list
3687      if (inSecondList) {
3688        vecRight->push_back((*itSrc));
3689      }
3690      // Check if the item belongs to the first list
3691      if (inFirstList) {
3692        // copy the content to the first list
3693        *itDest = *itSrc;
3694        itDest++; // and move the iterator
3695      }
3696    } // end of left boundary
3697  } // for
3698
3699  // Shorten the left array by resizing
3700  vec->resize(cntLout*2);
3701  assert(vecRight->size() == cntRout*2);
3702
3703#if 0
3704#ifdef _DEBUG
3705  for (SItemVec::iterator itt2 = vec->begin(); itt2 != vec->end(); itt2++) {
3706    if ((*itt2).obj->Flags() != 0) {
3707      cout << "DivideAxII ERROR 1 in final left list, flags = "
3708           << (*itt2).obj->Flags() << " obj = "
3709           << (void*)((*itt2).obj) << endl;
3710      abort();;
3711    }
3712  } // for
3713  for (SItemVec::iterator itt3 = vecRight->begin(); itt3 != vecRight->end(); itt3++) {
3714    if ((*itt3).obj->Flags() != 0) {
3715      cout << "DivideAxII ERROR 1 in final right list, flags = "
3716           << (*itt3).obj->Flags() << " obj = "
3717           << (void*)((*itt3).obj) << endl;
3718      abort();;
3719    }
3720  } // for
3721#endif // _DEBUG 
3722#endif
3723 
3724  return;
3725}
3726
3727#ifdef _VYPIS
3728void
3729CKTBABuildUp::PrintOut(CKTBAxes::Axes axis, int stat, SItem **first,
3730                       SItem **second, float /*position*/)
3731{
3732#if 1
3733  int cnt = 0;
3734  SItem *t = *first;
3735  if (stat & 1) {
3736    DEBUG << "----FIRST--- axis=" << (int)axis << "\n";
3737    while (t != NULL) {
3738#define _VYP
3739
3740#ifdef _VYP
3741      float v = t->value[axis];
3742      DEBUGW << "curr=" << t
3743#ifdef _BIDIRLISTS
3744             << " prev=" << t->prev[axis]
3745#endif
3746             << " next="
3747             << t->next[axis] << " val=" << v
3748             << " " << (t->IsRightBoundary() ? 1 : 0);
3749      if (t->IsRightBoundary())
3750        DEBUGW << "   obj= " << t->obj;
3751      else
3752        DEBUGW << "   pair= " << t->pairF;
3753     
3754      DEBUGW << endl;
3755#endif // _VYP
3756      t = t->next[axis];
3757      cnt++;
3758    }
3759   
3760    DEBUG<< "FIRST cnt= " << cnt <<"\n" << flush;
3761  }
3762
3763  if (stat & 2) {
3764    // the right list
3765    t = *second;
3766    cnt = 0;
3767    DEBUG << "----SECOND-----\n";
3768    while (t != NULL) {
3769#ifdef _VYP
3770      float v = t->value[axis];
3771      DEBUGW << "curr=" << t
3772#ifdef _BIDIRLISTS
3773             << " prev=" << t->prev[axis]
3774#endif
3775             << " next="
3776            << t->next[axis] << " val=" << v
3777            << " " << ((t->IsRightBoundary()) ? 1 : 0);
3778      if (t->IsRightBoundary())
3779        DEBUGW << "   obj= " << t->obj;
3780      else
3781        DEBUGW << "   pair= " << t->pairF;
3782      DEBUGW << "\n" << flush;
3783#endif // _VYP
3784      t = t->next[axis];
3785      cnt++;
3786    }
3787    DEBUGW << "SECOND cnt= " << cnt <<"\n" << flush;
3788  }
3789#endif
3790}
3791#endif // _VYPIS
3792
3793// --------------------------------------------------------------------
3794// class CKTBABuildUp::CTestAx
3795
3796// The initialization for the first axis to be tested, this time *****Z axis*******
3797void
3798CKTBABuildUp::SSplitState::InitXaxis(int cnt, const SBBox &boxN)
3799{
3800  // initialize the variables, mainly for surface area estimates
3801  cntAll = cnt;     // the number of all objects in the bounding box
3802  cntLeft = 0;      // the count of bounding boxes on the left
3803  cntRight = cnt;   // the count of bounding boxes on the right
3804  thickness = 0;    // the count of bounding boxes straddling the splitting plane
3805  axis = CKTBAxes::EE_X_axis; // the axis, where the splitting is proposed
3806  box = boxN;        // the box, that is subdivided
3807  //int topAxis = 1;        // axis that is considered in top .. depth
3808  //int frontAxis = 2;      // axis that is considered in front .. height 
3809  // the size of bounding box along the axis to be split
3810  width = box.Max().x - box.Min().x;
3811  // and along two other axes 
3812  topw = box.Max().y - box.Min().y;
3813  frontw = box.Max().z - box.Min().z;
3814  // surface area of the splitting plane .. one face !!
3815  areaSplitPlane = topw * frontw;
3816  // the sum of length of the boxes not to be subdivided
3817  areaSumLength = topw + frontw;
3818  // The surface are of the whole box
3819  areaWholeSA2 = width * areaSumLength + areaSplitPlane;
3820  // initial evaluation .. the worst cost
3821  bestCost = WorstEvaluation();
3822}
3823
3824// The initialization for the first axis to be tested, this time *****Y axis*******
3825// This can be run independently.
3826void
3827CKTBABuildUp::SSplitState::InitYaxis(int cnt, const SBBox &boxN)
3828{
3829  // initialize the variables, mainly for surface area estimates
3830  cntAll = cnt;     // the number of all objects in the bounding box
3831  cntLeft = 0;      // the count of bounding boxes on the left
3832  cntRight = cnt;   // the count of bounding boxes on the right
3833  thickness = 0;    // the count of bounding boxes straddling the splitting plane
3834  axis = CKTBAxes::EE_Y_axis; // the axis, where the splitting is proposed
3835  box = boxN;        // the box, that is subdivided
3836  //int frontAxis = oaxes[axis][0]; // axis that is considered in front .. height - x-axis
3837  //int topAxis = oaxes[axis][1]; // axis that is considered in top .. depth - y-axis
3838  // the size along the second axis - height
3839  frontw = box.Max().x - box.Min().x;
3840  // the size of bounding box along the axis to be split - width
3841  width = box.Max().y - box.Min().y;
3842  // The size along third axis - depth
3843  topw = box.Max().z - box.Min().z;
3844  // surface area of the splitting plane .. one face !!
3845  areaSplitPlane = topw * frontw;
3846  // the sum of length of the boxes not to be subdivided
3847  areaSumLength = topw + frontw;
3848  areaWholeSA2 = width * areaSumLength + areaSplitPlane;
3849  // initial evaluation .. the worst cost
3850  bestCost = WorstEvaluation();
3851}
3852
3853// The initialization for the first axis to be tested, this time *****Z axis*******
3854// This can be run independently.
3855void
3856CKTBABuildUp::SSplitState::InitZaxis(int cnt, const SBBox &boxN)
3857{
3858  // initialize the variables, mainly for surface area estimates
3859  cntAll = cnt;     // the number of all objects in the bounding box
3860  cntLeft = 0;      // the count of bounding boxes on the left
3861  cntRight = cnt;   // the count of bounding boxes on the right
3862  thickness = 0;    // the count of bounding boxes straddling the splitting plane
3863  axis = CKTBAxes::EE_Z_axis; // the axis, where the splitting is proposed
3864  box = boxN;       // the box, that is subdivided
3865  //int frontAxis = oaxes[axis][1]; // axis that is considered in front .. height - x axis
3866  //int topAxis = oaxes[axis][0]; // axis that is considered in top .. depth - y axis
3867  // The size along third axis - depth
3868  topw = box.Max().x - box.Min().x;
3869  // the size along the second axis - height
3870  frontw = box.Max().y - box.Min().y;
3871  // the size of bounding box along the axis to be split - width
3872  width = box.Max().z - box.Min().z;
3873  // surface area of the splitting plane .. one face !!
3874  areaSplitPlane = topw * frontw;
3875  // the sum of length of the boxes not to be subdivided
3876  areaSumLength = topw + frontw;
3877  areaWholeSA2 = width * areaSumLength + areaSplitPlane;
3878  // initial evaluation .. the worst cost
3879  bestCost = WorstEvaluation();
3880}
3881
3882// --------------------------------------------------------------------
3883// The initialization for the first axis to be tested, use only in case
3884// when all 3 axes are tested. It can be called only when InitXaxis was called before !!!!
3885void
3886CKTBABuildUp::SSplitState::ReinitYaxis(int cnt, const SBBox &boxN)
3887{
3888  // initialize the variables, mainly for surface area estimates
3889  assert(cntAll == cnt); // the number of all objects in the bounding box
3890  cntLeft = 0;      // the count of bounding boxes on the left
3891  cntRight = cnt;   // the count of bounding boxes on the right
3892  thickness = 0;    // the count of bounding boxes straddling the splitting plane
3893  axis = CKTBAxes::EE_Y_axis; // the axis, where the splitting is proposed
3894  //int topAxis = 2;      // axis that is considered in top .. depth
3895  //int frontAxis = 0;    // axis that is considered in front .. height 
3896  // Swap the width, height, and depth
3897  float tS = width;
3898  // the size of bounding box along the axis to be split
3899  width = topw;
3900  // and along two other axes
3901  topw = frontw;
3902  frontw = tS; // copy the temporary variable
3903  // surface area of the splitting plane .. one face !!
3904  areaSplitPlane = topw * frontw;
3905  // the sum of length of the boxes not to be subdivided
3906  areaSumLength = topw + frontw;
3907  //assert(areaWholeSA2 == width * SumLength + areaSplitPlane);
3908}
3909
3910// The initialization for the first axis to be tested, use only in case
3911// when all 3 axes are tested. It can be called only when ReinitYaxis was called before !!!
3912void
3913CKTBABuildUp::SSplitState::ReinitZaxis(int cnt, const SBBox &boxN)
3914{
3915  // initialize the variables, mainly for surface area estimates
3916  assert(cntAll == cnt); // the number of all objects in the bounding box
3917  cntLeft = 0;      // the count of bounding boxes on the left
3918  cntRight = cnt;   // the count of bounding boxes on the right
3919  thickness = 0;    // the count of bounding boxes straddling the splitting plane
3920  axis = CKTBAxes::EE_Z_axis; // the axis, where the splitting is proposed
3921  // Swap the width, height, and depth
3922  float tS = width;
3923  // the size of bounding box along the axis to be split
3924  width = topw;
3925  // and along two other axes 
3926  topw = frontw;
3927  frontw = tS; // copy the temporary variable
3928  // surface area of the splitting plane .. one face !!
3929  areaSplitPlane = topw * frontw;
3930  // the sum of length of the boxes not to be subdivided
3931  areaSumLength = topw + frontw;
3932  //assert(areaWholeSA2 == width * stateSumLength + areaSplitPlane);
3933}
3934
3935// This is the computation of the cost using surface area heuristics
3936// using LINEAR EXTIMATE
3937void
3938CKTBABuildUp::EvaluateCost(SSplitState &state)
3939{
3940  // the surface area of the bounding box for the left child
3941  assert(state.position > -1e-9);
3942  float areaLeft = state.position * state.areaSumLength + state.areaSplitPlane;
3943
3944  // the surface area of the right bounding box for the right child
3945  assert(state.width - state.position > -1e-9);
3946  float areaRight = (state.width - state.position) * state.areaSumLength +
3947                    state.areaSplitPlane;
3948
3949  // computation of the cost .. smaller is better
3950  float cost = areaLeft * state.cntLeft + areaRight * state.cntRight;
3951
3952#ifdef _DEBUG_COSTFUNCTION
3953  // Put there normalized position and cost also normalized somehow
3954  ReportCostStream(state.position / state.width,
3955                   cost / (state.areaWholeSA2 * state.cntAll));
3956#endif
3957 
3958  // This normalization is not in fact necessary here
3959  //cost //= state.areaWholeSA2;
3960  if (cost < state.bestCost) {
3961    state.bestCost = cost;
3962    // The iterator pointing to the best boundary
3963    state.bestIterator = state.it;
3964    // The number of objects whose boundaries are duplicated
3965    state.bestThickness = state.thickness;
3966  }
3967  return;
3968}
3969
3970// This is the computation of the cost using surface area heuristics
3971void
3972CKTBABuildUp::EvaluateCostFreeCut(SSplitState &state)
3973{
3974  // This is assumed to work for free cut (no object intersected)
3975  assert(state.thickness == 0);
3976 
3977  // the surface area of the bounding box for the left child
3978  assert(state.position > -1e-9);
3979  float areaLeft = state.position * state.areaSumLength + state.areaSplitPlane;
3980
3981  // the surface area of the right bounding box for the right child
3982  assert(state.width - state.position > -1e-9);
3983  float areaRight = (state.width - state.position) * state.areaSumLength +
3984                    state.areaSplitPlane;
3985
3986  // computation of the cost .. smaller is better
3987  float cost = biasFreeCuts *(areaLeft * state.cntLeft + areaRight * state.cntRight);
3988
3989#ifdef _DEBUG_COSTFUNCTION
3990  // Put there normalized position and cost also normalized somehow
3991  ReportCostStream(state.position / state.width,
3992                   cost / (state.areaWholeSA2 * state.cntAll));
3993#endif
3994 
3995  // This normalization is not in fact necessary here
3996  //cost //= state.areaWholeSA2;
3997  if (cost < state.bestCost) {
3998    state.bestCost = cost;
3999    // The iterator pointing to the best boundary
4000    state.bestIterator = state.it;
4001    // The number of objects whose boundaries are duplicated
4002    state.bestThickness = 0;
4003  }
4004  return;
4005}
4006
4007// -------------------------------------------------------------------------------
4008// TRIAL TO IMPROVE for a given X-axis search for the best splitting plane position.
4009// Starts from start item, the rest of the information is set to 'state' variable.
4010void
4011CKTBABuildUp::GetSplitPlaneOpt(SItemVec *vec, int axisToTest)
4012{
4013#ifdef _DEBUG
4014  if ((vec == NULL) || (vec->size() == 0)) {
4015    cerr << "Trying to subdivide some node without an object" << endl;
4016    abort();;
4017  }
4018#endif // _DEBUG 
4019
4020#if 0
4021  static int index = 0;
4022  cout << "index = " << index << endl;
4023  const int indexToFind = 8;
4024  if (index == indexToFind) {
4025    cout << "Index found = " << index << endl;
4026  } 
4027  index++;
4028#endif
4029 
4030  // some necessary space is required to cut off
4031  const float eps = 1.0e-6 * state.width;
4032  const float MinPosition = state.box.Min(axisToTest);
4033  const float MaxPosition = state.box.Max(axisToTest);
4034  //float MinPositionAccept = MinPosition + eps;
4035  //float MaxPositionAccept = MaxPosition - eps;
4036
4037  // wherefrom to start
4038  SItemVec::iterator curr = vec->begin();
4039  SItemVec::iterator next = vec->begin() + 1;
4040  // Set the type of splitting by this function, which is in this case
4041  // not overwritten in Evaluate() function
4042  state.bestTwoSplits = 1;
4043  float val = (*curr).pos;
4044  float nval;
4045  float pval = val;
4046
4047  // Evaluate the first possible position, cutting off empty
4048  // space on the left of the splitting plane
4049  state.position = val - MinPosition;
4050  if (state.position > 0.f) {
4051    state.position2 = next->pos - MinPosition;
4052    state.it = curr;
4053    // Here evaluate the cost of surface area heuristics
4054    EvaluateCostFreeCut(state);
4055  }
4056 
4057  //-------------------- TEST LOOP ---------------------------------------
4058  // make evaluation for each splitting position
4059  const int imax = vec->size()-1;
4060  int ii;
4061  for (ii = 0; ii < imax; ii++)
4062  {
4063    nval = next->pos;
4064
4065#if 0
4066    if (axisToTest == 2) {
4067      if ((val > 0.866) && (val < 0.870)) {       
4068        cout << "val = " << val;
4069        cout << "  nval = " << nval;
4070        if (curr->IsRightBoundary())
4071          cout << " R ";
4072        else
4073          cout << " L ";
4074        cout << " thickness = " << state.thickness << endl;
4075       
4076      }
4077    }
4078#endif   
4079    // update left box and rightbox properties
4080    if (curr->IsRightBoundary()) { // we are on the right boundary of an object
4081      state.cntRight--;
4082      state.thickness--;
4083      // Check possibly the position
4084      if ( (val != nval) ||
4085           ((curr->IsRightBoundary()) && (next->IsLeftBoundary())) )
4086        {
4087          state.position = val - MinPosition;
4088          state.position2 = nval - MinPosition;
4089          state.it = curr;
4090          // Here evaluate the cost of surface area heuristics
4091          EvaluateCost(state);
4092        }     
4093    }
4094    else {
4095      // the left boundary .. enter the object from left
4096
4097      // Check possibly the position
4098      if ( (pval != val)
4099           // || ((curr->IsRightBoundary()) && (next->IsLeftBoundary()))
4100           )
4101        {
4102          state.position = val - MinPosition;
4103          state.position2 = nval - MinPosition;
4104          state.it = curr;
4105          // Here evaluate the cost of surface area heuristics
4106          EvaluateCost(state);
4107        }   
4108     
4109      state.cntLeft++;
4110      state.thickness++;
4111    }
4112       
4113    curr = next; // pointer to the boundary
4114    next++;
4115    pval = val;
4116    val = nval; // value of splitting plane
4117  } // while
4118
4119  //cout << "i = " << i << endl;
4120 
4121  assert(state.cntLeft == state.cntAll);
4122  state.cntRight--;
4123  assert(state.cntRight == 0);
4124  state.thickness--;
4125  assert(state.thickness == 0); 
4126 
4127  // Evaluate last possible position
4128  state.position = val - MinPosition;
4129  if (state.position < MaxPosition) {
4130    state.position2 = state.width;
4131    state.it = curr;
4132    // Here evaluate the cost of surface area heuristics
4133    EvaluateCostFreeCut(state);
4134  }
4135   
4136  // In this case the best result is kept in 'state' variable
4137  return;
4138} // CTestAx::GetSplitPlaneOpt (for X, Y, Z)
4139
4140// -------------------------------------------------------------------------------
4141// TRIAL TO IMPROVE for a given X/Y/Z-axis search for the best splitting plane position.
4142// Starts from start item, the rest of the information is set to 'state' variable.
4143void
4144CKTBABuildUp::GetSplitPlaneOpt2(SItemVec *vec, int axisToTest)
4145{
4146#ifdef _DEBUG
4147  if ((vec == NULL) || (vec->size() == 0)) {
4148    cerr << "Trying to subdivide some node without an object" << endl;
4149    abort();;
4150  }
4151#endif // _DEBUG 
4152 
4153  // some necessary space is required to cut off
4154  const float eps = 1.0e-6 * state.width;
4155  const float MinPosition = state.box.Min(axisToTest);
4156  const float MaxPosition = state.box.Max(axisToTest);
4157  //float MinPositionAccept = MinPosition + eps;
4158  //float MaxPositionAccept = MaxPosition - eps;
4159
4160  // wherefrom to start
4161  SItemVec::iterator curr = vec->begin();
4162  SItemVec::iterator next = vec->begin() + 1;
4163  // Set the type of splitting by this function, which is in this case
4164  // not overwritten in Evaluate() function
4165  state.bestTwoSplits = 1;
4166  state.thickness = 0;
4167  float val = (*curr).pos;
4168  float nval;
4169  float pval = val;
4170
4171//#if 1
4172#ifdef _DEBUG_COSTFUNCTION
4173  static int indexToFind = 0;
4174  static int indexSS = 0;
4175  if ( (indexSS == indexToFindSS) &&
4176       (!_alreadyDebugged)) {
4177    cout << "Debugging cost function, index = " << indexSS << endl;
4178    InitCostStream(indexSS, state.cntAll,
4179                   (val-state.box.Min(axisToTest))/state.width,
4180                   //(val-MinPosition)/state.width,
4181                   //(MinPosition)/state.width,
4182                   (state.box.Max(axisToTest)-MinPosition)/state.width);
4183  }
4184  indexSS++;
4185#endif
4186 
4187  // Evaluate the first possible position, cutting off empty
4188  // space on the left of the splitting plane
4189  state.position = val - MinPosition;
4190  if (state.position > 0.f) {
4191    state.position2 = next->pos - MinPosition;
4192    state.it = curr;
4193    // Here evaluate the cost of surface area heuristics
4194    EvaluateCostFreeCut(state);
4195  }
4196 
4197  //-------------------- TEST LOOP ---------------------------------------
4198  // make evaluation for each splitting position
4199  const int imax = vec->size()-1;
4200  int ii;
4201  for (ii = 1; ii <= imax; ii++)
4202  {
4203    nval = (*vec)[ii].pos;
4204    // update left box and rightbox properties
4205    if (curr->IsRightBoundary()) {
4206      // the right boundary of an object
4207      state.cntRight--;
4208      state.thickness--;
4209     
4210      // Check possibly the position
4211      if ( (val != nval) ||
4212           ((curr->IsRightBoundary()) && (next->IsLeftBoundary())) )
4213        {
4214          if (state.thickness >= 0) {
4215            state.position = val - MinPosition;
4216            state.position2 = nval - MinPosition;
4217            state.it = curr;
4218            // Here evaluate the cost of surface area heuristics
4219            EvaluateCost(state);
4220          }
4221        }
4222    }
4223    else {
4224      // the left boundary .. enter the object from left
4225
4226      // Check possibly the position
4227      if (pval != val)
4228      {
4229        if (state.thickness >= 0) {
4230          state.position = val - MinPosition;
4231          state.position2 = nval - MinPosition;
4232          state.it = curr;
4233          // Here evaluate the cost of surface area heuristics
4234          EvaluateCost(state);
4235        }
4236      }     
4237      state.cntLeft++;
4238      state.thickness++;
4239    }
4240       
4241    curr = next; // pointer to the boundary
4242    pval = val;
4243    val = nval; // value of splitting plane
4244    next++;
4245  } // while
4246
4247  //cout << "i = " << i << endl;
4248 
4249  assert(state.cntLeft == state.cntAll);
4250  state.cntRight--;
4251  assert(state.cntRight == 0);
4252  state.thickness--;
4253  assert(state.thickness == 0); 
4254
4255#if 0
4256  if (state.thickness < 0) {
4257    cout << "SSSomething wrong happened here at end\n" << endl;
4258  }
4259#endif     
4260
4261#ifdef _DEBUG
4262  if (curr != vec->end()-1) {
4263    cerr << "Implementation bug in pointers" << endl;
4264    cout << "curr = " << curr - vec->begin()
4265         << " vec->end()-1 = " << vec->end()-1 - vec->begin() << endl;
4266    abort();;
4267  }
4268#endif
4269 
4270  // Evaluate last possible position
4271  state.position = val - MinPosition;
4272  if (state.position < MaxPosition) {
4273    state.position2 = state.width;
4274    state.it = curr;
4275    // Here evaluate the cost of surface area heuristics
4276    EvaluateCostFreeCut(state);
4277  }
4278 
4279#ifdef _DEBUG_COSTFUNCTION
4280  CloseCostStream();
4281#endif
4282 
4283  // In this case the best result is kept in 'state' variable
4284  return;
4285} // CTestAx::GetSplitPlaneOpt (for X, Y, Z)
4286
4287// -------------------------------------------------------------------------------
4288// TRIAL TO IMPROVE for a given X/Y/Z-axis search for the best splitting plane position.
4289// Using unrolling possibly understandable for inteligent compiler.
4290// Starts from start item, the rest of the information is set to 'state' variable.
4291void
4292CKTBABuildUp::GetSplitPlaneOpt3(SItemVec *vec, int axisToTest)
4293{
4294#ifdef _DEBUG
4295  if ((vec == NULL) || (vec->size() == 0)) {
4296    cerr << "Trying to subdivide some node without an object" << endl;
4297    abort();;
4298  }
4299#endif // _DEBUG 
4300
4301#if 0
4302  static int index = 0;
4303  cout << "index = " << index << endl;
4304  const int indexToFind = 27;
4305  if (index == indexToFind) {
4306    cout << "Index found = " << index << endl;
4307  } 
4308  index++;
4309#endif
4310 
4311  // some necessary space is required to cut off
4312  const float eps = 1.0e-6 * state.width;
4313  const float MinPosition = state.box.Min(axisToTest);
4314  const float MaxPosition = state.box.Max(axisToTest);
4315  //float MinPositionAccept = MinPosition + eps;
4316  //float MaxPositionAccept = MaxPosition - eps;
4317
4318  // wherefrom to start
4319  SItemVec::iterator curr = vec->begin();
4320  SItemVec::iterator next = vec->begin() + 1;
4321  // Set the type of splitting by this function, which is in this case
4322  // not overwritten in Evaluate() function
4323  state.bestTwoSplits = 1;
4324  float val = (*curr).pos;
4325  // value at previous and next position
4326  float pval, nval;
4327  pval = MinPosition;
4328
4329  // Evaluate the first possible position, cutting off empty
4330  // space on the left of the splitting plane
4331  state.position = val - MinPosition;
4332  // The first boundary must be left
4333  assert(curr->IsLeftBoundary());
4334  if (state.position > 0.f) {   
4335    state.position2 = next->pos - MinPosition;
4336    state.it = curr;
4337    // Here evaluate the cost of surface area heuristics
4338    EvaluateCostFreeCut(state);
4339  } // ------------------------------
4340  // Update for the next iteration
4341  state.cntLeft++;
4342  state.thickness++;
4343  pval = val;
4344 
4345  //-------------------- TEST LOOP ---------------------------------------
4346  // make evaluation for each splitting position, minus the first and
4347  // the last one
4348  int imax = vec->size()-2;
4349  // Setting unrolling
4350  const int MaxUNROLL = 4;
4351  const int loops = imax / MaxUNROLL;
4352  const int remainder = imax % MaxUNROLL;
4353  // There is some bug, when using this version, the boundaries are
4354  // incorrectly sorted.
4355  assert(remainder < MaxUNROLL);
4356  int ii = 1;
4357  if (loops > 0) {
4358    // The outer loop
4359    for (int l = 0; l < loops; l++, ii += MaxUNROLL)
4360    {
4361      // This has to be unrolled by compiler such as Intel Compiler or GCC4.1
4362      // The number of loops is known in advance !
4363      for (int qq = 0; qq < MaxUNROLL; qq++)
4364      {
4365        curr = vec->begin() + ii + qq;
4366        next = curr + 1;
4367        #if 1
4368        val = (*curr).pos;
4369        nval = (*next).pos;
4370        #else
4371        val = (*vec)[ii+qq].pos;
4372        nval = (*vec)[ii+qq+1].pos;
4373        #endif
4374       
4375        // update left box and rightbox properties
4376        if (curr->IsRightBoundary()){// we are on the right boundary of an object
4377          state.cntRight--;
4378          state.thickness--;
4379          // Check possibly the position
4380          if ( (val != nval) ||
4381               ((curr->IsRightBoundary()) && (next->IsLeftBoundary())) )
4382          {
4383            state.position = val - MinPosition;
4384            state.position2 = nval - MinPosition;
4385            state.it = curr;
4386            // Here evaluate the cost of surface area heuristics
4387            EvaluateCost(state);
4388            // It does make sense to change boundary value
4389            pval = (*vec)[ii+qq].pos;
4390          }
4391        }
4392        else {
4393          // the left boundary .. enter the object from left
4394       
4395          // Check possibly the position
4396          if (pval != val)
4397          {
4398            state.position = val - MinPosition;
4399            state.position2 = nval - MinPosition;
4400            state.it = curr;
4401            // Here evaluate the cost of surface area heuristics
4402            EvaluateCost(state);
4403            // It makes sense to change boundary value
4404            pval = (*vec)[ii+qq].pos;
4405          }
4406          state.cntLeft++;
4407          state.thickness++;
4408        } // left or right boundary     
4409      } // for qq
4410    } // for ii
4411
4412    // Set the correct pointer for the next iteration
4413    curr = next; next++;
4414    pval = val;
4415  }
4416  else {
4417    //cout << "No loop" << endl;
4418    pval = val;
4419    curr = next;
4420    next++;
4421  }
4422
4423  // The rest of the loop
4424  if (ii <= imax) {
4425    // Set the correct pointer 
4426    for ( ;ii <= imax; ii++) {
4427      nval = curr->pos;
4428      // update left box and rightbox properties
4429      if (curr->IsRightBoundary()) { // we are on the right boundary of an object
4430        state.cntRight--;
4431        state.thickness--;
4432        // Check possibly the position
4433        if ( (val != nval) ||
4434             ((curr->IsRightBoundary()) && (next->IsLeftBoundary())) )
4435        {
4436          state.position = val - MinPosition;
4437          state.position2 = nval - MinPosition;
4438          state.it = curr;
4439          // Here evaluate the cost of surface area heuristics
4440          EvaluateCost(state);
4441        }
4442      }
4443      else {
4444        // the left boundary .. enter the object from left
4445     
4446        // Check possibly the position
4447        if (pval != val)
4448        {
4449          state.position = val - MinPosition;
4450          state.position2 = nval - MinPosition;
4451          state.it = curr;
4452          // Here evaluate the cost of surface area heuristics
4453          EvaluateCost(state);
4454        }
4455        state.cntLeft++;
4456        state.thickness++;
4457      } // if right/left boundary
4458   
4459      curr = next; // pointer to the boundary
4460      pval = val;
4461      val = nval; // value of splitting plane
4462      next++;
4463    } // for qq
4464
4465    // For the last evaluation
4466  }
4467  else {
4468    if (!loops) {
4469      curr = next;
4470      val = nval;
4471    }
4472  }
4473
4474  //cout << "i = " << i << endl;
4475
4476#ifdef _DEBUG 
4477  if (curr != vec->end()-1) {
4478    cerr << "Implementation bug in pointers" << endl;
4479    cout << "curr = " << curr - vec->begin()
4480         << "vec->end()-1 = " << vec->end()-1 - vec->begin() << endl;
4481    abort();;
4482  }
4483#endif
4484 
4485  assert(state.cntLeft == state.cntAll);
4486  state.cntRight--;
4487  assert(state.cntRight == 0);
4488  state.thickness--;
4489  assert(state.thickness == 0);
4490 
4491  // Evaluate last possible position
4492  state.position = val - MinPosition;
4493  if (state.position < MaxPosition) {
4494    state.position2 = state.width;
4495    state.it = curr;
4496    // Here evaluate the cost of surface area heuristics
4497    EvaluateCostFreeCut(state);
4498  }
4499   
4500  // In this case the best result is kept in 'state' variable
4501  return;
4502} // CTestAx::GetSplitPlaneOpt3 (for X, Y, Z)
4503
4504#if 0
4505// Here is some bug, that is clearly visible for tree.nff !! The resulting
4506// tree is much less efficient. 3/1/2006 VH.
4507// -------------------------------------------------------------------------------
4508// TRIAL TO IMPROVE for a given X,Y,Z-axis search for the best splitting plane position.
4509// Starts from start item, the rest of the information is set to 'state' variable.
4510// It does not work for gcc compiler better than GetSplitPlaneOpt
4511void
4512CKTBABuildUp::GetSplitPlaneOptUnroll4(SItemVec *vec, int axisToTest)
4513{
4514#ifdef _DEBUG
4515  if ((vec == NULL) || (vec->size() == 0)) {
4516    cerr << "Trying to subdivide some node without an object" << endl;
4517    abort();;
4518  }
4519#endif // _DEBUG 
4520
4521#if 0
4522  static int index = 0;
4523  cout << "index = " << index << endl;
4524  const int indexToFind = 8;
4525  if (index == indexToFind) {
4526    cout << "Index found = " << index << endl;
4527  } 
4528  index++;
4529#endif
4530 
4531  // some necessary space is required to cut off
4532  const float eps = 1.0e-6 * state.width;
4533  const float MinPosition = state.box.Min(axisToTest);
4534  const float MaxPosition = state.box.Max(axisToTest);
4535  //float MinPositionAccept = MinPosition + eps;
4536  //float MaxPositionAccept = MaxPosition - eps;
4537
4538  // wherefrom to start
4539  // Set the type of splitting by this function, which is in this case
4540  // not overwritten in Evaluate() function
4541  state.bestTwoSplits = 1;
4542  float val4 = MinPosition;
4543 
4544  //-------------------- TEST LOOP ---------------------------------------
4545  // make evaluation for each splitting position
4546
4547  // the first position was already evaluated and the last is also a special case
4548  const int imax = vec->size() - 1;
4549  int cntDebug = 0;
4550  int imax4 = 4*(imax/4);
4551  int imax4rest = imax % 4;
4552  int i;
4553  // This is manual unrolling
4554  if (imax4 >= 4) {
4555  for(i = 0; i < imax4; i += 4)
4556  {
4557    // The evaluation I
4558    cntDebug++;
4559    float pval1 = val4;
4560    SItem* curr1 = &((*vec)[i+0]); float val1 = (*vec)[i+0].pos;
4561    SItem* next1 = &((*vec)[i+1]); float nval1 = (*vec)[i+1].pos;
4562    if (curr1->IsRightBoundary()) {
4563      state.cntRight--; state.thickness--;
4564      if ( (val1 != nval1) ||
4565           ((curr1->IsRightBoundary()) && (next1->IsLeftBoundary())) ) {
4566        state.position = val1 - MinPosition; state.position2 = nval1 - MinPosition;
4567        state.it = vec->begin() + i; EvaluateCost(state);
4568      }
4569    }
4570    else {
4571      if (pval1 != val1) {
4572          state.position = val1 - MinPosition; state.position2 = nval1 - MinPosition;
4573          state.it = vec->begin() + i; EvaluateCost(state);
4574      }
4575      state.cntLeft++; state.thickness++;
4576    }
4577    // The evaluation II
4578    cntDebug++;
4579    float pval2 = (*vec)[i+0].pos;
4580    SItem* curr2 = &((*vec)[i+1]); float val2 = (*vec)[i+1].pos;
4581    SItem* next2 = &((*vec)[i+2]); float nval2 = (*vec)[i+2].pos;
4582    if (curr2->IsRightBoundary()) {
4583      state.cntRight--; state.thickness--;
4584      if ( (val2 != nval2) ||
4585           ((curr2->IsRightBoundary()) && (next2->IsLeftBoundary())) ) {
4586        state.position = val2 - MinPosition; state.position2 = nval2 - MinPosition;
4587        state.it = vec->begin() + i + 1; EvaluateCost(state);
4588      }
4589    }
4590    else {
4591      if (pval2 != val2) {
4592          state.position = val2 - MinPosition; state.position2 = nval2 - MinPosition;
4593          state.it = vec->begin() + i + 1; EvaluateCost(state);
4594      }
4595      state.cntLeft++; state.thickness++;
4596    }
4597    // The evaluation III
4598    cntDebug++;
4599    float pval3 = (*vec)[i+1].pos;
4600    SItem* curr3 = &((*vec)[i+2]); float val3 = (*vec)[i+2].pos;
4601    SItem* next3 = &((*vec)[i+3]); float nval3 = (*vec)[i+3].pos;
4602    if (curr3->IsRightBoundary()) {
4603      state.cntRight--; state.thickness--;
4604      if ( (val3 != nval3) ||
4605           ((curr3->IsRightBoundary()) && (next3->IsLeftBoundary())) ) {
4606        state.position = val3 - MinPosition; state.position2 = nval3 - MinPosition;
4607        state.it = vec->begin() + i + 2; EvaluateCost(state);
4608      }
4609    }
4610    else {
4611      if (pval3 != val3) {
4612          state.position = val3 - MinPosition; state.position2 = nval3 - MinPosition;
4613          state.it = vec->begin() + i + 2; EvaluateCost(state);
4614      }
4615      state.cntLeft++; state.thickness++;
4616    }
4617    // The evaluation IV
4618    cntDebug++;
4619    float pval4 = (*vec)[i+2].pos;
4620    SItem* curr4 = &((*vec)[i+3]); val4 = (*vec)[i+3].pos;
4621    SItem* next4 = &((*vec)[i+4]); float nval4 = (*vec)[i+4].pos;
4622    if (curr4->IsRightBoundary()) {
4623      state.cntRight--; state.thickness--;
4624      if ( (val4 != nval4) ||
4625           ((curr4->IsRightBoundary()) && (next4->IsLeftBoundary())) ) {
4626        state.position = val4 - MinPosition; state.position2 = nval4 - MinPosition;
4627        state.it = vec->begin() + i + 3; EvaluateCost(state);
4628      }
4629    }
4630    else {
4631      if (pval4 != val4) {
4632          state.position = val4 - MinPosition; state.position2 = nval4 - MinPosition;
4633          state.it = vec->begin() + i + 3; EvaluateCost(state);
4634      }
4635      state.cntLeft++; state.thickness++;
4636    }
4637  } // for, unrolled loop by factor 4
4638  } // if more than 4 evaluations
4639 
4640  for(i = imax4; i < imax; i ++)
4641  {
4642    cntDebug++;
4643    float pvalI = val4;
4644    SItem* currI = &((*vec)[i+0]); float valI = (*vec)[i+0].pos;
4645    SItem* nextI = &((*vec)[i+1]); float nvalI = (*vec)[i+1].pos;
4646    if (currI->IsRightBoundary()) {
4647      state.cntRight--; state.thickness--;
4648      if ( (valI != nvalI) ||
4649           ((currI->IsRightBoundary()) && (nextI->IsLeftBoundary())) ) {
4650        state.position = valI - MinPosition; state.position2 = nvalI - MinPosition;
4651        state.it = vec->begin() + i; EvaluateCost(state);
4652      }
4653    }
4654    else {
4655      if (pvalI != valI) {
4656          state.position = valI - MinPosition; state.position2 = nvalI - MinPosition;
4657          state.it = vec->begin() + i; EvaluateCost(state);
4658      }
4659      state.cntLeft++; state.thickness++;
4660    }
4661    val4 = valI;
4662  } // for i
4663 
4664  //cout << "i = " << i << endl;
4665  assert(cntDebug == imax);
4666 
4667  assert(state.cntLeft == state.cntAll);
4668  state.cntRight--;
4669  assert(state.cntRight == 0);
4670  state.thickness--;
4671  assert(state.thickness == 0); 
4672 
4673  // Evaluate last possible position
4674  state.position = (*vec)[imax-1].pos - MinPosition;
4675  if (state.position < MaxPosition) {
4676    state.position2 = state.width;
4677    state.it = vec->end() - 1;
4678    // Here evaluate the cost of surface area heuristics
4679    EvaluateCost(state);
4680  }
4681   
4682  // In this case the best result is kept in 'state' variable
4683  return;
4684} // CTestAx::GetSplitPlaneX,Y,Z
4685#endif
4686
4687}
4688
Note: See TracBrowser for help on using the repository browser.