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

Revision 2592, 135.9 KB checked in by bittner, 17 years ago (diff)

havran ray caster update

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 = 50;
221  stackDepth = 2*maxTreeDepth;
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.2 + 2.0);
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
1208  if (verbose)
1209    cout << "MaxListLength = " << maxListLength << endl; 
1210
1211  return data;
1212}
1213
1214// interface function for building up CKTB tree
1215SKTBNodeT*
1216CKTBABuildUp::BuildUp(ObjectContainer &objlist,
1217                      const AxisAlignedBox3 &box,
1218                      bool verboseF)
1219{
1220  // check the number of objects
1221  if (objlist.size() == 0)
1222    return 0; // nothing
1223  // cerr<<"hh44"<<endl;
1224  // initialize the whole box of the kd-tree and sort
1225  // the boundary entries
1226  SInputData *d = Init(&objlist, box);
1227
1228  //cerr<<"hh45"<<endl;
1229
1230  // copy verbose to be used consistenly further on
1231  this->verbose = verboseF;
1232 
1233  if (verbose)
1234    cout << "Building up KTB tree for " << objlist.size()
1235                 << " objects." << endl << flush;
1236
1237  // Initialization of the first value to insert the min boxes
1238  d->lastMinBoxSA2 = box.SurfaceArea() * 10000.f;
1239  d->lastDepthForMinBoxes = -20; // for root you always put the node
1240  d->lastMinBoxNode = 0;
1241 
1242  // -----------------------------------------------
1243  // Start to build the tree
1244  if ( (cutEmptySpace) &&
1245       (initcnt <= maxListLength)) {
1246    // only cutting off empty space for initial leaf
1247    d->modeSubDiv = EE_SUBDIVCUTEMPTY;
1248    root = SubDiv(d);
1249  }
1250  else {
1251    // normal subdivision scheme, creating whole tree
1252    root = SubDiv(d);
1253  }
1254
1255#ifdef _DEBUG
1256  if (GetDepth() != 0) {
1257    cerr << "Using depth value does not work correctly, depth = "
1258          << GetDepth() << endl;
1259  }
1260#endif 
1261  assert(root); 
1262
1263  // ----------------------------------------------------
1264  // Deallocate all auxiliary data
1265  DeleteAuxiliaryData();
1266
1267  // delete the temporary data array
1268  solidArray.resize(0);
1269 
1270  // the pointer to the root node
1271  return GetRootNode();
1272}
1273
1274int
1275CKTBABuildUp::UpdateEvaluation(float &eval, const float &newEval)
1276{
1277  if (newEval < eval) {
1278    eval = newEval;
1279    return 1; // updated
1280  }
1281  return 0; // not updated
1282}
1283
1284// recursive function for subdividing the CKTB tree into two
1285// halves .. creates one node of CKTB tree
1286// list .. list of boundaries, count .. number of objects in current node
1287// bb .. current bounding box of node, (pars,reqGlobAxis, modeSubDiv) .. cut pref.
1288SKTBNodeT*
1289CKTBABuildUp::SubDiv(SInputData *d)
1290{
1291#ifdef _DEBUG
1292  Check3List(d);
1293#endif
1294
1295#if 0
1296  static int index = 0;
1297  cout << "SubDiv index = " << index << endl;
1298  const int indexToFind = 13916;
1299  if (index == indexToFind) {
1300    cout << " IndexToFind = " << indexToFind << endl;
1301    cout << " You should start debug" << endl;
1302  }
1303  index++;
1304#endif
1305
1306  // This is the node to return from this function
1307  SKTBNodeT *nodeToReturn = 0;
1308
1309  // if to make the interior node extended by the box here
1310  makeMinBoxHere = false;
1311
1312  // -----------------------------------------------------
1313  if (makeMinBoxes) {
1314#if 1
1315    if ( (d->count >= minObjectsToCreateMinBox) &&
1316         (GetDepth() - d->lastDepthForMinBoxes >= minDepthDistanceBetweenMinBoxes) )  {
1317        makeMinBoxHere = true;
1318        d->lastDepthForMinBoxes = GetDepth();
1319        d->lastMinBoxSA2 = d->box.SA2();
1320    }
1321#endif
1322#if 0
1323    if (d->count >= minObjectsToCreateMinBox) {
1324      if ( (d->lastMinBoxSA2 >= d->box.SA2() * minSA2ratioMinBoxes) &&
1325           (GetDepth() - d->lastDepthForMinBoxes >= minDepthDistanceBetweenMinBoxes) ) {
1326        makeMinBoxHere = true;
1327        d->lastDepthForMinBoxes = GetDepth();
1328        d->lastMinBoxSA2 = d->box.SA2();
1329      }
1330    }
1331#endif
1332#if 0
1333    if ( (d->count >= minObjectsToCreateMinBox) &&
1334         (GetDepth() % minDepthDistanceBetweenMinBoxes == 0) ) {
1335      makeMinBoxHere = true;
1336      d->lastDepthForMinBoxes = GetDepth();
1337      d->lastMinBoxSA2 = d->box.SA2();
1338    }
1339#endif
1340  }
1341
1342  if ( (makeMinBoxHere) && (makeTightMinBoxes) ) {
1343    // In case we should make the box of the current
1344    // interior node and this should be tight, we have
1345    // to first compute its extent before splitting
1346    SBBox tightbox;
1347    GetTightBox(*d, tightbox);
1348    // Here we decrease the box to the minimum size
1349    d->box = tightbox;
1350  }
1351 
1352  // ------------------------------------------------------------------
1353  // if we are forced to make subdivision at given point
1354  if (d->pars.reqPosition != Limits::Infinity) {
1355    // this code preceeds switch(mode) since 2-level evaluation
1356    d->position = d->pars.reqPosition;
1357    d->axis = d->pars.reqAxis;
1358    // next subdivison will not be forced
1359    d->pars.reqPosition = Limits::Infinity;
1360    d->pars.reqAxis = CKTBAxes::EE_Leaf;
1361    IncDepth();
1362    SKTBNodeT* n = MakeOneCut(d);
1363    if (!nodeToReturn)
1364      nodeToReturn = n;
1365    DecDepth();
1366    return nodeToReturn;
1367  }
1368 
1369  // -------------------------------------------------------------------------
1370  // determine between subdivision modes - different termination criteria
1371  switch (d->modeSubDiv) { // subdivision mode
1372    case EE_SUBDIVCUTEMPTY: { // cutting off empty space
1373      if ((GetDepth() >= absMaxAllowedDepth) ||
1374          (GetDepth() >= (startEmptyCutDepth + maxEmptyCutDepth)) ||
1375          (d->count == 0) ) {
1376        SKTBNodeT* n = MakeLeaf(d);
1377        if (!nodeToReturn)
1378          nodeToReturn = n;
1379        // cout << "LEC\n";
1380        return nodeToReturn;
1381      }
1382      // to require the axis has no meaning in cutting off
1383      d->pars.reqAxis = d->axis = CKTBAxes::EE_Leaf;
1384      break;
1385    }
1386    // both adhoc, clustering and automatic termination criteria
1387    case EE_SUBDIVADHOC:
1388    case EE_SUBDIVAUTOMATIC: {
1389      // automatic termination criteria are used in this phase as
1390      // the absolute threshold similarly to SUBDIVADHOC mode
1391      if ( (GetDepth() >= maxDepthAllowed) ||
1392           (d->count <= maxListLength)) {
1393        SKTBNodeT* n = MakeLeaf(d);
1394        if (!nodeToReturn)
1395          nodeToReturn = n;
1396        // cout << "LLA\n";
1397        return nodeToReturn;
1398      }
1399      break;
1400    }
1401    default: {
1402      cerr << " unknown subdivision mode in ibsp.cpp\n";
1403      abort();;
1404    }
1405  } // switch(modeSubDiv)
1406 
1407  // the axis where splitting should be done
1408  CKTBAxes::Axes axis = CKTBAxes::EE_Leaf;
1409  d->twoSplits = -1;
1410  d->bestCost = WorstEvaluation() * 0.5f;
1411  d->position = MAXFLOAT;
1412  d->cntThickness = 0;
1413
1414// which calls should be used, standard or optimized
1415#if 1
1416#define CallGetSplitPlaneX GetSplitPlaneOpt2
1417#define CallGetSplitPlaneY GetSplitPlaneOpt2
1418#define CallGetSplitPlaneZ GetSplitPlaneOpt2
1419#endif
1420#if 0
1421// This does not work for unknown reason... VH 2006/01/08
1422#define CallGetSplitPlaneX GetSplitPlaneOpt3
1423#define CallGetSplitPlaneY GetSplitPlaneOpt3
1424#define CallGetSplitPlaneZ GetSplitPlaneOpt3
1425#endif
1426#if 0
1427#define CallGetSplitPlaneX GetSplitPlaneOptUnroll4
1428#define CallGetSplitPlaneY GetSplitPlaneOptUnroll4
1429#define CallGetSplitPlaneZ GetSplitPlaneOptUnroll4
1430#endif
1431 
1432  // if the axis is required in pars structure for some reasons
1433  if (d->pars.useReqAxis) {
1434    axis = d->pars.reqAxis;
1435    // The next subdivision is not prescribed
1436    d->pars.useReqAxis = false;
1437    switch(axis) {
1438      case CKTBAxes::EE_X_axis : {
1439        state.InitXaxis(d->count, d->box); // init
1440        CallGetSplitPlaneX(d->xvec, 0); // evaluate
1441        break;
1442      }
1443      case CKTBAxes::EE_Y_axis : {
1444        state.InitYaxis(d->count, d->box); // init
1445        CallGetSplitPlaneY(d->yvec, 1); // evaluate
1446        break;
1447      }
1448      case CKTBAxes::EE_Z_axis : {
1449        state.InitZaxis(d->count, d->box); // init
1450        CallGetSplitPlaneZ(d->zvec, 2); // evaluate
1451        break;
1452      }
1453      default: {
1454        cerr << "No other option is allowed here" << endl;
1455        abort();;
1456      }
1457    }
1458    // compute the cost, normalized by surface area
1459    state.NormalizeCostBySA2();
1460    d->bestCost /= state.areaWholeSA2;
1461    if (UpdateEvaluation(d->bestCost, state.bestCost)) {
1462      // Compute correct position to be used for splitting
1463      d->position = (*(state.bestIterator)).pos;
1464      d->bestIterator = state.bestIterator;
1465      d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
1466      d->cntThickness = state.bestThickness;
1467      if (d->twoSplits > 1) {
1468        SItemVec::iterator it = state.bestIterator; it++;     
1469        if (it != d->GetItemVec(axis)->end())
1470          d->position2 = (*it).pos;
1471        else
1472          d->position2 = state.box.Max(axis);
1473      }
1474    }
1475  }
1476  else {
1477    int algorithmForAxisSel = _algorithmForAxisSelection;
1478    if (d->modeSubDiv == EE_SUBDIVCUTEMPTY)
1479      algorithmForAxisSel = 0;
1480
1481    // only a single axis will be tested
1482    bool useSingleAxis = true;
1483       
1484    // the required axis is by global function .. e.g. cyclic change x, y, z
1485    switch (algorithmForAxisSel) {
1486      case 0: {
1487        // all three axis will be used, starting from x, y, z
1488        useSingleAxis = false;
1489        break;
1490      }
1491      case 1: {
1492        // cyclic order of axes, x, y, z, x, y....
1493        axis = d->pars.reqAxis;
1494        // compute the axis for the next subdivision step
1495        d->pars.reqAxis = oaxes[axis][0];
1496        break;
1497      }
1498      case 2 : {
1499        // The next time will be determined the same way
1500        axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
1501        break;
1502      }
1503      case 3 : {
1504        float p = RandomValue(0.0f, 1.0f);
1505        const float thresholdAxisP = 0.8f;
1506        if (p < thresholdAxisP) {         
1507          axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
1508        }
1509        else {
1510          // Use the axis prescribed from previous step
1511          axis = d->pars.reqAxis;
1512        }
1513        // for the next time, different axis
1514        d->pars.reqAxis = oaxes[axis][0];
1515        break;
1516      }
1517      case 4 : {
1518        float p = RandomValue(0.0f, 1.0f);
1519        const float thresholdAxisP = 0.3f;
1520        if (p < thresholdAxisP) {         
1521          axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
1522        }
1523        else {
1524          // Use the axis prescribed from previous step
1525          axis = d->pars.reqAxis;
1526        }
1527        // for the next time, different axis
1528        d->pars.reqAxis = oaxes[axis][0];
1529        break;
1530      }
1531      case 5 : {
1532        float p = RandomValue(0.0f, 1.0f);
1533        const float thresholdAxisP = 0.2f;
1534        d->pars.reqAxis = (CKTBAxes::Axes)-1;
1535        if (p < thresholdAxisP) {
1536          axis = (CKTBAxes::Axes)(d->box.Diagonal().DrivingAxis());
1537          // for the next time, different axis
1538        }
1539        else {
1540          // Use the axis for which cost model gets minimum,
1541          // compute it for all axes x, y, z
1542          useSingleAxis = false;
1543        }
1544        break;
1545      }
1546      default: {
1547        cerr << "Selection algorithm = " << algorithmForAxisSel
1548              << " for axis is not implemented" << endl;
1549        abort();;
1550      }
1551    } // switch
1552
1553    // How the algorithm is used
1554    if (useSingleAxis) {
1555      // -------------------------------------------
1556      // Compute subdivision only for one axis
1557      switch(axis) {
1558        case CKTBAxes::EE_X_axis : {
1559          state.InitXaxis(d->count, d->box); // init
1560          CallGetSplitPlaneX(d->xvec, 0); // evaluate
1561          break;
1562        }
1563        case CKTBAxes::EE_Y_axis : {
1564          state.InitYaxis(d->count, d->box); // init
1565          CallGetSplitPlaneY(d->yvec, 1); // evaluate
1566          break;
1567        }
1568        case CKTBAxes::EE_Z_axis : {
1569          state.InitZaxis(d->count, d->box); // init
1570          CallGetSplitPlaneZ(d->zvec, 2); // evaluate
1571          break;
1572        }
1573        //case CKTBAxes::EE_Leaf : {
1574        //goto MAKE_LEAF;
1575        //}
1576        default: {
1577          cerr << "Selection algorithm = " << algorithmForAxisSel
1578                << " for axis = " << axis <<" is not implemented" << endl;
1579          abort();;
1580        }
1581      }
1582      state.NormalizeCostBySA2();
1583      d->bestCost /= state.areaWholeSA2;
1584      if (UpdateEvaluation(d->bestCost, state.bestCost)) {
1585        // Compute correct position to be used for splitting
1586        d->position = (*(state.bestIterator)).pos;
1587        d->bestIterator = state.bestIterator;
1588        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
1589        d->cntThickness = state.bestThickness;
1590        if (d->twoSplits > 1) {
1591          SItemVec::iterator it = state.bestIterator; it++;     
1592          if (it != d->GetItemVec(axis)->end())
1593            d->position2 = (*it).pos;
1594          else
1595            d->position2 = state.box.Max(axis);
1596        }
1597      }
1598    }
1599    else {
1600      // no axis is required, we can pickup any we like. We test all three axis
1601      // and we pickup the minimum cost for all three axes.
1602      state.InitXaxis(d->count, d->box); // init for x axis
1603      CallGetSplitPlaneX(d->xvec, 0); // evaluate for x axis
1604      if (UpdateEvaluation(d->bestCost, state.bestCost)) {     
1605        axis = CKTBAxes::EE_X_axis; // the x-axis should be used
1606        // Compute correct position to be used for splitting
1607        d->bestCost = state.bestCost;
1608        d->position = (*(state.bestIterator)).pos;
1609        d->bestIterator = state.bestIterator;
1610        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
1611        d->cntThickness = state.bestThickness;
1612        if (d->twoSplits > 1) {
1613          SItemVec::iterator it = state.bestIterator;
1614          it++;
1615          if (it != d->xvec->end())
1616            d->position2 = (*it).pos;
1617          else
1618            d->position2 = state.box.Max().x;
1619        }
1620      }
1621      // -----------------------
1622      // now test for y axis
1623      state.ReinitYaxis(d->count, d->box); // init for x axis
1624      //state.InitYaxis(d->count, d->box); // init for x axis
1625      CallGetSplitPlaneY(d->yvec, 1); // evaluate for x axis
1626      if (UpdateEvaluation(d->bestCost, state.bestCost)) {
1627        axis = CKTBAxes::EE_Y_axis; // the y-axis should be used
1628        // Compute correct position to be used for splitting
1629        d->position = (*(state.bestIterator)).pos;
1630        d->bestIterator = state.bestIterator;
1631        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
1632        d->cntThickness = state.bestThickness;
1633        if (d->twoSplits > 1) {
1634          SItemVec::iterator it = state.bestIterator;
1635          it++;
1636          if (it != d->yvec->end())
1637            d->position2 = (*it).pos;
1638          else
1639            d->position2 = state.box.Max().y;
1640        }
1641      }
1642      // -----------------------
1643      // now test for z axis
1644      state.ReinitZaxis(d->count, d->box); // init for x axis
1645      //state.InitZaxis(d->count, d->box); // init for x axis
1646      CallGetSplitPlaneZ(d->zvec, 2); // evaluate for x axis
1647      if (UpdateEvaluation(d->bestCost, state.bestCost)) {
1648        axis = CKTBAxes::EE_Z_axis; // the z-axis should be used
1649        // Compute correct position to be used for splitting
1650        d->position = (*(state.bestIterator)).pos;
1651        d->bestIterator = state.bestIterator;
1652        d->twoSplits = state.bestTwoSplits; // one or two splitting planes planed
1653        d->cntThickness = state.bestThickness;
1654        if (d->twoSplits > 1) {
1655          SItemVec::iterator it = state.bestIterator;
1656          it++;
1657          if (it != d->zvec->end())
1658            d->position2 = (*it).pos;
1659          else
1660            d->position2 = state.box.Max().z;
1661        }
1662      }
1663      // Now we have to renormalize the cost for purpose of termination the building
1664      d->bestCost /= state.areaWholeSA2;
1665    } // arbitrary order of axes
1666  } // for which axis to compute surface area heuristics
1667
1668  // ----------------------------------------------------------
1669  // when automatic termination criteria should be used
1670  if (d->modeSubDiv == EE_SUBDIVAUTOMATIC) {
1671    bool stopSubdivision = false;
1672    float ratio;
1673    // which algorithm to use for automatic termination criteria
1674    switch (algorithmAutoTermination) {
1675      // --------------------------------------------------
1676      case 0: { // This is described in my thesis.
1677        // This is the cost of unsubdivided node
1678        float leafCost = (float)(d->count);
1679        ratio = d->bestCost / leafCost;
1680        const float minRatio = 0.75;
1681        if (ratio > minRatio)
1682          stopSubdivision = true;
1683        break;
1684      }
1685      // --------------------------------------------------
1686      case 1: {
1687        // This is the cost of unsubdivided node, with uniformly distributed
1688        // objects and the spatial median subdivision, without any object straddling
1689        // the splitting plane. In some sense ideal setting for uniform distribution.
1690        Vector3 s = d->box.Diagonal();
1691        // Taking the widest side of the object to be subdivided
1692        int diagAxis = s.DrivingAxis();
1693        // Half of the width
1694        // ----------------------
1695        float w2 = s[diagAxis] * 0.5;
1696        float t = s[oaxes[diagAxis][0]];
1697        float h = s[oaxes[diagAxis][1]];
1698        float sah2 = w2*(t+h) + h*t; // half of the surface area of a child
1699        // This is the cost for case when the node is subdivided in the half and
1700        // the half of the object is on the left and half of the objects on the right
1701        // This is not the best cost possible !!!
1702        // = float subdividedCost = 2 * (sah2/state.areaWholeSA2 * count/2);
1703        float subdividedCost = sah2/state.areaWholeSA2 * (float)(d->count);
1704        ratio = d->bestCost / subdividedCost;
1705        // This is the max ratio allowed for splitting node subdivision
1706
1707        // The computed ratio = 1.0 corresponds to 'ideal case', so the threshold
1708        // must be higher than 1.0 !
1709        const float minRatio = 1.1;
1710        if (ratio > minRatio)
1711          stopSubdivision = true;
1712        break;
1713      }
1714      default : {
1715        cerr << "Uknown algorithm for automatic termination criteria = "
1716              << algorithmAutoTermination << endl;
1717      }
1718       
1719    } // switch
1720   
1721    //cout << "R=" << ratio << " ";
1722    if (stopSubdivision) {
1723      // cout << "F" << pars.failedSubDivCount << " ";
1724      // when small improvement by the subdivision is reached
1725      d->pars.failedSubDivCount++;
1726      if (d->pars.failedSubDivCount > maxCountTrials) {
1727        // make leaf and finish with this KD-tree branch
1728        // cout << "L, cnt=" << count << " ,d=" << GetDepth()()
1729        //     << " ,c=" << bestQ << " ,rat=" << ratio << "\n";
1730        // possibly cut empty space before leaf is created
1731        if (cutEmptySpace) {
1732          // setting initial depth if empty space cutting
1733          startEmptyCutDepth = GetDepth();
1734          IncDepth();
1735          d->modeSubDiv = EE_SUBDIVCUTEMPTY;
1736          // create the leaf first, with late empty cutting
1737          SKTBNodeT *n = SubDiv(d);
1738          if (!nodeToReturn)
1739            nodeToReturn = n;
1740          DecDepth();
1741          return nodeToReturn; // and finish by constructing a leaf
1742        }
1743
1744        // make leaf and finish
1745        SKTBNodeT *n = MakeLeaf(d);
1746        if (!nodeToReturn)
1747          nodeToReturn = n;
1748        return nodeToReturn;
1749      }
1750    }
1751
1752    // Do not stop subdivision now
1753   
1754    // update the progress in the parameters
1755    d->pars.ratioLastButOne = d->pars.ratioLast;
1756    d->pars.ratioLast = ratio;
1757
1758    // assure that axis is defined, if finding the
1759    // splitting plane has failed
1760    if ( (axis == CKTBAxes::EE_Leaf) ||
1761         (d->position == MAXFLOAT) ) {
1762      // compute the spatial median for arbitrary axis
1763      axis = (CKTBAxes::Axes)int(RandomValue(0.f, 2.98f));
1764      // set position based algorithm for the next cut and children
1765      d->algorithmBreakAx = 0;
1766      resetFlagsForBreakAx = true;
1767      d->position = (d->box.Min(axis) + d->box.Max(axis)) * 0.5f;
1768      d->twoSplits = 1;
1769    }
1770  } // subdivision automatic
1771
1772  // -----------------------------------------------
1773  // We have evaluated surface area heuristics above
1774  // if no winner for subdivision exists, make a leaf
1775  if ( (axis == CKTBAxes::EE_Leaf) ||
1776       (d->position == MAXFLOAT) )
1777  {
1778    //cerr << "This shall not happen" << endl;
1779    SKTBNodeT *n = MakeLeaf(d);
1780    if (!nodeToReturn)
1781      nodeToReturn = n;
1782    // cout << "LL\n";
1783    return nodeToReturn;
1784  }
1785 
1786  if (verbose) {
1787#ifdef _DEBUG
1788
1789//#define _KTBPRINTCUTS
1790#ifndef _KTBPRINTCUTS
1791    if (GetDepth() == 0)
1792#else
1793    if (_printCuts)
1794#endif
1795    {
1796      cout << "position: " << d->position << ", axis: "
1797             << (int)axis << ", depth=" << GetDepth() << ", box:" << endl;
1798      Describe(d->box, cout, 0);
1799      cout << endl;
1800    }
1801#endif // _DEBUG
1802  }
1803
1804#ifdef _DEBUG
1805  if (d->twoSplits == -1) {
1806    cerr << "Some problem in implementation" << endl;
1807    abort();;
1808  }
1809#endif
1810
1811  // copy the parameters to use for splitting
1812  d->axis = axis;
1813
1814  if (d->twoSplits == 1) {
1815    // create interior node with two successors
1816    // case (1)
1817    // DEBUG << "case 1 \n";
1818    // Make easy cut, using one position
1819    IncDepth();
1820    SKTBNodeT* n = MakeOneCut(d);
1821    if (!nodeToReturn)
1822      nodeToReturn = n;
1823    DecDepth();
1824    return nodeToReturn;
1825  }
1826
1827#if 1
1828  cerr << "Unexpected use in this implementation" << endl;
1829  cerr << "ABORTING " << endl;
1830  abort();;
1831#endif
1832 
1833//  case       case
1834//  (2)   or   (3) structure must be created
1835//   O          O                #
1836//  / \        / \               #
1837// L  RI      LI  R              #
1838//   /  \    / \                 #
1839//   RL RR  LL LR                #
1840
1841#ifdef _DEBUG
1842  if ( (d->box.Min()[axis] >= d->position) ||
1843       (d->position >= d->position2) ||
1844       (d->position2 >= d->box.Max()[axis]) ) {
1845    cerr << " the case 2 or 3 is defined incorrectly\n";
1846    abort();;
1847  }
1848#endif
1849
1850  // To be reworked !!!
1851  abort();;
1852  if (d->twoSplits == 2) {
1853    // -----------------------------------------
1854    // case (2)
1855    // DEBUG << "case 2 \n";
1856
1857    // make L part, linked in DFS order
1858    int cntLeft = 1;
1859    int cntRight = 0;
1860    IncDepth();
1861    SKTBNodeT *n = MakeOneCut(d);
1862    if (!nodeToReturn)
1863      nodeToReturn = n;
1864
1865    // make RI part, linked to the left node explictly
1866    SBBox rbb = d->box; // the bounding box
1867    rbb.Reduce(axis, 1, d->position);
1868    d->box = rbb;
1869    d->pars.reqAxis = axis;
1870    d->pars.reqPosition = d->position2;   
1871    SKTBNodeT *n2 = SubDiv(d);
1872    DecDepth();
1873    //Link the node
1874    SetInteriorNodeLinks(n, 0, n2);
1875    return nodeToReturn;
1876  }
1877
1878  // -----------------------------------------
1879  // case (3) ????????????????? I am not sure if this is correct implementation !!! VH
1880  // DEBUG << "case 3 \n";
1881 
1882  // make LI part
1883  SBBox lbb = d->box; // the bounding box
1884  lbb.Reduce(axis, 0, d->position2);
1885  int cntLeft = 0;
1886  int cntRight = 1;
1887  d->box = lbb;
1888  d->position = d->position2;
1889  d->axis = axis;
1890  // the next subdivision step for R part
1891  d->pars.reqAxis = axis;
1892  d->pars.reqPosition = d->position;
1893 
1894  IncDepth();
1895  SKTBNodeT *nl = SubDiv(d);
1896  if (!nodeToReturn)
1897    nodeToReturn = nl;
1898 
1899  // make R part
1900  SKTBNodeT *n = MakeOneCut(d);
1901  DecDepth();
1902  SetInteriorNodeLinks(nl, 0, n);
1903  // return left node, linked in DFS order to the previous one
1904  return nodeToReturn;
1905}
1906
1907// makes cutting on the given position on given axis
1908// at value position, bb is the bounding box of current node,
1909// count is the number of objects in the node
1910// flags LeftF and RightF declares if to continue down after the recursion
1911SKTBNodeT *
1912CKTBABuildUp::MakeOneCut(SInputData *d)
1913{
1914// for testing accuracy of setting position from the ideal position
1915#ifdef _RANDOMIZE_POSITION
1916  assert(d->algorithmBreakAx == 0); // position based
1917  RandomizePosition(d->position, d->box.Min(d->axis), d->box.Max(d->axis));
1918#endif //_RANDOMIZE_POSITION
1919
1920#ifdef _DEBUG
1921  assert(d->position >= d->box.Min(d->axis));
1922  assert(d->position <= d->box.Max(d->axis));
1923#endif
1924 
1925  SKTBNodeT *nodeToReturn = 0;
1926
1927#ifdef  _KTB_CONSTR_STATS
1928  // surface area of the interior nodes to be subdivided
1929  _sumSurfaceAreaInteriorNodes += d->box.SA2();
1930#endif // _KTB_CONSTR_STATS
1931 
1932  if (splitClip) {
1933    // empty object list to store the pointers to the split objects
1934    ObjectContainer objlist;
1935    cerr << "Not yet handled" << endl;
1936    abort();;
1937  }
1938
1939#if 0
1940  // check if the lists are correctly sorted .. in debug
1941  static int index = 0;
1942#if 1
1943  cout << "MakeOneCut index = " << index << " .. check axis = "
1944        << d->axis << " count = " << d->count
1945        << " depth = " << GetDepth() << endl;
1946  const int indexToFind = 3743;
1947  if (index == indexToFind) {
1948   
1949    cout << "Debug should start " << endl;
1950    cout << "index = indexToFind" << endl;
1951  }
1952#endif
1953  index++; 
1954#endif
1955 
1956#ifdef _DEBUG
1957  // We do not know the number of boundaries
1958  Check3List(d);
1959#endif
1960 
1961//#define _VYPIS
1962#ifdef _VYPIS
1963  DEBUG << "START: Subdivision at pos="
1964        << position << " axis = " << axis << "\n";
1965  PrintOut(axis, 1, list+axis, sec+axis, Limits::Infinity);
1966#endif
1967
1968  // axis .. the axis perpendicular to splitting plane positioned at position
1969  // break active axis = split the list into two pieces
1970  // list[axis] .. start of the first list in axis direction .. input
1971  // second[axis] .. start of the second list in axis direction .. output
1972  // objlist .. the list of objects that are duplicated
1973#ifdef _DEBUG
1974  Check1List(d, d->axis, d->count);
1975  //cout << "Axis = " << axis << endl;
1976#endif
1977
1978  // We have to allocate a new node for right child
1979  SInputData* rightData = AllocNewData(2*d->count);
1980  // Copy the basic data from parent
1981  rightData->CopyBasicData(d);
1982
1983  // The number of objects for left children and right children
1984  int cntL, cntR;
1985
1986//#ifdef _DEBUG
1987#if 1
1988  if (d->cntThickness < 0) {
1989    cerr << "Problem, d->cntThickness = " << d->cntThickness << endl;
1990    abort();;
1991  }
1992#endif
1993
1994  // Correct for cutting off empty space and many left
1995  // boundaries on the splitting plane.
1996  SItemVec *vec = d->GetItemVec(d->axis);
1997
1998  if (d->algorithmBreakAx) {
1999    // Use pointer-based break-ax algorithm
2000    if (d->bestIterator == vec->begin()) {
2001      assert(d->cntThickness == 0);
2002      d->bestIterator--;
2003    }
2004    else {
2005      if (d->count > 1) {
2006        int origThickness = d->cntThickness;
2007        SItemVec::iterator origIterator = d->bestIterator;
2008        SItemVec::iterator it = d->bestIterator;
2009        int caseBreak = 0;
2010        if ((*(it)).IsLeftBoundary()) {
2011          float pos = d->position;
2012          it--;
2013          int ir = 1;
2014          assert(it != vec->begin()-1);
2015          for ( ; true; ) {
2016            if ( (*(it)).pos < pos) {
2017              caseBreak = 0;
2018              d->bestIterator = it;
2019              if (ir > 1)
2020                d->cntThickness -= ir;
2021              break;
2022            }
2023            if (it == vec->begin()) {
2024              // cutting off, #objects on the left = 0
2025              d->cntThickness = 0;
2026              it--;
2027              d->bestIterator = it;
2028              caseBreak = 1;
2029              break;
2030            }
2031            if ((*(it)).IsRightBoundary()) {
2032              d->bestIterator = it;
2033              if (ir > 1)
2034                d->cntThickness -= ir;
2035              caseBreak = 2;
2036              break;
2037            }
2038            it--;
2039            ir++;
2040          } // for
2041#ifdef _DEBUG
2042          if (d->cntThickness < 0) {
2043            cerr << "Problem, d->cntThickness = " << d->cntThickness << endl;
2044            cerr << " origThickness = " << origThickness << " Depth = "
2045                  << GetDepth() << endl;
2046            for (SItemVec::iterator iq = it; iq != origIterator; iq++) {
2047              cout << "  pos = " << (*iq).pos;
2048              if ((*iq).IsLeftBoundary())
2049                cout << "  L  ";
2050              else
2051                cout << "  R  ";
2052              cout << (*iq).obj << endl;         
2053            } // for
2054            cout << "AAAA" << endl;
2055            abort();;
2056        }
2057#endif
2058        } // if left boundary
2059      } // if d->count > 1
2060    } // if d->vec is not the beginning of the list
2061  } // pointer based break-ax algorithm
2062   
2063#if 0
2064  DEBUG << "MakeOneCut cntL = " << cntL << " cntR = " << cntR << endl;
2065#endif
2066
2067  // Here is the breaking the arrays into two arrays, compute the number
2068  // of objects on the left and on the right of the splitting plane
2069  int dCountTagging = 0;
2070  // Use either
2071  if (d->algorithmBreakAx == 0)
2072    BreakAxPosition(d, d->axis, rightData, cntL, cntR);
2073  else
2074    BreakAx(d, d->axis, rightData, cntL, cntR);
2075   
2076#ifdef _DEBUG
2077  // the number of objects on the left and on the right
2078  int cntL_B = cntL;
2079  int cntR_B = cntR;
2080#if 0
2081  cout << "B - Count*2 = " << d->count*2
2082       << " countSum = " << cntL_B + cntR_B
2083       << " CntL_B = " << cntL_B
2084       << " cntR_B = " << cntR_B << endl;
2085#endif
2086  //Check1List(alist, axis, cntL_B, true, position);
2087  Check1List(d, d->axis, cntL_B);
2088  Check1List(rightData, d->axis, cntR_B);
2089  // cout << "AxisNext = " << oax0 << endl;
2090#endif
2091 
2092  // DEBUG << "THE CHECK at DEPTH=" << GetDepth() << endl;
2093 
2094#ifdef _VYPIS
2095  DEBUG << "After Break at pos="
2096        << position << " axis = " << axis << "\n";
2097  PrintOut(d, axis, d->position);
2098#endif
2099
2100#ifdef _VYPIS
2101  DEBUG << "BEFORE oax0\n";
2102  PrintOut(d, oax0, Limits::Infinity);
2103#endif
2104
2105  // --------------------------------------------------
2106  // break other two axes
2107  CKTBAxes::Axes oax0 = oaxes[d->axis][0];
2108#ifdef _DEBUG
2109  Check1List(d, oax0, d->count);
2110#endif
2111 
2112  SItemVec *vecc = d->GetItemVec(oax0);
2113  if (vecc->size() != 2 * d->count) {
2114    cout << "AAAA 1" << endl;
2115  }
2116 
2117  // Here break the list in the next axis
2118#ifdef _USE_OPTIMIZE_DIVIDE_AX
2119  DivideAx_I_opt(d, oax0, rightData, cntL, cntR);
2120#else
2121  DivideAx_I(d, oax0, rightData, cntL, cntR);
2122#endif
2123
2124  if (vecc->size() != 2 * cntL) {
2125    cout << "AAAA 3" << endl;
2126  }
2127  if (rightData->GetItemVec(oax0)->size() != 2 * cntR) {
2128    cout << "AAAA 4" << endl;
2129    cout << "cntR * 2 = " << 2 * cntR << endl;
2130    cout << "vecSize = " << rightData->GetItemVec(oax0)->size() << endl;
2131  }
2132 
2133#ifdef _DEBUG
2134  // the number of objects on the left and on the right
2135  int cntL_I = cntL;
2136  int cntR_I = cntR;
2137#if 0
2138  cout << "I - Count*2 = " << count*2
2139       << " countSum = " << cntL_I + cntR_I
2140       << " CntL_B = " << cntL_I
2141       << " cntR_B = " << cntR_I << endl;
2142#endif
2143 
2144  if (cntL_I != cntL_B) {
2145    cerr << "FirstList - Problem with DivideAx_I, cntL_B = " << cntL_B
2146          << " and cntL_I = " << cntL_I << endl;
2147  }
2148  if (cntR_I != cntR_B) {
2149    cerr << "SecondList = Problem with DivideAx_I, cntR_B = " << cntR_B
2150          << " and cntR_I = " << cntR_I << endl;
2151  }
2152 
2153  Check1List(d, oax0, cntL_B);
2154  Check1List(rightData, oax0, cntR_B);
2155
2156  //cout << "AxisNextNext = " << oax1 << endl;
2157
2158#endif
2159 
2160#ifdef _VYPIS
2161  DEBUG << "AFTER DivideAx oax0\n";
2162  PrintOut(oax0, 3, alist + oax0, sec+oax0, 0);
2163#endif
2164
2165#ifdef _VYPIS
2166  DEBUG << "BEFORE oax1\n";
2167  PrintOut(oax1, 1, alist + oax1, sec+oax1, Limits::Infinity);
2168#endif
2169
2170  // Her break the list in the next next axis
2171  CKTBAxes::Axes oax1 = oaxes[d->axis][1];
2172
2173  SItemVec *vec2 = d->GetItemVec(oax1);
2174  if (vec2->size() != 2 * d->count) {
2175    cout << "BBBB 1" << endl;
2176  }
2177 
2178#ifdef _USE_OPTIMIZE_DIVIDE_AX
2179  DivideAx_II_opt(d, oax1, rightData, cntL, cntR);
2180#else
2181  DivideAx_II(d, oax1, rightData, cntL, cntR);
2182#endif
2183
2184  if (vec2->size() != 2 * cntL) {
2185    cout << "BBBB 3" << endl;
2186    cout << "cntL * 2 = " << 2 * cntL << endl;
2187    cout << "vec2size = " << vec2->size() << endl;
2188  }
2189  if (rightData->GetItemVec(oax1)->size() != 2 * cntR) {
2190    cout << "BBBB 4" << endl;
2191    cout << "cntR * 2 = " << 2 * cntR << endl;
2192    cout << "vecSize = " << rightData->GetItemVec(oax1)->size() << endl;
2193  }
2194 
2195#ifdef _DEBUG
2196  int cntL_II = cntL;
2197  int cntR_II = cntR;
2198#if 0
2199  cout << "II - Count*2 = " << count*2
2200       << " countSum = " << cntL_II + cntR_II
2201       << " CntL_B = " << cntL_II
2202       << " cntR_B = " << cntR_II << endl;
2203#endif
2204  if (cntL_II != cntL_B) {
2205    cerr << "FirstList - Problem with DivideAx_II, cntL_B = " << cntL_B
2206          << " and cntL_I = " << cntL_I
2207          << " and cntL_II = " << cntL_II << endl;
2208  }
2209  if (cntR_II != cntR_B) {
2210    cerr << "Second List - Problem with DivideAx_II, cntR_B = " << cntR_B
2211          << " and cntR_I = " << cntR_I
2212          << " and cntR_II = " << cntR_II << endl;
2213  }
2214 
2215  Check1List(d, oax1, cntL_B);
2216  Check1List(rightData, oax1, cntR_B);
2217
2218  if ( (cntL_I != cntL_II) ||
2219       (cntR_I != cntR_II) ) {
2220    cout << "Problem with dividing axis II: cntL_I = " << cntL_I
2221         << " cntL_II = " << cntL_II
2222         << " cntR_I = " << cntR_I
2223         << " cntR_II = " << cntR_II << endl;     
2224  }
2225#endif
2226 
2227#ifdef _VYPIS
2228  DEBUG << "AFTER DivideAx oax1\n";
2229  PrintOut(oax1, 3, alist + oax1, sec+oax1, position);
2230#endif
2231 
2232  // Allocate the representation of the interior node
2233  SKTBNodeT *node = 0;
2234
2235  if (makeMinBoxHere) {
2236    // -------------------------------------------------------
2237    // interior node with the box
2238    node = AllocInteriorNodeWithBox(// interior node data
2239                                    d->axis, d->position, cntL, cntR,
2240                                    // box data
2241                                    d->box, d->lastMinBoxNode, d->lastDepthForMinBoxes);
2242
2243    //cout << "depth = " << GetDepth() << " node = " << node
2244    //   << " lastMinBoxNode = " << d->lastMinBoxNode << endl;
2245   
2246    d->lastMinBoxNode = node; // we have to remember the last node
2247    rightData->lastMinBoxNode = node;
2248    makeMinBoxHere = false; // reset for the next time
2249  }
2250  else
2251    // normal interior node (=without box)
2252    node = AllocInteriorNode(d->axis, d->position, cntL, cntR);
2253 
2254#ifdef _DEBUG
2255  if (!node) {
2256    cerr << "Allocation of Interior Node failed.\n";
2257    abort();;
2258  }
2259#endif // _DEBUG
2260  // set the node to be returned
2261  if (!nodeToReturn)
2262    nodeToReturn = nodeToLink;
2263  assert(node);
2264  assert(nodeToReturn);
2265
2266#if 0 
2267  // If we make split clipping, then we want to reduce the boxes for
2268  // the objects straddling the splitting plane now.
2269  if ((splitClip) && (objlist.ListCount() > 0 )) {
2270    // if the objects are split by splitting plane
2271    //DEBUG << "Reduces box ax=" << axis << " position="
2272    // << position << "#split_objs=" << objlist.ListCount() << "\n" << flush;
2273    ReduceBBoxes(d, d->axis, rightData, d->position);
2274  }
2275#endif
2276 
2277  // ---------------------------------------------
2278  // initialize the left bounding box 'bb'
2279  d->box.Reduce(d->axis, 0, d->position);
2280
2281  // initialize the right bounding box 'rbb'
2282  rightData->box.Reduce(d->axis, 1, d->position);
2283
2284  // Set correct number of objects
2285  d->count = cntL;
2286  rightData->count = cntR;
2287 
2288  // The children nodes to be created
2289  SKTBNodeT *nodeL=0, *nodeR=0;
2290
2291  // Recurse to the left child
2292  if (d->makeSubdivisionLeft) {
2293    // subdivide branches
2294    ESubdivMode modeLeft = d->modeSubDiv ;
2295    if ( (cutEmptySpace) &&
2296         (modeLeft != EE_SUBDIVCUTEMPTY) &&
2297         (d->count <= maxListLength) ) {
2298      // setting initial depth if empty space cutting
2299      startEmptyCutDepth = GetDepth();
2300      d->modeSubDiv = EE_SUBDIVCUTEMPTY;
2301    }
2302    // DEBUG << "Left cnt= " << (cnt[0] >> 1) << " depth= "
2303    //       << GetDepth() << "\n";
2304    nodeL = SubDiv(d);
2305  }
2306 
2307  // Recurse to the right child
2308  if (rightData->makeSubdivisionRight) {
2309    ESubdivMode modeRight = d->modeSubDiv ;
2310    if ( (cutEmptySpace) &&
2311         (modeRight != EE_SUBDIVCUTEMPTY) &&
2312         (rightData->count <= maxListLength) ) {
2313      // setting initial depth if empty space cutting
2314      startEmptyCutDepth = GetDepth();
2315      rightData->modeSubDiv = EE_SUBDIVCUTEMPTY;
2316    }
2317    // DEBUG << "Right cnt= " << (cnt[1] >> 1) << " depth= "
2318    // << GetDepth() << "\n";
2319    nodeR = SubDiv(rightData);
2320
2321  }
2322
2323  // Free the data to be reused further on
2324  FreeLastData();
2325
2326  // Set the node pointers to the children for created node
2327  SetInteriorNodeLinks(node, nodeL, nodeR);
2328
2329  return nodeToReturn; // return the node
2330}
2331
2332// returns a box enclosing all the objects in the node
2333void
2334CKTBABuildUp::GetTightBox(const SInputData &i, SBBox &tbox)
2335{
2336  // Index to the last boundary
2337  int cnt = 2*i.count-1; 
2338
2339  tbox.MinX() = (*(i.xvec))[0].pos;
2340  tbox.MinY() = (*(i.yvec))[0].pos;
2341  tbox.MinZ() = (*(i.zvec))[0].pos;
2342  tbox.MaxX() = (*(i.xvec))[cnt].pos;
2343  tbox.MaxY() = (*(i.yvec))[cnt].pos;
2344  tbox.MaxZ() = (*(i.zvec))[cnt].pos;
2345
2346  assert(tbox.MinX() >= i.box.MinX());
2347  assert(tbox.MaxX() <= i.box.MaxX());
2348  assert(tbox.MinY() >= i.box.MinY());
2349  assert(tbox.MaxY() <= i.box.MaxY());
2350  assert(tbox.MinZ() >= i.box.MinZ());
2351  assert(tbox.MaxZ() <= i.box.MaxZ());
2352}
2353
2354
2355// returns a box enclosing all the objects in the node
2356int
2357CKTBABuildUp::GetEBox(const SInputData &i, SBBox &tbox)
2358{
2359  // LEFT BOX
2360  // initialization
2361
2362  // Index to the last boundary
2363  int cnt = 2*i.count-1; 
2364  int changed = 0; // number of planes, that are changed by bounding box
2365
2366  float xMin = (*(i.xvec))[0].pos;
2367  changed = (xMin > i.box.Min(0)) ? changed+1 : changed;   
2368  float xMax = (*(i.xvec))[cnt].pos;
2369  changed = (xMax < i.box.Max(0)) ? changed+1 : changed;   
2370
2371  float yMin = (*(i.yvec))[0].pos;
2372  changed = (yMin > i.box.Min(1)) ? changed+1 : changed;   
2373  float yMax = (*(i.yvec))[cnt].pos;
2374  changed = (yMax < i.box.Max(1)) ? changed+1 : changed;   
2375
2376  float zMin = (*(i.zvec))[0].pos;
2377  changed = (zMin > i.box.Min(2)) ? changed+1 : changed;   
2378  float zMax = (*(i.zvec))[cnt].pos;
2379  changed = (zMax < i.box.Max(2)) ? changed+1 : changed;   
2380
2381  // Set the bounding (reduced) box
2382  tbox.Min() = Vector3(xMin, yMin, zMin);
2383  tbox.Max() = Vector3(xMax, yMax, zMax);
2384
2385  // Return the number of bounding splitting planes with respect to original box
2386  return changed;
2387}
2388
2389
2390// This is only for debugging purposes, to be removed!
2391static
2392CKTBABuildUp::SSolid* objtf = (CKTBABuildUp::SSolid*)0x0;
2393
2394// removes the objects specified in "tagobjlist" from boundary "list"
2395// it supposes the "tagobjlist" is ordered in the left boundary order
2396// in CKTBAxes::EE_X_axis.
2397void
2398CKTBABuildUp::RemoveObjects(SItemVec *vec, int cntObjects)
2399{
2400  assert(vec);
2401  assert(cntObjects);
2402 
2403  // number of removed left boundaries
2404  int cnt = 0;
2405
2406  SItemVec::iterator curr;
2407#ifdef _DEBUG
2408  int cnt2 = 0;
2409  for (curr = vec->begin(); curr != vec->end(); curr++)
2410  {
2411    if (curr->obj->ToBeRemoved()) {
2412      //cout << cnt2 << "  " << curr - vec->begin() << endl;
2413      cnt2++;
2414    }
2415    //if (curr->obj == objtf)
2416    //  cout << "Index found cnt2 = " << cnt2 << endl;
2417  } // for
2418  if (cnt2 != 2*cntObjects) {
2419    cerr << "Some problem with data marking" << endl;
2420    abort();;
2421  }
2422#endif
2423
2424  // Find the first boundary
2425  for (curr = vec->begin(); curr != vec->end(); curr++)
2426  {
2427    //if (curr->obj == objtf)
2428    //cout << "2 Index found cnt2 = " << cnt2 << endl;
2429
2430    if (curr->obj->ToBeRemoved()) {
2431      cnt++; // the number of removed boundaries
2432      assert(curr->IsLeftBoundary());
2433      //cout << cnt << "  objToRemove=" << curr->obj << endl;
2434      break;
2435    }
2436  }
2437  // the first element to be overwritten
2438  SItemVec::iterator it = curr;
2439  curr++;
2440
2441  for ( ; curr != vec->end(); curr++) {
2442    //if (curr->obj == objtf) {
2443    //  cout << "3 Index found cnt2 = " << cnt2 << endl;
2444
2445    if (curr->obj->ToBeRemoved()) {
2446      cnt++; // the number of removed boundaries
2447      //cout << cnt << "  " << curr->obj << endl;
2448      if (cnt == cntObjects*2) {
2449        curr++;
2450        break;
2451      }
2452    }
2453    else {
2454      // copy the element to the right place
2455      *it = *curr;
2456      it++;
2457    }
2458  } // for
2459
2460#ifdef _DEBUG
2461  if (cnt != 2*cntObjects) {
2462    cerr << "Problem with removing objects implementation" << endl;
2463    abort();;
2464  }
2465#endif
2466
2467  // copy the rest of the list to the right place
2468  for ( ; curr != vec->end(); curr++, it++)
2469  {
2470    //if (curr->obj == objtf) {
2471    //  cout << "4 Index found cnt2 = " << cnt2 << endl;
2472
2473    // copy the element to the right place
2474    *it = *curr;
2475  }
2476     
2477  // resize the vector correctly
2478  int oldSize = vec->size();
2479  assert(oldSize >= cnt);
2480  if (oldSize == cnt)
2481    vec->erase(vec->begin(), vec->end());
2482  else
2483    vec->resize(oldSize - cnt);
2484  assert(vec->size() % 2 == 0);
2485  assert(vec->size() == oldSize - cnt);
2486
2487#ifdef _DEBUG   
2488  for (curr = vec->begin(); curr != vec->end(); curr++)
2489  {
2490    if (curr->obj->ToBeRemoved()) {
2491      cerr << "Implementation bug" << endl;
2492      abort();;
2493    } // if left boundary
2494  } // for
2495#endif
2496
2497  //cout << "size = " << vec->size() << endl;
2498
2499  return;
2500}
2501
2502// removes the objects specified in "tagobjlist" from boundary "list"
2503// it supposes the "tagobjlist" is ordered in the left boundary order
2504// in CKTBAxes::EE_X_axis.
2505void
2506CKTBABuildUp::RemoveObjectsReset(SItemVec *vec, int cntObjects)
2507{
2508  assert(vec);
2509  assert(cntObjects);
2510 
2511  // number of removed left boundaries
2512  int cnt = 0;
2513  SItemVec::iterator curr;
2514
2515#ifdef _DEBUG
2516  int cnt2 = 0;
2517  for (curr = vec->begin(); curr != vec->end(); curr++)
2518  {
2519    if (curr->obj->ToBeRemoved()) {
2520      //cout << cnt2 << "  " << curr - vec->begin() << endl;
2521      cnt2++;
2522    }
2523  } // for
2524  if (cnt2 != 2*cntObjects) {
2525    cerr << "Some problem with data marking" << endl;
2526    abort();;
2527  }
2528#endif
2529
2530  // Find the first boundary
2531  for (curr = vec->begin(); curr != vec->end(); curr++)
2532  {
2533    if (curr->obj->ToBeRemoved()) {
2534      cnt++; // the number of removed boundaries
2535      assert(curr->IsLeftBoundary());
2536      curr++;
2537      break;
2538    }
2539  }
2540  SItemVec::iterator it = curr;
2541  it--; // the first element to be overwritten
2542
2543  for ( ; curr != vec->end(); curr++) {
2544    if (curr->obj->ToBeRemoved()) {
2545      cnt++; // the number of removed boundaries
2546      if (curr->IsRightBoundary()) {
2547        curr->obj->ResetFlags();
2548        assert(curr->obj->flags == 0);
2549      }
2550      if (cnt == cntObjects*2) {
2551        curr++;
2552        break;
2553      }
2554    }
2555    else {
2556      // copy the element to the right place
2557      *it = *curr;
2558      it++;
2559    }
2560  } // for
2561 
2562#ifdef _DEBUG
2563  if (cnt != 2*cntObjects) {
2564    cerr << "Problem with removing objects implementation" << endl;
2565    abort();;
2566  }
2567#endif
2568
2569  // copy the rest of the list to the right place
2570  for ( ; curr != vec->end(); curr++, it++)
2571  {
2572    // copy the element to the right place
2573    *it = *curr;
2574  }
2575     
2576  // resize the vector correctly
2577  int oldSize = vec->size();
2578  assert(oldSize >= cnt);
2579  if (oldSize == cnt)
2580    vec->erase(vec->begin(), vec->end());
2581  else
2582    vec->resize(oldSize - cnt);
2583  assert(vec->size() % 2 == 0);
2584  assert(vec->size() == oldSize - cnt);
2585
2586#ifdef _DEBUG   
2587  for (curr = vec->begin(); curr != vec->end(); curr++)
2588  {
2589    if (curr->obj->ToBeRemoved()) {
2590      cerr << "Implementation bug" << endl;
2591      abort();;
2592    } // if left boundary
2593  } // for
2594#endif
2595
2596  //cout << "size = " << vec->size() << endl;
2597 
2598  return;
2599}
2600
2601
2602// make the full leaf from current node
2603SKTBNodeT*
2604CKTBABuildUp::MakeLeaf(SInputData *d)
2605{
2606#ifdef  _KTBPRINTCUTS
2607  cout << " Creating leaf, depth = " << GetDepth()
2608       << " cntObjs = " << count << endl;
2609  //exit(3);
2610#endif
2611 
2612#ifdef _KTB_CONSTR_STATS
2613  // update the statistics
2614  float sa2bb;
2615  _sumSurfaceAreaLeaves += (sa2bb = d->box.SA2());
2616  _sumSurfaceAreaMULcntLeaves += ((float)d->count * sa2bb);
2617#endif // _KTB_CONSTR_STATS
2618 
2619  // ------------------------------------------------------
2620  // The version with allocating empty leaves. This makes
2621  // sense since in DFS order it works correctly
2622  SKTBNodeT *node = AllocLeaf(d->count);
2623
2624  if (d->count > 0) {
2625    // copy the list of objects
2626    ObjectContainer *objlist = NULL;
2627    objlist = new ObjectContainer;
2628
2629    SItemVec *vec = d->GetItemVec(0);
2630    SItemVec::iterator curr;   
2631    for (curr = vec->begin(); curr != vec->end(); curr++) {
2632      if ( (*curr).IsLeftBoundary()) // skip all right boundaries
2633        objlist->push_back( (*curr).obj->obj);
2634    }
2635
2636    assert(objlist->size() != 0);
2637    // Do not delete the object list, it is used as allocated above
2638    SetFullLeaf(node, objlist);
2639  }
2640
2641  // We assume that the allocation of the single leaf cannot
2642  // fail, since this one shall be the last in the sequence.
2643  assert(node == nodeToLink);
2644 
2645  // Return the node to be used in the linking
2646  return nodeToLink;
2647}
2648
2649// split the list along the required axis
2650// first is the begin of the SItem list, axis is the axis, where
2651// to split, val is the splitting value in integer notation.
2652// The output of the function are two lists, first and second
2653// Optionally, objlist store the objects that straddle the splitting plane
2654void
2655CKTBABuildUp::BreakAx(SInputData *d, int axis,
2656                      SInputData *rightData,
2657                      int &cntLout, int &cntRout)
2658{
2659  // the number of object boundaries for left array
2660  int cntLeft = (d->bestIterator) - d->GetItemVec(axis)->begin() + 1;
2661  assert(cntLeft >= 0);
2662  assert(cntLeft <= 2*d->count); 
2663  // the number of object boundaries for right array
2664  int cntRight = d->count*2 - cntLeft;
2665  assert(cntRight >= 0);
2666  assert(cntRight <= 2*d->count);
2667  assert(d->cntThickness >= 0);
2668  cntLeft += d->cntThickness; // duplicated boundaries for left array
2669  cntRight += d->cntThickness; // duplicated boundaries for right array
2670  assert(cntLeft + cntRight == (d->count + d->cntThickness) * 2);
2671
2672#ifdef _DEBUG
2673  if ( (cntLeft % 2 != 0) ||
2674       (cntRight % 2 != 0) ) {
2675    int i = 0;
2676    SItemVec *vec = d->GetItemVec(axis);
2677    cout << "d->count = " << d->count ;
2678    cout << "cntLeft = " << cntLeft << " cntRight" << cntRight << endl;
2679    cout << "bestIterator = " << (int)(d->bestIterator - vec->begin());
2680    cout << " cntThickness= " << d->cntThickness << " position=" << d->position << endl;
2681    for (SItemVec::iterator itt = vec->begin(); itt != vec->end(); itt++, i++)
2682    {
2683      cout << i << " = ";
2684      if ((*itt).IsLeftBoundary())
2685        cout << "L"; else cout << "R";
2686      cout << " obj = " << (*itt).obj << " pos = " << (*itt).pos << endl;
2687    }
2688    abort();;
2689  }
2690#endif
2691 
2692  // A special case, right list is empty
2693  if (cntRight == 0) {
2694    assert(d->cntThickness == 0);
2695    rightData->Alloc(0);
2696    rightData->count = 0;
2697    // The final setting of number of objects
2698    cntLout = cntLeft / 2;
2699    cntRout = cntRight / 2;   
2700    return; // nothing to do
2701  }
2702
2703  // We have some problem that is difficult to detect, when tagging objects
2704  // are allowed, the flags are not always reset to zero in some cases.
2705  // VH 2/1/2006. The flags are always left boundaries (=1).
2706  // Therefore here it is a hack to assure that the flags are definitely zero.
2707  if (resetFlagsForBreakAx) {
2708    SItemVec *vecd = d->GetItemVec(axis);
2709    for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
2710      (*itt).obj->ResetFlags();
2711    } // for
2712  } // if
2713 
2714#ifdef _DEBUG
2715  SItemVec *vecd = d->GetItemVec(axis);
2716  for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
2717    if ((*itt).obj->Flags() != 0) {
2718      cout << "Init ERROR 1 in original list, flags = "
2719           << (*itt).obj->Flags() << " obj = " << (void*)((*itt).obj) << endl;
2720    }
2721  } // for
2722#endif // _DEBUG 
2723
2724  assert(cntRight > 0);
2725 
2726  // Allocate the appropriate size for right array, for all 3 axes
2727  rightData->Alloc(cntRight);
2728 
2729  // First find create all left boundaries for right array
2730  SItem boundary; // create left boundary
2731  SItemVec *vecRight = rightData->GetItemVec(axis);
2732
2733#define _USE_PUSHBACK
2734#ifdef _USE_PUSHBACK
2735  vecRight->resize(0);
2736#else
2737  vecRight->resize(cntRight);
2738#endif 
2739  if (d->cntThickness) {
2740    boundary.pos = d->position;
2741    boundary.axis = axis;
2742    boundary.SetLeftBoundary();
2743#ifdef _USE_PUSHBACK
2744    int i;
2745    for (i = 0; i < d->cntThickness; i++) {
2746      vecRight->push_back(boundary); // Note that the objects pointed to are not set !!!
2747    } // for i;
2748#else
2749    for (int i = 0; i < d->cntThickness; i++) {
2750      (*vecRight)[i] = boundary; // Note that the objects pointed to are not set !!!
2751    } // for i;
2752#endif
2753  }
2754  const SItemVec::iterator itendLeft = d->bestIterator + 1;
2755  SItemVec *vec = d->GetItemVec(axis);
2756  SItemVec::iterator it;
2757  //int cntSetLeft = 0;
2758  for (it = vec->begin(); it != itendLeft; it++) {
2759    // Mark all the items to be in the first list
2760    (*it).obj->SetInFirstList();
2761    //cntSetLeft++;
2762  } // for
2763  // The boundaries that belong definitely to the second list
2764  //int cntSetRight = 0;
2765#ifdef _USE_PUSHBACK
2766  const SItemVec::iterator itend = vec->end();
2767  for (; it != itend; it++) {
2768    // Mark all the items to be in the second list
2769    (*it).obj->SetInSecondList();
2770    // and copy the item to the right array
2771    vecRight->push_back(*it);
2772    //cntSetRight++;
2773  } // for
2774#else
2775  const int imax = vec->end() - (d->bestIterator + 1);
2776  int ir = d->cntThickness;
2777  for (int i = 0; i < imax; i++) {
2778    (*it).obj->SetInSecondList();
2779    // and copy the item to the right array
2780    (*vecRight)[ir] = (*it);
2781    ir++; it++;
2782  }
2783  assert(ir == cntRight);
2784#endif
2785 
2786#ifdef _DEBUG
2787  if (vecRight->size() != cntRight) {
2788    cerr << "Implementation problem vecRight->size() = "
2789          << vecRight->size() << "  cntRight = " << cntRight
2790          << " d->cntThickness = " << d->cntThickness << endl;
2791    abort();;
2792  }
2793#endif 
2794  // Only check, if right boundary is set correctly
2795  assert(vecRight->size() == cntRight);
2796 
2797  if (d->cntThickness) {
2798    // Now go only through the items on the left side of the splitting plane
2799    SItemVec::iterator itR = vecRight->begin(); // to set left boundaries in right array
2800    SItemVec::iterator itL = itendLeft; // to set new right boundaries in left array
2801    boundary.SetRightBoundary();
2802#ifdef _USE_PUSHBACK
2803    for (it = vec->begin(); it != itendLeft; it++) {
2804      if ((*it).obj->InBothLists()) {
2805        // set new right boundaries in left array
2806        (*itL) = boundary;
2807        (*itL).obj = (*it).obj; // set the right object
2808        itL++;
2809        // and set new left boundary in right array, only object is set
2810        (*itR).obj = (*it).obj;
2811        itR++;
2812      }
2813    } // for
2814#else
2815    int imax = itendLeft - vec->begin();
2816    int i;
2817    for (i = 0; i < imax ; i++) {
2818      if ((*vec)[i].obj->InBothLists()) {
2819        // set new right boundaries in left array
2820        (*itL) = boundary;
2821        (*itL).obj = (*vec)[i].obj; // set the right object
2822        itL++;
2823        // and set new left boundary in right array, only object is set
2824        (*itR).obj = (*vec)[i].obj;
2825        itR++;
2826      }
2827    } // for
2828#endif
2829#ifdef _DEBUG
2830    int sizeL = itL - vec->begin();
2831    assert(sizeL == cntLeft);
2832    int sizeR = itR - vecRight->begin();
2833    assert(sizeR == d->cntThickness);
2834#endif
2835  }
2836#ifdef _DEBUG
2837  else {
2838    vec->resize(cntLeft);   
2839    assert(d->cntThickness == 0);
2840    int cntLL = 0;
2841    int cntRR = 0;
2842    for (it = vec->begin(); it != vec->end(); it++) {
2843      if ((*it).obj->InBothLists()) {
2844        cntLL++;
2845        cout << "ERROR in left list " << cntLL << endl;
2846      }
2847    } // for   
2848    for (it = vecRight->begin(); it != vecRight->end(); it++) {
2849      if ((*it).obj->InBothLists()) {
2850        cntRR++;
2851        cout << "ERROR in right list " << cntRR << endl;
2852      }
2853    } // for   
2854  } 
2855#endif // _DEBUG 
2856 
2857  // The boundary on the left, unused items are not used
2858  vec->resize(cntLeft);
2859
2860  // The final setting of size
2861  cntLout = cntLeft / 2;
2862  cntRout = cntRight / 2;
2863
2864  return;
2865}
2866
2867
2868void
2869CKTBABuildUp::BreakAxPosition(SInputData *d, int axis,
2870                              SInputData *rightData,
2871                              int &cntLout, int &cntRout)
2872{
2873  // the number of object boundaries for left array is not known
2874  int cntLeft = 0;
2875  int cntRight = 0;
2876  d->cntThickness = 0;
2877  const float positionToDecide = d->position;
2878
2879  // We have some problem that is difficult to detect, when tagging objects
2880  // are allowed, the flags are not always reset to zero in some cases.
2881  // VH 2/1/2006. The flags are always left boundaries (=1).
2882  // Therefore here it is a hack to assure that the flags are definitely zero.
2883  if (resetFlagsForBreakAx) {
2884    SItemVec *vecd = d->GetItemVec(axis);
2885    for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
2886      (*itt).obj->ResetFlags();
2887    } // for
2888  } // if
2889 
2890#ifdef _DEBUG
2891  SItemVec *vecd = d->GetItemVec(axis);
2892  for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
2893    if ((*itt).obj->Flags() != 0) {
2894      cout << "Init ERROR 2 in original list, flags = "
2895           << (*itt).obj->Flags() << " obj = " << (void*)((*itt).obj) << endl;
2896    }
2897  } // for
2898#endif // _DEBUG 
2899
2900  // Allocate the appropriate size for right array, for all 3 axes
2901  // This is only estimated at this algorithm.
2902  int estimatedSize = (int)(0.6f * (float)d->count);
2903  if (estimatedSize < 10)
2904    estimatedSize = 10;
2905  rightData->Alloc(estimatedSize * 2);
2906 
2907  // First find create all left boundaries for right array
2908  SItemVec *vecRight = rightData->GetItemVec(axis);
2909
2910  // Resize to 0 to start from the beginning
2911  vecRight->resize(0);
2912  SItemVec *vec = d->GetItemVec(axis);
2913  SItemVec::iterator it;
2914  const SItemVec::iterator itEnd = vec->end();
2915  for (it = vec->begin(); it != itEnd; it++) {
2916    const float &pos = (*it).pos; 
2917    if (pos < positionToDecide) {
2918      cntLeft++;
2919      // Mark the item to be in the first list
2920      (*it).obj->SetInFirstList();
2921    }
2922    else {
2923      if (pos > positionToDecide) {
2924        // It belongs to the right list
2925        break;
2926      }
2927      else {
2928        // pos == positionToDecide
2929        if ((*it).IsRightBoundary()) {
2930          // it belongs to the left list
2931          cntLeft++;
2932          (*it).obj->SetInFirstList();
2933        }
2934        else {
2935          // this is the left boundary
2936          // it belongs already to the right list
2937          break;
2938        }
2939      }
2940    }
2941  } // for it
2942
2943  // Remember the position of the iterator, where right
2944  // list has to start
2945  SItemVec::iterator itStart = it;
2946  int cntDups = 0;
2947  SItem boundary;
2948
2949  // create left boundaries for right list
2950  boundary.pos = d->position;
2951  boundary.axis = axis;
2952  boundary.SetLeftBoundary();
2953  for (; it != itEnd; it++) {
2954    assert((*it).pos >= positionToDecide);
2955    // Mark all the items to be in the second list
2956    if ((*it).obj->InFirstList()) {
2957      // and copy the newly created boundary to the right array
2958      boundary.obj = (*it).obj;
2959      vecRight->push_back(boundary);
2960      cntDups++;
2961      cntRight++;
2962    }
2963  } // for
2964
2965  // Now go through the right part of the list again
2966  // and copy the rest of the list to the right list
2967  for (it = itStart; it != itEnd; it++) {
2968    (*it).obj->SetInSecondList();
2969    assert((*it).pos >= positionToDecide);
2970    vecRight->push_back((*it));
2971    cntRight++;
2972  } // for
2973
2974  // in final pass create the new boundaries for the left list
2975  int i;
2976  SItemVec::iterator srcIt = vecRight->begin();
2977  for (it = itStart, i = 0 ; i < cntDups; it++, srcIt++, i++) {
2978    cntLeft++;
2979    // first copy the boundary from the beginning of the right list
2980    *it = *srcIt;
2981    // then set the right boundary in copied item
2982    (*it).SetRightBoundary();
2983    assert( (*srcIt).IsLeftBoundary());
2984  } // for it
2985 
2986  // The boundary on the left, unused items are not used
2987  vec->resize(cntLeft);
2988  assert(cntLeft % 2 == 0);
2989  assert(cntRight % 2 == 0);
2990  assert(cntLeft+cntRight == 2*(d->count+cntDups));
2991 
2992  // The final setting of size
2993  cntLout = cntLeft / 2;
2994  cntRout = cntRight / 2;
2995  // Set correctly also the number of objects
2996  d->cntThickness = cntDups;
2997
2998#ifdef _DEBUG
2999  SItemVec *vecdL = d->GetItemVec(axis);
3000  int cntL = 0; int cntL_LB = 0; int cntL_RB = 0; int cntL_LR = 0;
3001  for (SItemVec::iterator itl = vecdL->begin(); itl != vecdL->end(); itl++) {
3002    if ((*itl).obj->InBothLists())
3003      cntL_LR++;
3004    if ((*itl).IsLeftBoundary()) cntL_LB++;
3005    if ((*itl).IsRightBoundary()) cntL_RB++;
3006    cntL++;
3007  } // for
3008  SItemVec *vecdR = rightData->GetItemVec(axis);
3009  int cntR = 0; int cntR_LB = 0; int cntR_RB = 0; int cntR_LR = 0;
3010  for (SItemVec::iterator itr = vecdR->begin(); itr != vecdR->end(); itr++) {
3011    if ((*itr).obj->InBothLists())
3012      cntR_LR++;
3013    if ((*itr).IsLeftBoundary()) cntR_LB++;
3014    if ((*itr).IsRightBoundary()) cntR_RB++;
3015    cntR++;
3016  } // for
3017  assert(cntL_LR == cntR_LR);
3018  assert(cntL_LB == cntL_RB);
3019  assert(cntR_LB == cntR_RB);
3020#endif
3021 
3022  return;
3023}
3024
3025#if 0
3026// This below was an attempt to do it differently in Vienna, Dec 2007.
3027// VH
3028
3029// ---------------------------------------------------------
3030// split the list along the required axis
3031// first is the begin of the SItem list, axis is the axis, where
3032// to split, val is the splitting value in integer notation.
3033// The output of the function are two lists, first and second
3034// Optionally, objlist store the objects that straddle the splitting plane
3035void
3036CKTBABuildUp::BreakAxPosition(SInputData *d, int axis,
3037                              SInputData *rightData,
3038                              int &cntLout, int &cntRout)
3039{
3040  // the number of object boundaries for left array is not known
3041  int cntLeft = 0;
3042  int cntRight = 0;
3043  d->cntThickness = 0;
3044  const float positionToDecide = d->position;
3045
3046  // We have some problem that is difficult to detect, when tagging objects
3047  // are allowed, the flags are not always reset to zero in some cases.
3048  // VH 2/1/2006. The flags are always left boundaries (=1).
3049  // Therefore here it is a hack to assure that the flags are definitely zero.
3050  if (resetFlagsForBreakAx) {
3051    SItemVec *vecd = d->GetItemVec(axis);
3052    for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
3053      (*itt).obj->ResetFlags();
3054    } // for
3055  } // if
3056 
3057#ifdef _DEBUG
3058  SItemVec *vecd = d->GetItemVec(axis);
3059  for (SItemVec::iterator itt = vecd->begin(); itt != vecd->end(); itt++) {
3060    if ((*itt).obj->Flags() != 0) {
3061      cout << "Init ERROR 2 in original list, flags = "
3062           << (*itt).obj->Flags() << " obj = " << (void*)((*itt).obj) << endl;
3063    }
3064  } // for
3065#endif // _DEBUG 
3066
3067  // Allocate the appropriate size for right array, for all 3 axes
3068  // This is only estimated at this algorithm.
3069  int estimatedSize = (int)(0.6f * (float)d->count);
3070  if (estimatedSize < 10)
3071    estimatedSize = 10;
3072  assert(estimatedSize > 0);
3073  rightData->Alloc(estimatedSize * 2);
3074 
3075  // First find create all left boundaries for right array
3076  SItemVec *vecRight = rightData->GetItemVec(axis);
3077
3078  // Resize to 0 to start from the beginning
3079  vecRight->resize(0);
3080  SItemVec *vec = d->GetItemVec(axis);
3081  SItemVec::iterator it;
3082  SItemVec::iterator itLeft;
3083  const SItemVec::iterator itEnd = vec->end();
3084  bool foundBoundary = false; 
3085  for (it = vec->begin(); it != itEnd; it++) {
3086    const float &pos = (*it).pos; 
3087    if (pos < positionToDecide) {
3088      cntLeft++;
3089      // Mark the item to be in the first list
3090      (*it).obj->SetInFirstList();
3091    }
3092    else {
3093      if (pos > positionToDecide) {
3094        cntRight++;
3095        (*it).obj->SetInSecondList();
3096        if (!foundBoundary) {
3097          foundBoundary = true;
3098          itLeft = it;
3099        }
3100        // It belongs to the right list
3101      }
3102      else {
3103        // pos == positionToDecide
3104        if ((*it).IsRightBoundary()) {
3105          // it belongs to the left list
3106          cntLeft++;
3107          (*it).obj->SetInFirstList();
3108        }
3109        else {
3110          // this is the left boundary
3111          // it belongs already to the right list
3112          cntRight++;
3113          (*it).obj->SetInSecondList();
3114          if (!foundBoundary) {
3115            foundBoundary = true;
3116            itLeft = it;
3117          }
3118        }
3119      }
3120    }
3121  } // for it
3122
3123  // Remember the position of the iterator, where right
3124  // list has to start
3125  SItem boundary;
3126  int cntDups = 0;
3127  cntLeft = 0;
3128  cntRight = 0;
3129 
3130  // create left boundaries for right list
3131  boundary.pos = d->position;
3132  boundary.axis = axis;
3133  boundary.SetLeftBoundary();
3134  for (it = itLeft; it != itEnd; it++) {
3135    assert((*it).pos >= positionToDecide);
3136    // Mark all the items to be in the second list
3137    if ((*it).obj->InFirstList()) {
3138      // and copy the newly created boundary to the right array
3139      boundary.obj = (*it).obj;
3140      vecRight->push_back(boundary);
3141      cntDups++;
3142      cntRight++;
3143    }
3144  } // for
3145
3146  // Now go through the right part of the list again
3147  // and copy the rest of the list to the right list
3148  for (it = vec->begin(); it != itEnd; it++) {
3149    const float &pos = (*it).pos;
3150    if (pos < positionToDecide) {
3151      // we copy the boundary to the left list
3152      *itLeft = *it;
3153      itLeft++;
3154      cntLeft++;
3155    }   
3156    if (pos > positionToDecide) {
3157      // we copy the boundary to the right list
3158      vecRight->push_back(*it);
3159      cntRight++;
3160    }
3161    else {
3162      if (pos == positionToDecide) {
3163        if ((*it).obj->InBothLists()) {
3164          if ((*it).IsRightBoundary()) {
3165            // right boundary is added to the right as
3166            // we already created left boundaries
3167            vecRight->push_back(*it);
3168            cntRight++;
3169          }
3170          else {
3171            // we copy the left boundary to the left list
3172            *itLeft = *it;
3173            itLeft++;
3174            cntLeft++;
3175          }
3176        }
3177        else {
3178          // This boundary belongs to only a single list
3179          if ((*it).obj->InFirstList()) {
3180            // we copy this boundary  to the left list
3181            *itLeft = *it;
3182            itLeft++;
3183            cntLeft++;
3184          }
3185          else
3186          if ((*it).obj->InSecondList()) {
3187            // we copy this boundary to the right list
3188            vecRight->push_back(*it);
3189            cntRight++;     
3190          }
3191        }
3192      }
3193    }
3194  } // for
3195
3196 
3197  // in final pass create the new boundaries for the left list
3198  int i;
3199  SItemVec::iterator srcIt = vecRight->begin();
3200  for (it = itLeft, i = 0 ; i < cntDups; it++, srcIt++, i++) {
3201    cntLeft++;
3202    // first copy the boundary from the beginning of the right list
3203    assert( (*srcIt).IsLeftBoundary());
3204    *it = *srcIt;
3205    // then set the right boundary in copied item
3206    (*it).SetRightBoundary();
3207    assert( (*it).IsRightBoundary());
3208  } // for it
3209 
3210  // The boundary on the left, unused items are not used
3211  vec->resize(cntLeft);
3212  assert(cntLeft % 2 == 0);
3213  assert(cntRight % 2 == 0);
3214  assert(cntLeft+cntRight == 2*(d->count+cntDups));
3215 
3216  // The final setting of size
3217  cntLout = cntLeft / 2;
3218  cntRout = cntRight / 2;
3219  // Set correctly also the number of objects
3220  d->cntThickness = cntDups;
3221
3222#ifdef _DEBUG
3223  SItemVec *vecdL = d->GetItemVec(axis);
3224  int cntL = 0; int cntL_LB = 0; int cntL_RB = 0; int cntL_LR = 0;
3225  for (SItemVec::iterator itl = vecdL->begin(); itl != vecdL->end(); itl++) {
3226    if ((*itl).obj->InBothLists())
3227      cntL_LR++;
3228    if ((*itl).IsLeftBoundary()) cntL_LB++;
3229    if ((*itl).IsRightBoundary()) cntL_RB++;
3230    cntL++;
3231  } // for
3232  SItemVec *vecdR = rightData->GetItemVec(axis);
3233  int cntR = 0; int cntR_LB = 0; int cntR_RB = 0; int cntR_LR = 0;
3234  for (SItemVec::iterator itr = vecdR->begin(); itr != vecdR->end(); itr++) {
3235    if ((*itr).obj->InBothLists())
3236      cntR_LR++;
3237    if ((*itr).IsLeftBoundary()) cntR_LB++;
3238    if ((*itr).IsRightBoundary()) cntR_RB++;
3239    cntR++;
3240  } // for
3241  assert(vecdL->size() == 2 * cntLout);
3242  assert(vecdR->size() == 2 * cntRout);
3243  assert(vecdL->size() == cntL_LB + cntL_RB);
3244  assert(vecdR->size() == cntR_LB + cntR_RB); 
3245  assert(cntL_LR == cntR_LR);
3246  assert(cntL_LB == cntL_RB);
3247  assert(cntR_LB == cntR_RB);
3248#endif
3249 
3250  return;
3251}
3252#endif
3253
3254// ------------------------------------------------------------------
3255// break the given list of SItem into two parts by reference axis
3256// and reference value, the output is in the first and second SItem
3257void
3258CKTBABuildUp::DivideAx_I(SInputData *d, int axis,
3259                         SInputData *rightData,
3260                         int &cntLout, int &cntRout)
3261{
3262  // First check, if this is not only singular case
3263  if (cntRout == 0) {
3264    assert(d->cntThickness == 0);
3265    return; // nothing to do
3266  }
3267 
3268  // the number of found boundaries on the left and on the right
3269  int cntLeft = 0;
3270  int cntRight = 0;
3271 
3272  // vector for left and right part
3273  SItemVec *vec = d->GetItemVec(axis);
3274#ifdef _DEBUG
3275  if (d->cntThickness == 0) {
3276    SItemVec::iterator it;
3277    for (it = vec->begin(); it != vec->end(); it++) {
3278      if ((*it).obj->InBothLists()) {
3279        cout << "ERROR in left 3 list" << endl;
3280      }
3281      if ((*it).obj->flags > 3) {
3282        cout << "Problem with 4 flags = " << (*it).obj->flags << endl;   
3283      }
3284    } // for   
3285  } 
3286#endif // _DEBUG 
3287 
3288  assert(vec->size() == d->count*2);
3289  SItemVec *vecRight = rightData->GetItemVec(axis);
3290  vecRight->resize(0);
3291 
3292  // Now go through source list
3293  SItemVec::iterator itDest = vec->begin();
3294  SItemVec::iterator itSrc = vec->begin();
3295  // traverse all boundaries, belonging only to the left part
3296  int i = 0;
3297  const int imax = d->count*2;
3298  for (; i != imax; itSrc++, i++) {   
3299    // Check if the item belongs to the second list
3300    if ((*itSrc).obj->InSecondList()) {
3301      break; // we can copy the
3302    }
3303    // Check if the item belongs to the first list
3304    if ((*itSrc).obj->InFirstList()) {
3305      // copy the content to the first list
3306      // Note that *itDest = *itSrc is not necessary, since itDest = itSrc
3307      cntLeft++;
3308    }
3309  } // for
3310
3311  // for newly copied object in the next step
3312  itDest = itSrc;
3313
3314  // the boundaries can belong to both parts
3315  for (; i != imax; itSrc++, i++) {   
3316    // Check if the item belongs to the second list
3317    if ((*itSrc).obj->InSecondList()) {
3318      vecRight->push_back((*itSrc));
3319      cntRight++;
3320    }
3321    // Check if the item belongs to the first list
3322    if ((*itSrc).obj->InFirstList()) {
3323      // copy the content to the first list
3324      *itDest = *itSrc;
3325      itDest++; // and move the iterator
3326      cntLeft++;
3327    }
3328  } // for
3329
3330  assert(cntLeft % 2 == 0);
3331  assert(cntRight % 2 == 0);
3332  assert(cntLeft == 2 * cntLout);
3333  assert(cntRight == 2 * cntRout);
3334
3335  vec->resize(cntLeft); 
3336
3337  return;
3338}
3339
3340// break the given list of SItem into two parts by reference axis
3341// and reference value, the output is in the first and second SItem
3342void
3343CKTBABuildUp::DivideAx_II(SInputData *d, int axis,
3344                          SInputData *rightData,
3345                          int &cntLout, int &cntRout)
3346{
3347  // First check, if this is not only singular case
3348  if (cntRout == 0) {
3349    assert(d->cntThickness == 0);
3350    return; // nothing to do
3351  }
3352 
3353  // the number of found boundaries on the left and on the right
3354  int cntLeft = 0;
3355  int cntRight = 0;
3356 
3357  // vector for left and right part
3358  SItemVec *vec = d->GetItemVec(axis);
3359#ifdef _DEBUG
3360  if (d->cntThickness == 0) {
3361    SItemVec::iterator it;
3362    for (it = vec->begin(); it != vec->end(); it++) {
3363      if ((*it).obj->InBothLists()) {
3364        cout << "ERROR in left list" << endl;
3365      }
3366      if ((*it).obj->flags > 3) {
3367        cout << "Problem with flags = " << (*it).obj->flags << endl;     
3368      }
3369    } // for   
3370  }
3371#endif // _DEBUG 
3372 
3373  assert(vec->size() == d->count*2);
3374  SItemVec *vecRight = rightData->GetItemVec(axis);
3375  vecRight->resize(0);
3376 
3377  // Now go through the source list
3378  SItemVec::iterator itDest = vec->begin();
3379  SItemVec::iterator itSrc = vec->begin();
3380  // traverse all boundaries, belonging only to the left part
3381  int i = 0;
3382  const int imax = d->count*2;
3383  for (; i != imax; itSrc++, i++) {   
3384    // Check if the item belongs to the second list
3385    if ((*itSrc).obj->InSecondList()) {
3386      break; // we can copy the
3387    }
3388    // Check if the item belongs to the first list
3389    if ((*itSrc).obj->InFirstList()) {
3390      // copy the content to the first list
3391      // Note that *itDest = *itSrc is not necessary, since itDest = itSrc
3392      cntLeft++;
3393    }
3394    // Reset flags for the next subdivision step
3395    if ((*itSrc).IsRightBoundary())
3396      (*itSrc).obj->ResetFlags();
3397  } // for
3398
3399  // for newly copied object in the next step
3400  itDest = itSrc;
3401
3402  // !!!!!!!!!!!!!!
3403  // for newly copied object in the next step
3404  itDest = itSrc;
3405
3406  // the boundaries can belong to both parts
3407  for (; i != imax; itSrc++, i++) {   
3408    //if ((*itSrc).obj == objtf) {
3409    //  cout << "DivideAxII - obj found2 , i = " << i << endl;
3410
3411    // Reset flags for the next subdivision step
3412    if ((*itSrc).IsRightBoundary()) {
3413      // This is the right boundary, we have to reset flags
3414      bool inFirstList = (*itSrc).obj->InFirstList();
3415      bool inSecondList = (*itSrc).obj->InSecondList();
3416      (*itSrc).obj->ResetFlags();
3417      if (inSecondList) {
3418        cntRight++;
3419        vecRight->push_back((*itSrc));
3420      }
3421       
3422      // Check if the item belongs to the first list
3423      if (inFirstList) {
3424        cntLeft++;
3425        // copy the content to the first list
3426        *itDest = *itSrc;
3427        itDest++; // and move the iterator
3428      }
3429    }
3430    else { // left boundary
3431      // Check if the item belongs to the second list
3432      if ((*itSrc).obj->InSecondList()) {
3433        cntRight++;
3434        vecRight->push_back((*itSrc));
3435      }
3436      // Check if the item belongs to the first list
3437      if ((*itSrc).obj->InFirstList()) {
3438        cntLeft++;
3439        // copy the content to the first list
3440        *itDest = *itSrc;
3441        itDest++; // and move the iterator
3442      }
3443    } // end of left boundary
3444  } // for
3445 
3446  assert(cntLeft % 2 == 0);
3447  assert(cntRight % 2 == 0);
3448  assert(cntLeft == 2 * cntLout);
3449  assert(cntRight == 2 * cntRout);
3450  vec->resize(cntLeft); 
3451
3452#ifdef _DEBUG
3453  for (SItemVec::iterator itt = vec->begin(); itt != vec->end(); itt++) {
3454    if ((*itt).obj->Flags() != 0) {
3455      cout << "DivideAxII ERROR 1 in final left list, flags = "
3456           << (*itt).obj->Flags() << " obj = " << (void*)((*itt).obj) << endl;
3457      abort();;
3458    }
3459  } // for
3460  for (SItemVec::iterator itt = vecRight->begin(); itt != vecRight->end(); itt++) {
3461    if ((*itt).obj->Flags() != 0) {
3462      cout << "DivideAxII ERROR 1 in final right list, flags = "
3463           << (*itt).obj->Flags() << " obj = " << (void*)((*itt).obj) << endl;
3464      abort();;
3465    }
3466  } // for
3467#endif // _DEBUG 
3468 
3469  return;
3470}
3471
3472
3473// break the given list of SItem into two parts by reference axis
3474// and reference value, the output is in the first and second SItem
3475void
3476CKTBABuildUp::DivideAx_I_opt(SInputData *d, int axis,
3477                             SInputData *rightData,
3478                             int cntLout, int cntRout)
3479{
3480  // First check, if this is not only singular case
3481  if (cntRout == 0) {
3482    assert(d->cntThickness == 0);
3483    return; // nothing to do
3484  }
3485 
3486  // vector for left and right part
3487  SItemVec *vec = d->GetItemVec(axis);
3488
3489#ifdef _DEBUG
3490  if (vec->size() != d->count*2) {
3491    cout << "vec->size() = " << vec->size()
3492         << " d->count*2 = " << d->count*2 << endl;
3493  }
3494
3495  if (d->cntThickness == 0) {
3496    SItemVec::iterator it;
3497    for (it = vec->begin(); it != vec->end(); it++) {
3498      if ((*it).obj->InBothLists()) {
3499        cout << "ERROR in left list" << endl;
3500      }
3501      if ((*it).obj->flags > 3) {
3502        cout << "Problem with flags = " << (*it).obj->flags << endl;     
3503      }
3504    } // for   
3505  } 
3506#endif // _DEBUG 
3507 
3508  assert(vec->size() == d->count*2);
3509  SItemVec *vecRight = rightData->GetItemVec(axis);
3510  vecRight->resize(0);
3511 
3512  // Now go through source list
3513  SItemVec::iterator itDest = vec->begin();
3514  SItemVec::iterator itSrc = vec->begin();
3515#if 1
3516  if (vec->size() != d->count * 2) {
3517    cout << "Something wrong HERE vec->size="
3518         << vec->size() << " d->count*2 = "
3519         << d->count*2 << endl;
3520    abort();
3521  }
3522#endif 
3523
3524  // traverse all boundaries, belonging only to the left part
3525  int i = 0;
3526  const int imax = d->count*2;
3527  for (; i != imax; itSrc++, i++) {   
3528    // Check if the item belongs to the second list
3529    if ((*itSrc).obj->InSecondList()) {
3530      break; // we can copy the
3531    }
3532  } // for
3533
3534  // for newly copied object in the next step
3535  itDest = itSrc;
3536 
3537  // the boundaries can belong to both parts
3538  for (; i != imax; itSrc++, i++) {
3539    SItem &itt = *itSrc;
3540    // Check if the item belongs to the second list   
3541    if (itt.obj->InSecondList()) {
3542      // and first copy the right boundary to a new list
3543      vecRight->push_back((*itSrc));
3544    }
3545    // Check if the item belongs to the first list
3546    if (itt.obj->InFirstList()) {
3547      // copy the content to the first list
3548      *itDest = *itSrc;
3549      itDest++; // and move the iterator
3550    }
3551  } // for
3552 
3553  // Shorten the left array by resizing
3554  vec->resize(cntLout*2);
3555  assert(vecRight->size() == cntRout*2);
3556 
3557  return;
3558}
3559
3560// break the given list of SItem into two parts by reference axis
3561// and reference value, the output is in the first and second SItem
3562void
3563CKTBABuildUp::DivideAx_II_opt(SInputData *d, int axis,
3564                              SInputData *rightData,
3565                              int cntLout, int cntRout)
3566{
3567  // First check, if this is not only singular case
3568  if (cntRout == 0) {
3569    assert(d->cntThickness == 0);
3570    return; // nothing to do
3571  }
3572   
3573  // vector for left and right part
3574  SItemVec *vec = d->GetItemVec(axis);
3575  assert(vec->size() == d->count*2);
3576#ifdef _DEBUG
3577  if (d->cntThickness == 0) {
3578    SItemVec::iterator it;
3579    for (it = vec->begin(); it != vec->end(); it++) {     
3580      if ((*it).obj->InBothLists()) {
3581        cout << "ERROR in left list" << endl;
3582      }
3583      if ((*it).obj->flags > 3) {
3584        cout << "Problem with flags = " << (*it).obj->flags << endl;     
3585      }
3586    } // for   
3587  }
3588  {
3589    SItemVec::iterator ittt;
3590    for (ittt = vec->begin(); ittt != vec->end(); ittt++) {
3591      if ((*ittt).obj->ToBeRemoved()) {
3592        cout << "ERROR - to be removed in the left list" << endl;
3593      }
3594    }
3595  }
3596#endif // _DEBUG 
3597 
3598  assert(vec->size() == d->count*2);
3599  SItemVec *vecRight = rightData->GetItemVec(axis);
3600  vecRight->resize(0);
3601 
3602  // Now go through the source list
3603  SItemVec::iterator itDest = vec->begin();
3604  SItemVec::iterator itSrc = vec->begin();
3605#if 1
3606  if (vec->size() != d->count * 2) {
3607    cout << "Something wrong HERE II vec->size="
3608         << vec->size() << " d->count*2 = "
3609         << d->count*2 << endl;
3610    abort();
3611  }
3612#endif 
3613  // traverse all boundaries, belonging only to the left part
3614  int i = 0;
3615  const int imax = d->count*2;
3616  for (; i != imax; itSrc++, i++) {
3617    //if ((*itSrc).obj == objtf) {
3618    //  cout << "DivideAxII - obj found1 , i = " << i << endl;
3619
3620    // Check if the item belongs to the second list
3621    if ((*itSrc).obj->InSecondList()) {
3622      break; // we can copy the
3623    }
3624    // Reset flags for the next subdivision step
3625    if ((*itSrc).IsRightBoundary())
3626      (*itSrc).obj->ResetFlags();
3627  } // for
3628
3629  // for newly copied object in the next step
3630  itDest = itSrc;
3631
3632  // the boundaries can belong to both parts
3633  for (; i != imax; itSrc++, i++) {   
3634    //if ((*itSrc).obj == objtf) {
3635    //  cout << "DivideAxII - obj found2 , i = " << i << endl;
3636
3637    bool inFirstList = (*itSrc).obj->InFirstList();
3638    bool inSecondList = (*itSrc).obj->InSecondList();
3639
3640    // Reset flags for the next subdivision step
3641    if ((*itSrc).IsRightBoundary()) {
3642      // This is the right boundary, we have to reset flags
3643      if (inSecondList)
3644        vecRight->push_back((*itSrc));
3645      // Check if the item belongs to the first list
3646      if (inFirstList) {
3647        // copy the content to the first list
3648        *itDest = *itSrc;
3649        itDest++; // and move the iterator
3650      }
3651    }
3652    else { // left boundary
3653      // Check if the item belongs to the second list
3654      if (inSecondList) {
3655        vecRight->push_back((*itSrc));
3656      }
3657      // Check if the item belongs to the first list
3658      if (inFirstList) {
3659        // copy the content to the first list
3660        *itDest = *itSrc;
3661        itDest++; // and move the iterator
3662      }
3663    } // end of left boundary
3664  } // for
3665
3666  // Shorten the left array by resizing
3667  vec->resize(cntLout*2);
3668  assert(vecRight->size() == cntRout*2);
3669
3670#if 0
3671#ifdef _DEBUG
3672  for (SItemVec::iterator itt2 = vec->begin(); itt2 != vec->end(); itt2++) {
3673    if ((*itt2).obj->Flags() != 0) {
3674      cout << "DivideAxII ERROR 1 in final left list, flags = "
3675           << (*itt2).obj->Flags() << " obj = "
3676           << (void*)((*itt2).obj) << endl;
3677      abort();;
3678    }
3679  } // for
3680  for (SItemVec::iterator itt3 = vecRight->begin(); itt3 != vecRight->end(); itt3++) {
3681    if ((*itt3).obj->Flags() != 0) {
3682      cout << "DivideAxII ERROR 1 in final right list, flags = "
3683           << (*itt3).obj->Flags() << " obj = "
3684           << (void*)((*itt3).obj) << endl;
3685      abort();;
3686    }
3687  } // for
3688#endif // _DEBUG 
3689#endif
3690 
3691  return;
3692}
3693
3694#ifdef _VYPIS
3695void
3696CKTBABuildUp::PrintOut(CKTBAxes::Axes axis, int stat, SItem **first,
3697                       SItem **second, float /*position*/)
3698{
3699#if 1
3700  int cnt = 0;
3701  SItem *t = *first;
3702  if (stat & 1) {
3703    DEBUG << "----FIRST--- axis=" << (int)axis << "\n";
3704    while (t != NULL) {
3705#define _VYP
3706
3707#ifdef _VYP
3708      float v = t->value[axis];
3709      DEBUGW << "curr=" << t
3710#ifdef _BIDIRLISTS
3711             << " prev=" << t->prev[axis]
3712#endif
3713             << " next="
3714             << t->next[axis] << " val=" << v
3715             << " " << (t->IsRightBoundary() ? 1 : 0);
3716      if (t->IsRightBoundary())
3717        DEBUGW << "   obj= " << t->obj;
3718      else
3719        DEBUGW << "   pair= " << t->pairF;
3720     
3721      DEBUGW << endl;
3722#endif // _VYP
3723      t = t->next[axis];
3724      cnt++;
3725    }
3726   
3727    DEBUG<< "FIRST cnt= " << cnt <<"\n" << flush;
3728  }
3729
3730  if (stat & 2) {
3731    // the right list
3732    t = *second;
3733    cnt = 0;
3734    DEBUG << "----SECOND-----\n";
3735    while (t != NULL) {
3736#ifdef _VYP
3737      float v = t->value[axis];
3738      DEBUGW << "curr=" << t
3739#ifdef _BIDIRLISTS
3740             << " prev=" << t->prev[axis]
3741#endif
3742             << " next="
3743            << t->next[axis] << " val=" << v
3744            << " " << ((t->IsRightBoundary()) ? 1 : 0);
3745      if (t->IsRightBoundary())
3746        DEBUGW << "   obj= " << t->obj;
3747      else
3748        DEBUGW << "   pair= " << t->pairF;
3749      DEBUGW << "\n" << flush;
3750#endif // _VYP
3751      t = t->next[axis];
3752      cnt++;
3753    }
3754    DEBUGW << "SECOND cnt= " << cnt <<"\n" << flush;
3755  }
3756#endif
3757}
3758#endif // _VYPIS
3759
3760// --------------------------------------------------------------------
3761// class CKTBABuildUp::CTestAx
3762
3763// The initialization for the first axis to be tested, this time *****Z axis*******
3764void
3765CKTBABuildUp::SSplitState::InitXaxis(int cnt, const SBBox &boxN)
3766{
3767  // initialize the variables, mainly for surface area estimates
3768  cntAll = cnt;     // the number of all objects in the bounding box
3769  cntLeft = 0;      // the count of bounding boxes on the left
3770  cntRight = cnt;   // the count of bounding boxes on the right
3771  thickness = 0;    // the count of bounding boxes straddling the splitting plane
3772  axis = CKTBAxes::EE_X_axis; // the axis, where the splitting is proposed
3773  box = boxN;        // the box, that is subdivided
3774  //int topAxis = 1;        // axis that is considered in top .. depth
3775  //int frontAxis = 2;      // axis that is considered in front .. height 
3776  // the size of bounding box along the axis to be split
3777  width = box.Max().x - box.Min().x;
3778  // and along two other axes 
3779  topw = box.Max().y - box.Min().y;
3780  frontw = box.Max().z - box.Min().z;
3781  // surface area of the splitting plane .. one face !!
3782  areaSplitPlane = topw * frontw;
3783  // the sum of length of the boxes not to be subdivided
3784  areaSumLength = topw + frontw;
3785  // The surface are of the whole box
3786  areaWholeSA2 = width * areaSumLength + areaSplitPlane;
3787  // initial evaluation .. the worst cost
3788  bestCost = WorstEvaluation();
3789}
3790
3791// The initialization for the first axis to be tested, this time *****Y axis*******
3792// This can be run independently.
3793void
3794CKTBABuildUp::SSplitState::InitYaxis(int cnt, const SBBox &boxN)
3795{
3796  // initialize the variables, mainly for surface area estimates
3797  cntAll = cnt;     // the number of all objects in the bounding box
3798  cntLeft = 0;      // the count of bounding boxes on the left
3799  cntRight = cnt;   // the count of bounding boxes on the right
3800  thickness = 0;    // the count of bounding boxes straddling the splitting plane
3801  axis = CKTBAxes::EE_Y_axis; // the axis, where the splitting is proposed
3802  box = boxN;        // the box, that is subdivided
3803  //int frontAxis = oaxes[axis][0]; // axis that is considered in front .. height - x-axis
3804  //int topAxis = oaxes[axis][1]; // axis that is considered in top .. depth - y-axis
3805  // the size along the second axis - height
3806  frontw = box.Max().x - box.Min().x;
3807  // the size of bounding box along the axis to be split - width
3808  width = box.Max().y - box.Min().y;
3809  // The size along third axis - depth
3810  topw = box.Max().z - box.Min().z;
3811  // surface area of the splitting plane .. one face !!
3812  areaSplitPlane = topw * frontw;
3813  // the sum of length of the boxes not to be subdivided
3814  areaSumLength = topw + frontw;
3815  areaWholeSA2 = width * areaSumLength + areaSplitPlane;
3816  // initial evaluation .. the worst cost
3817  bestCost = WorstEvaluation();
3818}
3819
3820// The initialization for the first axis to be tested, this time *****Z axis*******
3821// This can be run independently.
3822void
3823CKTBABuildUp::SSplitState::InitZaxis(int cnt, const SBBox &boxN)
3824{
3825  // initialize the variables, mainly for surface area estimates
3826  cntAll = cnt;     // the number of all objects in the bounding box
3827  cntLeft = 0;      // the count of bounding boxes on the left
3828  cntRight = cnt;   // the count of bounding boxes on the right
3829  thickness = 0;    // the count of bounding boxes straddling the splitting plane
3830  axis = CKTBAxes::EE_Z_axis; // the axis, where the splitting is proposed
3831  box = boxN;       // the box, that is subdivided
3832  //int frontAxis = oaxes[axis][1]; // axis that is considered in front .. height - x axis
3833  //int topAxis = oaxes[axis][0]; // axis that is considered in top .. depth - y axis
3834  // The size along third axis - depth
3835  topw = box.Max().x - box.Min().x;
3836  // the size along the second axis - height
3837  frontw = box.Max().y - box.Min().y;
3838  // the size of bounding box along the axis to be split - width
3839  width = box.Max().z - box.Min().z;
3840  // surface area of the splitting plane .. one face !!
3841  areaSplitPlane = topw * frontw;
3842  // the sum of length of the boxes not to be subdivided
3843  areaSumLength = topw + frontw;
3844  areaWholeSA2 = width * areaSumLength + areaSplitPlane;
3845  // initial evaluation .. the worst cost
3846  bestCost = WorstEvaluation();
3847}
3848
3849// --------------------------------------------------------------------
3850// The initialization for the first axis to be tested, use only in case
3851// when all 3 axes are tested. It can be called only when InitXaxis was called before !!!!
3852void
3853CKTBABuildUp::SSplitState::ReinitYaxis(int cnt, const SBBox &boxN)
3854{
3855  // initialize the variables, mainly for surface area estimates
3856  assert(cntAll == cnt); // the number of all objects in the bounding box
3857  cntLeft = 0;      // the count of bounding boxes on the left
3858  cntRight = cnt;   // the count of bounding boxes on the right
3859  thickness = 0;    // the count of bounding boxes straddling the splitting plane
3860  axis = CKTBAxes::EE_Y_axis; // the axis, where the splitting is proposed
3861  //int topAxis = 2;      // axis that is considered in top .. depth
3862  //int frontAxis = 0;    // axis that is considered in front .. height 
3863  // Swap the width, height, and depth
3864  float tS = width;
3865  // the size of bounding box along the axis to be split
3866  width = topw;
3867  // and along two other axes
3868  topw = frontw;
3869  frontw = tS; // copy the temporary variable
3870  // surface area of the splitting plane .. one face !!
3871  areaSplitPlane = topw * frontw;
3872  // the sum of length of the boxes not to be subdivided
3873  areaSumLength = topw + frontw;
3874  //assert(areaWholeSA2 == width * SumLength + areaSplitPlane);
3875}
3876
3877// The initialization for the first axis to be tested, use only in case
3878// when all 3 axes are tested. It can be called only when ReinitYaxis was called before !!!
3879void
3880CKTBABuildUp::SSplitState::ReinitZaxis(int cnt, const SBBox &boxN)
3881{
3882  // initialize the variables, mainly for surface area estimates
3883  assert(cntAll == cnt); // the number of all objects in the bounding box
3884  cntLeft = 0;      // the count of bounding boxes on the left
3885  cntRight = cnt;   // the count of bounding boxes on the right
3886  thickness = 0;    // the count of bounding boxes straddling the splitting plane
3887  axis = CKTBAxes::EE_Z_axis; // the axis, where the splitting is proposed
3888  // Swap the width, height, and depth
3889  float tS = width;
3890  // the size of bounding box along the axis to be split
3891  width = topw;
3892  // and along two other axes 
3893  topw = frontw;
3894  frontw = tS; // copy the temporary variable
3895  // surface area of the splitting plane .. one face !!
3896  areaSplitPlane = topw * frontw;
3897  // the sum of length of the boxes not to be subdivided
3898  areaSumLength = topw + frontw;
3899  //assert(areaWholeSA2 == width * stateSumLength + areaSplitPlane);
3900}
3901
3902// This is the computation of the cost using surface area heuristics
3903// using LINEAR EXTIMATE
3904void
3905CKTBABuildUp::EvaluateCost(SSplitState &state)
3906{
3907  // the surface area of the bounding box for the left child
3908  assert(state.position > -1e-9);
3909  float areaLeft = state.position * state.areaSumLength + state.areaSplitPlane;
3910
3911  // the surface area of the right bounding box for the right child
3912  assert(state.width - state.position > -1e-9);
3913  float areaRight = (state.width - state.position) * state.areaSumLength +
3914                    state.areaSplitPlane;
3915
3916  // computation of the cost .. smaller is better
3917  float cost = areaLeft * state.cntLeft + areaRight * state.cntRight;
3918
3919#ifdef _DEBUG_COSTFUNCTION
3920  // Put there normalized position and cost also normalized somehow
3921  ReportCostStream(state.position / state.width,
3922                   cost / (state.areaWholeSA2 * state.cntAll));
3923#endif
3924 
3925  // This normalization is not in fact necessary here
3926  //cost //= state.areaWholeSA2;
3927  if (cost < state.bestCost) {
3928    state.bestCost = cost;
3929    // The iterator pointing to the best boundary
3930    state.bestIterator = state.it;
3931    // The number of objects whose boundaries are duplicated
3932    state.bestThickness = state.thickness;
3933  }
3934  return;
3935}
3936
3937// This is the computation of the cost using surface area heuristics
3938void
3939CKTBABuildUp::EvaluateCostFreeCut(SSplitState &state)
3940{
3941  // This is assumed to work for free cut (no object intersected)
3942  assert(state.thickness == 0);
3943 
3944  // the surface area of the bounding box for the left child
3945  assert(state.position > -1e-9);
3946  float areaLeft = state.position * state.areaSumLength + state.areaSplitPlane;
3947
3948  // the surface area of the right bounding box for the right child
3949  assert(state.width - state.position > -1e-9);
3950  float areaRight = (state.width - state.position) * state.areaSumLength +
3951                    state.areaSplitPlane;
3952
3953  // computation of the cost .. smaller is better
3954  float cost = biasFreeCuts *(areaLeft * state.cntLeft + areaRight * state.cntRight);
3955
3956#ifdef _DEBUG_COSTFUNCTION
3957  // Put there normalized position and cost also normalized somehow
3958  ReportCostStream(state.position / state.width,
3959                   cost / (state.areaWholeSA2 * state.cntAll));
3960#endif
3961 
3962  // This normalization is not in fact necessary here
3963  //cost //= state.areaWholeSA2;
3964  if (cost < state.bestCost) {
3965    state.bestCost = cost;
3966    // The iterator pointing to the best boundary
3967    state.bestIterator = state.it;
3968    // The number of objects whose boundaries are duplicated
3969    state.bestThickness = 0;
3970  }
3971  return;
3972}
3973
3974// -------------------------------------------------------------------------------
3975// TRIAL TO IMPROVE for a given X-axis search for the best splitting plane position.
3976// Starts from start item, the rest of the information is set to 'state' variable.
3977void
3978CKTBABuildUp::GetSplitPlaneOpt(SItemVec *vec, int axisToTest)
3979{
3980#ifdef _DEBUG
3981  if ((vec == NULL) || (vec->size() == 0)) {
3982    cerr << "Trying to subdivide some node without an object" << endl;
3983    abort();;
3984  }
3985#endif // _DEBUG 
3986
3987#if 0
3988  static int index = 0;
3989  cout << "index = " << index << endl;
3990  const int indexToFind = 8;
3991  if (index == indexToFind) {
3992    cout << "Index found = " << index << endl;
3993  } 
3994  index++;
3995#endif
3996 
3997  // some necessary space is required to cut off
3998  const float eps = 1.0e-6 * state.width;
3999  const float MinPosition = state.box.Min(axisToTest);
4000  const float MaxPosition = state.box.Max(axisToTest);
4001  //float MinPositionAccept = MinPosition + eps;
4002  //float MaxPositionAccept = MaxPosition - eps;
4003
4004  // wherefrom to start
4005  SItemVec::iterator curr = vec->begin();
4006  SItemVec::iterator next = vec->begin() + 1;
4007  // Set the type of splitting by this function, which is in this case
4008  // not overwritten in Evaluate() function
4009  state.bestTwoSplits = 1;
4010  float val = (*curr).pos;
4011  float nval;
4012  float pval = val;
4013
4014  // Evaluate the first possible position, cutting off empty
4015  // space on the left of the splitting plane
4016  state.position = val - MinPosition;
4017  if (state.position > 0.f) {
4018    state.position2 = next->pos - MinPosition;
4019    state.it = curr;
4020    // Here evaluate the cost of surface area heuristics
4021    EvaluateCostFreeCut(state);
4022  }
4023 
4024  //-------------------- TEST LOOP ---------------------------------------
4025  // make evaluation for each splitting position
4026  const int imax = vec->size()-1;
4027  int ii;
4028  for (ii = 0; ii < imax; ii++)
4029  {
4030    nval = next->pos;
4031
4032#if 0
4033    if (axisToTest == 2) {
4034      if ((val > 0.866) && (val < 0.870)) {       
4035        cout << "val = " << val;
4036        cout << "  nval = " << nval;
4037        if (curr->IsRightBoundary())
4038          cout << " R ";
4039        else
4040          cout << " L ";
4041        cout << " thickness = " << state.thickness << endl;
4042       
4043      }
4044    }
4045#endif   
4046    // update left box and rightbox properties
4047    if (curr->IsRightBoundary()) { // we are on the right boundary of an object
4048      state.cntRight--;
4049      state.thickness--;
4050      // Check possibly the position
4051      if ( (val != nval) ||
4052           ((curr->IsRightBoundary()) && (next->IsLeftBoundary())) )
4053        {
4054          state.position = val - MinPosition;
4055          state.position2 = nval - MinPosition;
4056          state.it = curr;
4057          // Here evaluate the cost of surface area heuristics
4058          EvaluateCost(state);
4059        }     
4060    }
4061    else {
4062      // the left boundary .. enter the object from left
4063
4064      // Check possibly the position
4065      if ( (pval != val)
4066           // || ((curr->IsRightBoundary()) && (next->IsLeftBoundary()))
4067           )
4068        {
4069          state.position = val - MinPosition;
4070          state.position2 = nval - MinPosition;
4071          state.it = curr;
4072          // Here evaluate the cost of surface area heuristics
4073          EvaluateCost(state);
4074        }   
4075     
4076      state.cntLeft++;
4077      state.thickness++;
4078    }
4079       
4080    curr = next; // pointer to the boundary
4081    next++;
4082    pval = val;
4083    val = nval; // value of splitting plane
4084  } // while
4085
4086  //cout << "i = " << i << endl;
4087 
4088  assert(state.cntLeft == state.cntAll);
4089  state.cntRight--;
4090  assert(state.cntRight == 0);
4091  state.thickness--;
4092  assert(state.thickness == 0); 
4093 
4094  // Evaluate last possible position
4095  state.position = val - MinPosition;
4096  if (state.position < MaxPosition) {
4097    state.position2 = state.width;
4098    state.it = curr;
4099    // Here evaluate the cost of surface area heuristics
4100    EvaluateCostFreeCut(state);
4101  }
4102   
4103  // In this case the best result is kept in 'state' variable
4104  return;
4105} // CTestAx::GetSplitPlaneOpt (for X, Y, Z)
4106
4107// -------------------------------------------------------------------------------
4108// TRIAL TO IMPROVE for a given X/Y/Z-axis search for the best splitting plane position.
4109// Starts from start item, the rest of the information is set to 'state' variable.
4110void
4111CKTBABuildUp::GetSplitPlaneOpt2(SItemVec *vec, int axisToTest)
4112{
4113#ifdef _DEBUG
4114  if ((vec == NULL) || (vec->size() == 0)) {
4115    cerr << "Trying to subdivide some node without an object" << endl;
4116    abort();;
4117  }
4118#endif // _DEBUG 
4119 
4120  // some necessary space is required to cut off
4121  const float eps = 1.0e-6 * state.width;
4122  const float MinPosition = state.box.Min(axisToTest);
4123  const float MaxPosition = state.box.Max(axisToTest);
4124  //float MinPositionAccept = MinPosition + eps;
4125  //float MaxPositionAccept = MaxPosition - eps;
4126
4127  // wherefrom to start
4128  SItemVec::iterator curr = vec->begin();
4129  SItemVec::iterator next = vec->begin() + 1;
4130  // Set the type of splitting by this function, which is in this case
4131  // not overwritten in Evaluate() function
4132  state.bestTwoSplits = 1;
4133  state.thickness = 0;
4134  float val = (*curr).pos;
4135  float nval;
4136  float pval = val;
4137
4138//#if 1
4139#ifdef _DEBUG_COSTFUNCTION
4140  static int indexToFind = 0;
4141  static int indexSS = 0;
4142  if ( (indexSS == indexToFindSS) &&
4143       (!_alreadyDebugged)) {
4144    cout << "Debugging cost function, index = " << indexSS << endl;
4145    InitCostStream(indexSS, state.cntAll,
4146                   (val-state.box.Min(axisToTest))/state.width,
4147                   //(val-MinPosition)/state.width,
4148                   //(MinPosition)/state.width,
4149                   (state.box.Max(axisToTest)-MinPosition)/state.width);
4150  }
4151  indexSS++;
4152#endif
4153 
4154  // Evaluate the first possible position, cutting off empty
4155  // space on the left of the splitting plane
4156  state.position = val - MinPosition;
4157  if (state.position > 0.f) {
4158    state.position2 = next->pos - MinPosition;
4159    state.it = curr;
4160    // Here evaluate the cost of surface area heuristics
4161    EvaluateCostFreeCut(state);
4162  }
4163 
4164  //-------------------- TEST LOOP ---------------------------------------
4165  // make evaluation for each splitting position
4166  const int imax = vec->size()-1;
4167  int ii;
4168  for (ii = 1; ii <= imax; ii++)
4169  {
4170    nval = (*vec)[ii].pos;
4171    // update left box and rightbox properties
4172    if (curr->IsRightBoundary()) {
4173      // the right boundary of an object
4174      state.cntRight--;
4175      state.thickness--;
4176     
4177      // Check possibly the position
4178      if ( (val != nval) ||
4179           ((curr->IsRightBoundary()) && (next->IsLeftBoundary())) )
4180        {
4181          if (state.thickness >= 0) {
4182            state.position = val - MinPosition;
4183            state.position2 = nval - MinPosition;
4184            state.it = curr;
4185            // Here evaluate the cost of surface area heuristics
4186            EvaluateCost(state);
4187          }
4188        }
4189    }
4190    else {
4191      // the left boundary .. enter the object from left
4192
4193      // Check possibly the position
4194      if (pval != val)
4195      {
4196        if (state.thickness >= 0) {
4197          state.position = val - MinPosition;
4198          state.position2 = nval - MinPosition;
4199          state.it = curr;
4200          // Here evaluate the cost of surface area heuristics
4201          EvaluateCost(state);
4202        }
4203      }     
4204      state.cntLeft++;
4205      state.thickness++;
4206    }
4207       
4208    curr = next; // pointer to the boundary
4209    pval = val;
4210    val = nval; // value of splitting plane
4211    next++;
4212  } // while
4213
4214  //cout << "i = " << i << endl;
4215 
4216  assert(state.cntLeft == state.cntAll);
4217  state.cntRight--;
4218  assert(state.cntRight == 0);
4219  state.thickness--;
4220  assert(state.thickness == 0); 
4221
4222#if 0
4223  if (state.thickness < 0) {
4224    cout << "SSSomething wrong happened here at end\n" << endl;
4225  }
4226#endif     
4227
4228#ifdef _DEBUG
4229  if (curr != vec->end()-1) {
4230    cerr << "Implementation bug in pointers" << endl;
4231    cout << "curr = " << curr - vec->begin()
4232         << " vec->end()-1 = " << vec->end()-1 - vec->begin() << endl;
4233    abort();;
4234  }
4235#endif
4236 
4237  // Evaluate last possible position
4238  state.position = val - MinPosition;
4239  if (state.position < MaxPosition) {
4240    state.position2 = state.width;
4241    state.it = curr;
4242    // Here evaluate the cost of surface area heuristics
4243    EvaluateCostFreeCut(state);
4244  }
4245 
4246#ifdef _DEBUG_COSTFUNCTION
4247  CloseCostStream();
4248#endif
4249 
4250  // In this case the best result is kept in 'state' variable
4251  return;
4252} // CTestAx::GetSplitPlaneOpt (for X, Y, Z)
4253
4254// -------------------------------------------------------------------------------
4255// TRIAL TO IMPROVE for a given X/Y/Z-axis search for the best splitting plane position.
4256// Using unrolling possibly understandable for inteligent compiler.
4257// Starts from start item, the rest of the information is set to 'state' variable.
4258void
4259CKTBABuildUp::GetSplitPlaneOpt3(SItemVec *vec, int axisToTest)
4260{
4261#ifdef _DEBUG
4262  if ((vec == NULL) || (vec->size() == 0)) {
4263    cerr << "Trying to subdivide some node without an object" << endl;
4264    abort();;
4265  }
4266#endif // _DEBUG 
4267
4268#if 0
4269  static int index = 0;
4270  cout << "index = " << index << endl;
4271  const int indexToFind = 27;
4272  if (index == indexToFind) {
4273    cout << "Index found = " << index << endl;
4274  } 
4275  index++;
4276#endif
4277 
4278  // some necessary space is required to cut off
4279  const float eps = 1.0e-6 * state.width;
4280  const float MinPosition = state.box.Min(axisToTest);
4281  const float MaxPosition = state.box.Max(axisToTest);
4282  //float MinPositionAccept = MinPosition + eps;
4283  //float MaxPositionAccept = MaxPosition - eps;
4284
4285  // wherefrom to start
4286  SItemVec::iterator curr = vec->begin();
4287  SItemVec::iterator next = vec->begin() + 1;
4288  // Set the type of splitting by this function, which is in this case
4289  // not overwritten in Evaluate() function
4290  state.bestTwoSplits = 1;
4291  float val = (*curr).pos;
4292  // value at previous and next position
4293  float pval, nval;
4294  pval = MinPosition;
4295
4296  // Evaluate the first possible position, cutting off empty
4297  // space on the left of the splitting plane
4298  state.position = val - MinPosition;
4299  // The first boundary must be left
4300  assert(curr->IsLeftBoundary());
4301  if (state.position > 0.f) {   
4302    state.position2 = next->pos - MinPosition;
4303    state.it = curr;
4304    // Here evaluate the cost of surface area heuristics
4305    EvaluateCostFreeCut(state);
4306  } // ------------------------------
4307  // Update for the next iteration
4308  state.cntLeft++;
4309  state.thickness++;
4310  pval = val;
4311 
4312  //-------------------- TEST LOOP ---------------------------------------
4313  // make evaluation for each splitting position, minus the first and
4314  // the last one
4315  int imax = vec->size()-2;
4316  // Setting unrolling
4317  const int MaxUNROLL = 4;
4318  const int loops = imax / MaxUNROLL;
4319  const int remainder = imax % MaxUNROLL;
4320  // There is some bug, when using this version, the boundaries are
4321  // incorrectly sorted.
4322  assert(remainder < MaxUNROLL);
4323  int ii = 1;
4324  if (loops > 0) {
4325    // The outer loop
4326    for (int l = 0; l < loops; l++, ii += MaxUNROLL)
4327    {
4328      // This has to be unrolled by compiler such as Intel Compiler or GCC4.1
4329      // The number of loops is known in advance !
4330      for (int qq = 0; qq < MaxUNROLL; qq++)
4331      {
4332        curr = vec->begin() + ii + qq;
4333        next = curr + 1;
4334        #if 1
4335        val = (*curr).pos;
4336        nval = (*next).pos;
4337        #else
4338        val = (*vec)[ii+qq].pos;
4339        nval = (*vec)[ii+qq+1].pos;
4340        #endif
4341       
4342        // update left box and rightbox properties
4343        if (curr->IsRightBoundary()){// we are on the right boundary of an object
4344          state.cntRight--;
4345          state.thickness--;
4346          // Check possibly the position
4347          if ( (val != nval) ||
4348               ((curr->IsRightBoundary()) && (next->IsLeftBoundary())) )
4349          {
4350            state.position = val - MinPosition;
4351            state.position2 = nval - MinPosition;
4352            state.it = curr;
4353            // Here evaluate the cost of surface area heuristics
4354            EvaluateCost(state);
4355            // It does make sense to change boundary value
4356            pval = (*vec)[ii+qq].pos;
4357          }
4358        }
4359        else {
4360          // the left boundary .. enter the object from left
4361       
4362          // Check possibly the position
4363          if (pval != val)
4364          {
4365            state.position = val - MinPosition;
4366            state.position2 = nval - MinPosition;
4367            state.it = curr;
4368            // Here evaluate the cost of surface area heuristics
4369            EvaluateCost(state);
4370            // It makes sense to change boundary value
4371            pval = (*vec)[ii+qq].pos;
4372          }
4373          state.cntLeft++;
4374          state.thickness++;
4375        } // left or right boundary     
4376      } // for qq
4377    } // for ii
4378
4379    // Set the correct pointer for the next iteration
4380    curr = next; next++;
4381    pval = val;
4382  }
4383  else {
4384    //cout << "No loop" << endl;
4385    pval = val;
4386    curr = next;
4387    next++;
4388  }
4389
4390  // The rest of the loop
4391  if (ii <= imax) {
4392    // Set the correct pointer 
4393    for ( ;ii <= imax; ii++) {
4394      nval = curr->pos;
4395      // update left box and rightbox properties
4396      if (curr->IsRightBoundary()) { // we are on the right boundary of an object
4397        state.cntRight--;
4398        state.thickness--;
4399        // Check possibly the position
4400        if ( (val != nval) ||
4401             ((curr->IsRightBoundary()) && (next->IsLeftBoundary())) )
4402        {
4403          state.position = val - MinPosition;
4404          state.position2 = nval - MinPosition;
4405          state.it = curr;
4406          // Here evaluate the cost of surface area heuristics
4407          EvaluateCost(state);
4408        }
4409      }
4410      else {
4411        // the left boundary .. enter the object from left
4412     
4413        // Check possibly the position
4414        if (pval != val)
4415        {
4416          state.position = val - MinPosition;
4417          state.position2 = nval - MinPosition;
4418          state.it = curr;
4419          // Here evaluate the cost of surface area heuristics
4420          EvaluateCost(state);
4421        }
4422        state.cntLeft++;
4423        state.thickness++;
4424      } // if right/left boundary
4425   
4426      curr = next; // pointer to the boundary
4427      pval = val;
4428      val = nval; // value of splitting plane
4429      next++;
4430    } // for qq
4431
4432    // For the last evaluation
4433  }
4434  else {
4435    if (!loops) {
4436      curr = next;
4437      val = nval;
4438    }
4439  }
4440
4441  //cout << "i = " << i << endl;
4442
4443#ifdef _DEBUG 
4444  if (curr != vec->end()-1) {
4445    cerr << "Implementation bug in pointers" << endl;
4446    cout << "curr = " << curr - vec->begin()
4447         << "vec->end()-1 = " << vec->end()-1 - vec->begin() << endl;
4448    abort();;
4449  }
4450#endif
4451 
4452  assert(state.cntLeft == state.cntAll);
4453  state.cntRight--;
4454  assert(state.cntRight == 0);
4455  state.thickness--;
4456  assert(state.thickness == 0);
4457 
4458  // Evaluate last possible position
4459  state.position = val - MinPosition;
4460  if (state.position < MaxPosition) {
4461    state.position2 = state.width;
4462    state.it = curr;
4463    // Here evaluate the cost of surface area heuristics
4464    EvaluateCostFreeCut(state);
4465  }
4466   
4467  // In this case the best result is kept in 'state' variable
4468  return;
4469} // CTestAx::GetSplitPlaneOpt3 (for X, Y, Z)
4470
4471#if 0
4472// Here is some bug, that is clearly visible for tree.nff !! The resulting
4473// tree is much less efficient. 3/1/2006 VH.
4474// -------------------------------------------------------------------------------
4475// TRIAL TO IMPROVE for a given X,Y,Z-axis search for the best splitting plane position.
4476// Starts from start item, the rest of the information is set to 'state' variable.
4477// It does not work for gcc compiler better than GetSplitPlaneOpt
4478void
4479CKTBABuildUp::GetSplitPlaneOptUnroll4(SItemVec *vec, int axisToTest)
4480{
4481#ifdef _DEBUG
4482  if ((vec == NULL) || (vec->size() == 0)) {
4483    cerr << "Trying to subdivide some node without an object" << endl;
4484    abort();;
4485  }
4486#endif // _DEBUG 
4487
4488#if 0
4489  static int index = 0;
4490  cout << "index = " << index << endl;
4491  const int indexToFind = 8;
4492  if (index == indexToFind) {
4493    cout << "Index found = " << index << endl;
4494  } 
4495  index++;
4496#endif
4497 
4498  // some necessary space is required to cut off
4499  const float eps = 1.0e-6 * state.width;
4500  const float MinPosition = state.box.Min(axisToTest);
4501  const float MaxPosition = state.box.Max(axisToTest);
4502  //float MinPositionAccept = MinPosition + eps;
4503  //float MaxPositionAccept = MaxPosition - eps;
4504
4505  // wherefrom to start
4506  // Set the type of splitting by this function, which is in this case
4507  // not overwritten in Evaluate() function
4508  state.bestTwoSplits = 1;
4509  float val4 = MinPosition;
4510 
4511  //-------------------- TEST LOOP ---------------------------------------
4512  // make evaluation for each splitting position
4513
4514  // the first position was already evaluated and the last is also a special case
4515  const int imax = vec->size() - 1;
4516  int cntDebug = 0;
4517  int imax4 = 4*(imax/4);
4518  int imax4rest = imax % 4;
4519  int i;
4520  // This is manual unrolling
4521  if (imax4 >= 4) {
4522  for(i = 0; i < imax4; i += 4)
4523  {
4524    // The evaluation I
4525    cntDebug++;
4526    float pval1 = val4;
4527    SItem* curr1 = &((*vec)[i+0]); float val1 = (*vec)[i+0].pos;
4528    SItem* next1 = &((*vec)[i+1]); float nval1 = (*vec)[i+1].pos;
4529    if (curr1->IsRightBoundary()) {
4530      state.cntRight--; state.thickness--;
4531      if ( (val1 != nval1) ||
4532           ((curr1->IsRightBoundary()) && (next1->IsLeftBoundary())) ) {
4533        state.position = val1 - MinPosition; state.position2 = nval1 - MinPosition;
4534        state.it = vec->begin() + i; EvaluateCost(state);
4535      }
4536    }
4537    else {
4538      if (pval1 != val1) {
4539          state.position = val1 - MinPosition; state.position2 = nval1 - MinPosition;
4540          state.it = vec->begin() + i; EvaluateCost(state);
4541      }
4542      state.cntLeft++; state.thickness++;
4543    }
4544    // The evaluation II
4545    cntDebug++;
4546    float pval2 = (*vec)[i+0].pos;
4547    SItem* curr2 = &((*vec)[i+1]); float val2 = (*vec)[i+1].pos;
4548    SItem* next2 = &((*vec)[i+2]); float nval2 = (*vec)[i+2].pos;
4549    if (curr2->IsRightBoundary()) {
4550      state.cntRight--; state.thickness--;
4551      if ( (val2 != nval2) ||
4552           ((curr2->IsRightBoundary()) && (next2->IsLeftBoundary())) ) {
4553        state.position = val2 - MinPosition; state.position2 = nval2 - MinPosition;
4554        state.it = vec->begin() + i + 1; EvaluateCost(state);
4555      }
4556    }
4557    else {
4558      if (pval2 != val2) {
4559          state.position = val2 - MinPosition; state.position2 = nval2 - MinPosition;
4560          state.it = vec->begin() + i + 1; EvaluateCost(state);
4561      }
4562      state.cntLeft++; state.thickness++;
4563    }
4564    // The evaluation III
4565    cntDebug++;
4566    float pval3 = (*vec)[i+1].pos;
4567    SItem* curr3 = &((*vec)[i+2]); float val3 = (*vec)[i+2].pos;
4568    SItem* next3 = &((*vec)[i+3]); float nval3 = (*vec)[i+3].pos;
4569    if (curr3->IsRightBoundary()) {
4570      state.cntRight--; state.thickness--;
4571      if ( (val3 != nval3) ||
4572           ((curr3->IsRightBoundary()) && (next3->IsLeftBoundary())) ) {
4573        state.position = val3 - MinPosition; state.position2 = nval3 - MinPosition;
4574        state.it = vec->begin() + i + 2; EvaluateCost(state);
4575      }
4576    }
4577    else {
4578      if (pval3 != val3) {
4579          state.position = val3 - MinPosition; state.position2 = nval3 - MinPosition;
4580          state.it = vec->begin() + i + 2; EvaluateCost(state);
4581      }
4582      state.cntLeft++; state.thickness++;
4583    }
4584    // The evaluation IV
4585    cntDebug++;
4586    float pval4 = (*vec)[i+2].pos;
4587    SItem* curr4 = &((*vec)[i+3]); val4 = (*vec)[i+3].pos;
4588    SItem* next4 = &((*vec)[i+4]); float nval4 = (*vec)[i+4].pos;
4589    if (curr4->IsRightBoundary()) {
4590      state.cntRight--; state.thickness--;
4591      if ( (val4 != nval4) ||
4592           ((curr4->IsRightBoundary()) && (next4->IsLeftBoundary())) ) {
4593        state.position = val4 - MinPosition; state.position2 = nval4 - MinPosition;
4594        state.it = vec->begin() + i + 3; EvaluateCost(state);
4595      }
4596    }
4597    else {
4598      if (pval4 != val4) {
4599          state.position = val4 - MinPosition; state.position2 = nval4 - MinPosition;
4600          state.it = vec->begin() + i + 3; EvaluateCost(state);
4601      }
4602      state.cntLeft++; state.thickness++;
4603    }
4604  } // for, unrolled loop by factor 4
4605  } // if more than 4 evaluations
4606 
4607  for(i = imax4; i < imax; i ++)
4608  {
4609    cntDebug++;
4610    float pvalI = val4;
4611    SItem* currI = &((*vec)[i+0]); float valI = (*vec)[i+0].pos;
4612    SItem* nextI = &((*vec)[i+1]); float nvalI = (*vec)[i+1].pos;
4613    if (currI->IsRightBoundary()) {
4614      state.cntRight--; state.thickness--;
4615      if ( (valI != nvalI) ||
4616           ((currI->IsRightBoundary()) && (nextI->IsLeftBoundary())) ) {
4617        state.position = valI - MinPosition; state.position2 = nvalI - MinPosition;
4618        state.it = vec->begin() + i; EvaluateCost(state);
4619      }
4620    }
4621    else {
4622      if (pvalI != valI) {
4623          state.position = valI - MinPosition; state.position2 = nvalI - MinPosition;
4624          state.it = vec->begin() + i; EvaluateCost(state);
4625      }
4626      state.cntLeft++; state.thickness++;
4627    }
4628    val4 = valI;
4629  } // for i
4630 
4631  //cout << "i = " << i << endl;
4632  assert(cntDebug == imax);
4633 
4634  assert(state.cntLeft == state.cntAll);
4635  state.cntRight--;
4636  assert(state.cntRight == 0);
4637  state.thickness--;
4638  assert(state.thickness == 0); 
4639 
4640  // Evaluate last possible position
4641  state.position = (*vec)[imax-1].pos - MinPosition;
4642  if (state.position < MaxPosition) {
4643    state.position2 = state.width;
4644    state.it = vec->end() - 1;
4645    // Here evaluate the cost of surface area heuristics
4646    EvaluateCost(state);
4647  }
4648   
4649  // In this case the best result is kept in 'state' variable
4650  return;
4651} // CTestAx::GetSplitPlaneX,Y,Z
4652#endif
4653
4654}
4655
Note: See TracBrowser for help on using the repository browser.