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

// ============================================================================
// $Id: $
//
// ktbai.h
//     classes for building up the different KD-trees
//
// Class: CKTBBuildUp, CKTBBuildUp_new
//
// REPLACEMENT_STRING
//
// Initial coding by Vlasta Havran, February 2007

#ifndef __KTBAI_H__
#define __KTBAI_H__

// GOLEM headers
#include "configh.h"
#include "ktbconf.h"
#include "ktb.h"
#include "ktb8b.h"
#include "Containers.h"
#include "IntersectableWrapper.h"

namespace GtpVisibilityPreprocessor {

// forward declarations
class SKTBNode;

#ifndef _KTB8Bytes
// Use 12 Bytes representation
#define CKTBAllocManPredecessor CKTBAllocMan
#undef  SKTBNodeT
#define SKTBNodeT CKTBNodeAbstract::SKTBNode
#else
// Use 8 Bytes representation per node
#define CKTBAllocManPredecessor CKTB8BAllocMan
#undef  SKTBNodeT
#define SKTBNodeT CKTB8BNodeAbstract::SKTBNode
#endif

#ifndef INFINITY
#define INFINITY 10e10
#endif

class AxisAlignedBox3Intersectable:
    public IntersectableWrapper<AxisAlignedBox3>
{
public:
  AxisAlignedBox3Intersectable(const AxisAlignedBox3 &item):
    IntersectableWrapper<AxisAlignedBox3>(item) { }

  AxisAlignedBox3 GetBox() const { return mItem;}

  int Type() const
  {
    // This is not ture, but for our purposes it is OK
    return Intersectable::TRIANGLE_INTERSECTABLE;
  }

};

// ---------------------------------------------------------------  
// The base class for KD-tree with irregular change of axes, where
// the splitting plane can be positioned.
class CKTBABuildUp:
  public CKTBAllocManPredecessor
{
public:
  // The definition of flags
  enum EBoundaryType {
    EE_LeftBoundary = 1,
    EE_InLeftList = 1,
    EE_RightBoundary = 2,
    EE_InRightList = 2,
    EE_BothBoundaries = 3,
    EE_ToBeRemoved = 4
  };
  
  // the item in the list for all objects
  struct SSolid
  {
    Intersectable *obj; // pointer to the object itself
    unsigned int flags; // the flags to be set, they are common for all boundaries

    // query functions
    inline bool InFirstList() const { return (flags & EE_InLeftList); }
    inline bool InSecondList() const { return (flags & EE_InRightList); }
    inline bool InBothLists() const { return (flags == EE_BothBoundaries); }
    inline bool ToBeRemoved() const { return (flags & EE_ToBeRemoved); }
    inline bool ToBeRemovedOnly() const { return (flags == EE_ToBeRemoved); }
    inline unsigned int Flags() const { return flags;}
    // Setting functions
    inline void SetInFirstList() { flags |= 1; }
    inline void SetInSecondList() { flags |= 2; }
    inline void SetToRemove() { flags |= 4; }
    inline void SetToRemoveOnly() { flags = 4; }
    
    inline void ResetFlags() { flags = 0;}
    SSolid() { ResetFlags(); }
  };

  // The container of the object entries
  typedef vector<SSolid> SSolidVec;
  
  // the array of all objects in the scene
  SSolidVec solidArray;

  // the item of the boundary list - either left or right boundary
  // of the axis-aligned bounding box of the object. This structure
  // is intentionally of small size, namely 12 or 16 Bytes.

  struct SItem
  {
    float pos;  // boundary values for all three axes
    struct SSolid *obj; // the pointer to the object with flags

    // The axis represented by the item (CKTBAxes::Axes)
    uint1 axis; // = (X=0, Y=1, Z=2)
    // the type of boundary (low, high)
    uint1 typeLoHi;
    // only allignment to 12 Bytes
    //uint2 dummy;
    // -------------------------------------------------
    // some basic functions
    SItem(float posN, SSolid *objN, int axisN, EBoundaryType LoHiN) {
      pos = posN; obj = objN; axis = axisN; typeLoHi = (uint1)LoHiN;
    }
    // Simply constructor, just initializing flags
    SItem() {obj = 0; axis = 255; typeLoHi = 0; }
    //SItem& operator=(const SItem &src) {
    //  this->pos = src.pos; this->obj = src.obj;
    //  this->axis = src.axis; this->typeLoHi = src.typeLoHi;
    //  return *this;
    // }

    // Simple query functions
    inline bool IsLeftBoundary() const {
      return (typeLoHi ==  EE_LeftBoundary);
    }
    inline bool IsRightBoundary() const {
      return (typeLoHi == EE_RightBoundary);
    }    
    // Simple set functions
    void SetLeftBoundary() { typeLoHi = EE_LeftBoundary; }
    void SetRightBoundary() { typeLoHi = EE_RightBoundary; }
    // For quicksort
    friend bool operator<(const SItem &a, const SItem &b) {
      if (a.pos < b.pos)
        return -1;
      else
        if (a.pos > b.pos)
          return 1;
        else
          // the coordinates are equal
          if ( (a.IsRightBoundary()) &&
               (b.IsLeftBoundary()) )
            return -1; // right_boundary < left_boundary
          else
            if ( (a.IsLeftBoundary()) &&
                 (b.IsRightBoundary()) )
              return 1; // left_boundary > right_boundary
      // coordinates are equal, the same value and type, order is correct
      return 0; 
    }
  };

  // Here is the extended element for RadixSort
  struct SItemRadix:
    public SItem
  {
    // the pointer needed to chain the data during sorting
    SItemRadix *next;
    // Basic operations
    SItemRadix(): SItem() { next = 0;}
    // This is necessary constructor
    SItemRadix(const SItem &it) {
      memcpy(this, &it, sizeof(SItem));
      next = 0;
    }
    // This is necessary copy operator
    SItemRadix& operator=(const SItem &it) {
      memcpy(this, &it, sizeof(SItem));
      next = 0;
      return *this;
    }
  };
  

  // ---------------------------------------------------------
  // The declaration of container with object boundaries
  typedef vector<SItem> SItemVec;
  typedef vector<SItemRadix> SItemVecRadix;

  // ---------------------------------------------------------
  // Sorting by QuickSort and RadixSort

  // QuickSort
  // compare function for SItem*
  static int Compare(const SItem *p, const SItem *q);
  // bounding box sorting by Quick Sort
  void SortOneAxis(SItemVec &itemvec, int cnt, int * const stackQuickSort);

  // -------------------------------------
  // For Radix Sort
  // Radix sort int 2^8=256 classes and three
  // passes .. 4x8 bits=32 bits
  bool _useRadixSort;
  enum { RXBITS30 = 11 }; // the number of bits used for one phase of Radix Sort
  enum { // the number of buckets for RadixSort
         RXBUFS30 = 1 << RXBITS30,
         RXBUFS30_2 = 1 << (RXBITS30-1)};
  // This is one bucket of radix soft
  struct SRadix {
    SItemRadix *beg, *end;
  };
  // for 3-passes radix sort over the vectored data
  void CopyToAuxArray(const SItemVec &bounds, SItemVecRadix &aux);
  void RadixPassHoffset11(SItemVecRadix &bounds, int bit, SRadix *rb,
                          float offset, SItemRadix **start);
  void RadixPass11(SItemRadix **start, int cnt, int bit, SRadix *rb);
  void RadixPassOffset10(SItemRadix **start, int cnt, int bit, SRadix *rb,
                         float offset);
  void CopyFromAuxArray(SItemRadix *aux, SItemVec &bounds);

  // forward declaration
  struct SInputData;
  // sorts all three axes, cnt is the number of elems
  void SortAxes(SInputData *data);
  // initialization of the bounding box for a given object
  void LoadBB(SBBox &bb, SSolid *obj);

  // test if the lists are correctly sorted
  void Check3List(SInputData *data);
  void Check1List(SItemVec *vec, int axis, int countExpected);
  void Check1List(SInputData *data, int axis, int countExpected);
    
  //----------------------------------------------------------------
  // Termination criteria and fixing the splitting plane orientation

  // structure for prefered and required params for evaluation functions
  // and the termination criteria
  struct SReqPrefParams
  {
    //if any position on required axis is preferred for next subdivision step
    float reqPosition; // then reqPosition>0
    // if any axis is prefered for next step
    bool  useReqAxis;
    // the prescribed axis for the next subdivision
    CKTBAxes::Axes reqAxis;

    // -------------- AUTOMATIC TERMINATION CRITERIA ----------------
    // the ratio of improvement for the cost by subdivision and not-subdividing
    // for the previous subdivision
    float ratioLast;
    // the ratio of improvement for the subdivision in the previous step
    float ratioLastButOne;
    // the number of subdivision from the root node, where the improvement
    // in the cost failed
    int   failedSubDivCount;
    void Init() {
      reqPosition = Limits::Infinity;
      useReqAxis = false;
      reqAxis = CKTBAxes::EE_Leaf;
  
      ratioLast = 1000.0;
      ratioLastButOne = 1000.0;
      failedSubDivCount = 0;      
    }
  };

  // initialize required and preferenced parameters before first subdivision
  void InitReqPref(SReqPrefParams *pars);

  // ------------------------------------------------------
  // A structure for a single step of subdivision
  struct SInputData {
    // the traversal bounding box of the scene (not necessarily tight)
    SBBox box;
    // the number of objects in the node (= number_of_boundaries/2)
    int count;
    // the number of reserved boundaries in the node (>=2*count)
    int cntReserved;

    // The list of x-boundaries, y-boundaries, z-boundaries
    SItemVec *xvec;
    SItemVec *yvec;
    SItemVec *zvec;

    // only for allignment, it can be used for different purpose
    int algorithmBreakAx;

    // ----------------------------
    // The mode of subdivision
    ESubdivMode modeSubDiv;

    // Some prescribed parameters to be used
    SReqPrefParams pars;

    // ----------------------------------
    // Axis to be used if prescribed
    CKTBAxes::Axes axis;
    // the position to be used for MakeOneCut
    float position;
    float position2;
    // the number of objects to be duplicated
    int cntThickness;
    // the iterator to be used for splitting
    SItemVec::iterator bestIterator;
    // if 1 or 2 splits
    int twoSplits;
    // the best cost
    float bestCost;

    // if to make subdivision on the left node
    int makeSubdivisionLeft;
    // if to make subdivision on the right node
    int makeSubdivisionRight;

    // -----------------------------------
    // When the min boxes was inserted as the first one
    int lastDepthForMinBoxes;
    // The surface area for the last minimum box inserted
    float lastMinBoxSA2;
    // The pointer to the last inserted minimum box
    SKTBNodeT* lastMinBoxNode;
  private:
    void Init() {
      box.Initialize();
      algorithmBreakAx = 0;
      count = 0; cntReserved = 0;
      xvec = yvec = zvec = 0;
      pars.Init();
      cntThickness = 0;
      makeSubdivisionLeft = makeSubdivisionRight = 1;
      lastDepthForMinBoxes = 0;
      lastMinBoxSA2 = INFINITY;
      lastMinBoxNode = 0;
    }
  public:

    // -----------------------------------
    // Implicit constructor
    SInputData() {
      Init();
    }
    ~SInputData() {
      Free();
    }
    
    // Allocate at least for one object
    void Alloc(int sizeN = 2) {
      if (!xvec)
        xvec = new GALIGN16 vector<SItem>;
      assert(xvec);
      if (!yvec)
        yvec = new GALIGN16 vector<SItem>;
      assert(yvec);
      if (!zvec)
        zvec = new GALIGN16 vector<SItem>;
      assert(zvec);
      cntReserved = sizeN;
      // cout << "SizeN = " << sizeN << endl;
      xvec->reserve(sizeN); xvec->resize(0); 
      yvec->reserve(sizeN); yvec->resize(0); 
      zvec->reserve(sizeN); zvec->resize(0);
      count = 0;
    } // Alloc
    void Free() {
      delete xvec; xvec = 0;
      delete yvec; yvec = 0;
      delete zvec; zvec = 0;
      count = cntReserved = 0;
    } // Free
    void Reserve(int sizeN) {
      assert(xvec);
      assert(yvec);
      assert(zvec);
      if (sizeN > cntReserved) {
        xvec->reserve(sizeN);
        yvec->reserve(sizeN);
        zvec->reserve(sizeN);
        cntReserved = sizeN;
      }
    } // Reserve
    void Resize(int sizeN) {
      assert(sizeN >= 0);
      assert(xvec);
      assert(yvec);
      assert(zvec);
      xvec->resize(sizeN);
      yvec->resize(sizeN);
      zvec->resize(sizeN);
      count = sizeN*2;
    } // Reserve

    // Return the item using the index
    SItemVec* GetItemVec(int i) {
      assert((i >= 0) && (i < 3));
      return (&xvec)[i];
    }
    void CopyBasicData(SInputData *d) {
      box = d->box;
      count = 0;
      algorithmBreakAx = d->algorithmBreakAx;
      modeSubDiv = d->modeSubDiv;
      pars = d->pars;
      axis = d->axis;
      position = d->position;
      // added 12/2007 VH
      position2 = d->position2;
      cntThickness = d->cntThickness;
      //
      makeSubdivisionLeft = d->makeSubdivisionLeft;
      makeSubdivisionRight = d->makeSubdivisionRight;      
      // Intentionally, do not copy vectors of items
      lastDepthForMinBoxes = d->lastDepthForMinBoxes;
      lastMinBoxSA2 = d->lastMinBoxSA2;
      lastMinBoxNode = d->lastMinBoxNode;
    }
  };

  // Stack of data to be used
  SInputData *stackID;
  // current index
  int         stackIndex;
  // the maximum depth of tree
  int         maxTreeDepth;
  // the depth of the stack
  int         stackDepth;
  // Return the new data to be used
  SInputData* AllocNewData() {
    int i = stackIndex;
    stackIndex++;
    return &(stackID[i]);
  }
  SInputData* AllocNewData(int cnt) {
    int i = stackIndex;
    stackID[i].Alloc(cnt);
    stackIndex++;
    return &(stackID[i]);
  }
  // Free the last data allocated
  void FreeLastData() { stackIndex--; }
    
  // ---------------------------------------------------------------------
  // upper-level function for building up CKTB tree

  // creates all the auxiliary structures for building up CKTB tree
  SInputData* Init(ObjectContainer *objlist, const AxisAlignedBox3 &box);

  void DeleteAuxiliaryData() {
    for (int i = 0; i < stackDepth; i++) {
      stackID[i].Free();
    }
  }
  
  // ---------------------------------------------------------------------
  // Working with boundaries of objects
  
  // make the full leaf from current node
  SKTBNodeT* MakeLeaf(SInputData *i);

  // breaks the list into two list for a given axis and value
  void BreakAx(SInputData *i, int axis,
               SInputData *right,
               int &cntL, int &cntR);
  // breaks the list into two list for a given axis and value
  void BreakAxPosition(SInputData *i, int axis,
                       SInputData *right,
                       int &cntL, int &cntR);

  // split the list in the other than splitting axis into two lists
  void DivideAx_I(SInputData *i, int axis,
                  SInputData *right,
                  int &cntL, int &cntR);
  // also split and set the boundaries to be only in the first list
  void DivideAx_II(SInputData *i, int axis,
                   SInputData *right,
                   int &cntL, int &cntR);
  // split the list in the other than splitting axis into two lists
  void DivideAx_I_opt(SInputData *i, int axis,
                      SInputData *right,
                      int cntL, int cntR);
  // also split and set the boundaries to be only in the first list
  void DivideAx_II_opt(SInputData *i, int axis,
                       SInputData *right,
                       int cntL, int cntR);
  // reduce bounding boxes of objects split by the splitting plane
  void ReduceBBoxes(SInputData *i, int axis,
                    SInputData *right,
                    const float &position);
  // Remove the objects from the containter
  void RemoveObjects(SItemVec *, int cntObjects);
  void RemoveObjectsReset(SItemVec *, int cntObjects);

  // Computes the tight bounding box and the number of changed planes
  // when the tight box is used
  int GetEBox(const SInputData &i, SBBox &tbox);

  // returns a box enclosing all the objects in the node
  void GetTightBox(const SInputData &i, SBBox &tbox);
  
  // creates one cut inside CKTB tree
  SKTBNodeT* MakeOneCut(SInputData *i);
  // recursive function for creation of CKTB tree
  SKTBNodeT* SubDiv(SInputData *i);

  // ------ Methods for building up CKTB tree ------------------
  // returns 1 to supress to call the following criteria
  struct SSplitState
  {
    // counts
    int   cntAll;       // the number of all objects in the bounding box
    int   cntLeft;      // the count of bounding boxes on the left
    int   cntRight;     // the count of bounding boxes on the right
    int   thickness;    // the count of bounding boxes straddling the splitting plane

    CKTBAxes::Axes  axis;  // the axis, where the splitting is proposed
    float sizeb[3]; // the size of the box for x, y, and z
    SBBox box; // the box, that is subdivided
    
    // derived values from basic ones
    float width;        // the size of bounding box along the axis
    float frontw;       // the size of the bounding box in another axis (depth)
    float topw;         // the size of the bounding box in next next axis (height)
    float areaSplitPlane; // the area of the splitting plane
    float areaSumLength; // the size of the bounding as sum of height and depth
    float areaWholeSA2; // the half of the surface area of the whole box for this node

    // The iterator valid for current position
    SItemVec::iterator it;    
    // The position for this splitting plane to be evaluated
    float position; // the distance from the left boundary of the box for this node
    // The position for the next position, makes sense only for free interval (thickness=0)
    float position2; // the distance from the left boundary of the box for this node
    
    // The evaluation best cost until now
    float bestCost;
    // The position to be used
    SItemVec::iterator bestIterator;
    // The number of objects stradling the spliting plane for best position
    int   bestThickness;
    // Which mechanism to be used for splitting, either 0,1, or 2 splitting planes
    int   bestTwoSplits;

    // setting the evaluation for split cases that must not be done
    float WorstEvaluation() const { return MAXFLOAT;}
    
    // The initialization for the first axis to be tested.
    void InitXaxis(int cnt, const SBBox &box);
    void InitYaxis(int cnt, const SBBox &box);
    void InitZaxis(int cnt, const SBBox &box);
    // This function can be called only if InitXaxis was called before
    void ReinitYaxis(int cnt, const SBBox &box);
    // This function can be called only if InitXaxis was called, and subsequently
    // the function ReinitYaxis was called.
    void ReinitZaxis(int cnt, const SBBox &box);
    // Normalize the best cost by surface area of the box
    void NormalizeCostBySA2() { bestCost /= areaWholeSA2;}
  };

  // splitting state for current search
  SSplitState state;

  // Evaluating the cost, given the state and the values of splitting
  void EvaluateCost(SSplitState &state);

  // Evaluating the cost, given the state and the values of splitting
  // for free cuts
  void EvaluateCostFreeCut(SSplitState &state);
  
  // ----- statistical data ---------
  int cntDuplicate;     // count of duplicated objects until now
  bool resetFlagsForBreakAx; 

  // ------ debugging data ----------
  // if to print out the tree during construction
  bool _printCuts;
  
protected:

  // ---------------------------------
  // The selection of the axis
  int _algorithmForAxisSelection;
  
  //  ---- termination criteria -----
  int algorithmAutoTermination; // the algorithm for automatic termination criteria
  int maxDepthAllowed; // maximal depth of CKTB tree
  int maxListLength; // maximal list length of CKTB tree
  int maxCountTrials; // maximum number of trials for automatic termination criteria
  // the cutting off empty space in leaves
  bool cutEmptySpace; // if to cut off empty space in leaves in postprocessing
  int absMaxAllowedDepth; // maximal depth from the root - mut not be surpassed
  // maximal depth allowed for cutting within the leaf .. cut off empty space
  int maxEmptyCutDepth; // must be <0,1,2,3,4,5,6> since six planes are enough
  // This is working variable, denoting the depth of the leaf to be created.
  int startEmptyCutDepth;

  // Biasing the empty cuts (no objects are split). The cost is multiplied
  // by the coefficient which is assumed to be 0.8-0.9
  float  biasFreeCuts;
  
  // ---------- Special improvements on the kd-tree construction --------
  // flag if to split bounding boxes during splitting
  bool    splitClip;
  // flag if to put minimum enclosing boxes sparsely during the construction
  bool    makeMinBoxes;
  // if we make tight boxes if we put min box !
  bool    makeTightMinBoxes;
  // parameters to drive the minboxes construction
  int     minObjectsToCreateMinBox, minDepthDistanceBetweenMinBoxes;
  float   minSA2ratioMinBoxes;
  // Make min box here
  bool   makeMinBoxHere;
  
  // two next axes are stored in oaxes for each axis
  static const CKTBAxes::Axes oaxes[3][2];

  // ------ data to create the tree --------------------
  int         initcnt;      // initial number of objects
  SBBox       wBbox; // the box of the world in float values
  Vector3   boxSize;   // the size of world bounding box in float
  float       wholeBoxArea; // the surface area of the scene bounding box
  
  // for some functions it is necessary to have determined the following costs
  float       Ct; // traversal cost - going in given direction != decision
  float       Ci; // intersection cost - average intersection cost with object

  // just if to be verbose
  bool verbose;
  
  // the main function of this class .. returns the best splitting plane
  // in X axis, requires InitXaxis() to be called before or InitYaxis() or InitZaxis()
  // optimized version
  void GetSplitPlaneOpt(SItemVec *vec, int axisToTest);
  void GetSplitPlaneOpt2(SItemVec *vec, int axisToTest);
  void GetSplitPlaneOpt3(SItemVec *vec, int axisToTest);
  void GetSplitPlaneOptUnroll4(SItemVec *vec, int axisToTest);
public:
  // setting the evaluation for split cases that must not be done
  float WorstEvaluation() const { return MAXFLOAT;}

  // update the best value for evaluation
  int UpdateEvaluation(float &eval, const float &newEval);

public:
  // default constructor
  CKTBABuildUp();

  // default destructor
  virtual ~CKTBABuildUp();

  // provide info about construction method
  virtual void ProvideID(ostream &app);  
  
  // constructs the KD-tree for given objectlist and given bounding box
  // returns NULL in case of failure, in case of success returns
  // the pointer to the root node of constructed KD-tree.
  virtual SKTBNodeT* BuildUp(ObjectContainer &objlist,
                             const AxisAlignedBox3 &box,
                             bool verbose = true);
};

}

#endif // __KTBAI_H__