source: orxonox.OLD/orxonox/trunk/src/lib/collision_detection/obb_tree_node.cc @ 4627

/*
   orxonox - the future of 3D-vertical-scrollers

   Copyright (C) 2004 orx

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2, or (at your option)
   any later version.

### File Specific:
   main-programmer: Patrick Boenzli
   co-programmer: ...
*/

#define DEBUG_SPECIAL_MODULE DEBUG_MODULE_COLLISION

#include "obb_tree_node.h"
#include "list.h"
#include "obb.h"
#include "obb_tree.h"
#include "vector.h"
#include "abstract_model.h"

#include <math.h>


#define WANT_STREAM
#define WANT_MATH
#define WANT_FSTREAM


#include "include.h"
#include "newmat.h"
#include "newmatap.h"
#include "newmatio.h"

#include "lin_alg.h"




using namespace std;

OBBTree*  OBBTreeNode::obbTree = NULL;

/**
   \brief standard constructor
 */
OBBTreeNode::OBBTreeNode ()
{
  this->setClassID(CL_OBB_TREE_NODE, "OBBTreeNode");
  this->nodeLeft = NULL;
  this->nodeRight = NULL;
}


/**
   \brief standard deconstructor
 */
OBBTreeNode::~OBBTreeNode ()
{
  // delete what has to be deleted here
}



/**
   \brief creates a new BVTree or BVTree partition
   \param depth: how much more depth-steps to go: if == 1 don't go any deeper!
   \param verticesList: the list of vertices of the object - each vertices triple is interpreted as a triangle
 */
void OBBTreeNode::spawnBVTree(const int depth, sVec3D *verticesList, const int length)
{
  printf("OBB Depth: %i\n", depth);
  this->depth = depth;

  this->bvElement = this->createBox();
  PRINTF(0)("Created OBBox\n");
  this->calculateBoxAttributes(this->bvElement, verticesList, length);
  PRINTF(0)("Calculated attributes\n");

  if( likely( this->depth > 0))
  {
    this->forkBox(this->bvElement);

  }
}


OBB* OBBTreeNode::createBox()
{
  return new OBB();
}


void OBBTreeNode::calculateBoxAttributes(OBB* box, sVec3D* verticesList, int length)
{
  float     facelet[length];                         //!< surface area of the i'th triangle of the convex hull
  float     face;                                    //!< surface area of the entire convex hull
  Vector    centroid[length];                        //!< centroid of the i'th convex hull
  Vector    center;                                  //!< the center of the entire hull
  Vector    p, q, r;                                 //!< holder of the polygon data, much more conveniant to work with Vector than sVec3d
  Vector    t1, t2;                                  //!< temporary values
  double    covariance[3][3];                        //!< the covariance matrix

  this->numOfVertices = length;
  this->vertices = verticesList;
  box->vertices = verticesList;
  box->numOfVertices = length;


  /* fist compute all the convex hull face/facelets and centroids */
  for(int i = 0; i < length; i+=3)          /* FIX-ME-QUICK: hops of 3, array indiscontinuity*/
  {
    p = verticesList[i];
    q = verticesList[i +1];
    r = verticesList[i + 2];

    t1 = p - q; t2 = p - r;

    /* finding the facelet surface via cross-product */
    facelet[i] = 0.5f * fabs( t1.cross(t2).len() );
    /* update the entire convex hull surface */
    face += facelet[i];

    /* calculate the cetroid of the hull triangles */
    centroid[i] = (p + q + r) * 1/3;
    /* now calculate the centroid of the entire convex hull, weighted average of triangle centroids */
    center += centroid[i] * facelet[i];
  }
  /* take the average of the centroid sum */
  center /= face;
  PRINTF(0)("-- Calculated Center\n");


  /* now calculate the covariance matrix - if not written in three for-loops, it would compute faster: minor */
  for(int j = 0; j < 3; ++j)
  {
    for(int k = 0; k < 3; ++k)
    {
      for(int i = 0; i < length; i+=3)
      {
        p = verticesList[i];
        q = verticesList[i +1];
        r = verticesList[i + 2];

        covariance[j][k] = facelet[i] / (12.0f * face) * (9.0f * centroid[i][j] * centroid[i][k] + p[j]* p[k] +
            q[j] * q[k] + r[j]*r[k]) - center[j] * center[k];
      }
    }
  }
  PRINTF(0)("-- Calculated Covariance\n");

//   printf("\nVertex Data:\n");
//   for(int i = 0; i < length; i++)
//   {
//     printf("vertex %i: %f, %f, %f\n", i, verticesList[i][0], verticesList[i][1], verticesList[i][2]);
//   }

  printf("\nCovariance Matrix:\n");
  for(int j = 0; j < 3; ++j)
  {
    printf(" |");
    for(int k = 0; k < 3; ++k)
    {
      printf(" \b%f ", covariance[j][k]);
    }
    printf(" |\n");
  }
  printf("center: %f, %f, %f\n\n", center.x, center.y, center.z);


  for(int i = 0; i < 3; ++i)
  {

    box->covarianceMatrix[i][0] = covariance[i][0];
    box->covarianceMatrix[i][1] = covariance[i][1];
    box->covarianceMatrix[i][3] = covariance[i][2];
  }
  *box->center = center;
  PRINTF(0)("-- Written Result to obb\n");

  /* now getting spanning vectors of the sub-space:
  the eigenvectors of a symmertric matrix, such as the
  covarience matrix are mutually orthogonal.
  after normalizing them, they can be used as a the basis
  vectors
  */
  Matrix                V(3,3);                               //!< for eigenvectors
  DiagonalMatrix        D(3);                                 //!< for eigenvalues
  SymmetricMatrix       C(3);                                 //!< for the covariance symmetrical matrix
  Vector**              axis = new Vector*[3];                //!< the references to the obb axis


  double a[4][4];

  a[0][0] = C(1,1) = covariance[0][0];
  a[0][1] = C(1,2) = covariance[0][1];
  a[0][2] = C(1,3) = covariance[0][2];
  a[1][0] = C(2,1) = covariance[1][0];
  a[1][1] = C(2,2) = covariance[1][1];
  a[1][2] = C(2,3) = covariance[1][2];
  a[2][0] = C(3,1) = covariance[2][0];
  a[2][1] = C(3,2) = covariance[2][1];
  a[2][2] = C(3,3) = covariance[2][2];

  Jacobi(C, D, V);                                            /* do the jacobi decomposition */
  PRINTF(0)("-- Done Jacobi Decomposition\n");


  /* new jacobi tests */
  double eigenvectors[3][3];
  double eigval[3];

  EVJacobi jac;
  jac.setMatrix(2, covariance, 0, 0);
  jac.getEigenVector(eigenvectors);


  printf("Old Jacobi\n");
  for(int j = 1; j < 4; ++j)
  {
    printf(" |");
    for(int k = 1; k < 4; ++k)
    {
      printf(" \b%f ", V(j, k));
    }
    printf(" |\n");
  }

  printf("New Jacobi\n");
  for(int j = 0; j < 3; ++j)
  {
    printf(" |");
    for(int k = 0; k < 3; ++k)
    {
      printf(" \b%f ", eigenvectors[j][k]);
    }
    printf(" |\n");
  }

  axis[0] = new Vector(V(1, 1), V(2, 1), V(3, 1));
  axis[1] = new Vector(V(1, 2), V(2, 2), V(3, 2));
  axis[2] = new Vector(V(1, 3), V(2, 3), V(3, 3));
  box->axis = axis;
  PRINTF(0)("-- Got Axis\n");

//   delete &V;
//   delete &D;
//   delete &V;

//   printf("\neigenvector: %f, %f, %f\n", box->axis[0]->x, box->axis[0]->y, box->axis[0]->z);
//   printf("eigenvector: %f, %f, %f\n", box->axis[1]->x, box->axis[1]->y, box->axis[1]->z);
//   printf("eigenvector: %f, %f, %f\n", box->axis[2]->x, box->axis[2]->y, box->axis[2]->z);


  /* now get the axis length */
  Line                ax[3];                                 //!< the axis
  float*              halfLength = new float[3];             //!< half length of the axis
  float               tmpLength;                             //!< tmp save point for the length
  Plane               p0(*box->axis[0], *box->center);       //!< the axis planes
  Plane               p1(*box->axis[1], *box->center);
  Plane               p2(*box->axis[2], *box->center);

  halfLength[0] = -1.0f;
  for(int j = 0; j < length; ++j)
  {
    tmpLength = fabs(p0.distancePoint(vertices[j]));
    if( tmpLength > halfLength[0])
      halfLength[0] = tmpLength;
  }

  halfLength[1] = -1.0f;
  for(int j = 0; j < length; ++j)
  {
    tmpLength = fabs(p1.distancePoint(vertices[j]));
    if( tmpLength > halfLength[1])
      halfLength[1] = tmpLength;
  }

  halfLength[2] = -1.0f;
  for(int j = 0; j < length; ++j)
  {
    tmpLength = fabs(p2.distancePoint(vertices[j]));
    if( tmpLength > halfLength[2])
      halfLength[2] = tmpLength;
  }

  box->halfLength = halfLength;
  PRINTF(0)("-- Written Axis to obb\n");


//   printf("\nwe got length: \n");
//   for(int i = 0; i < 3; ++i)
//     printf("length[%i] = %f\n", i, box->halfLength[i]);
}



/**
  \brief this separates an ob-box in the middle
  \param box: the box to separate

  this will separate the box into to smaller boxes. the separation is done along the middle of the longest axis
 */
void OBBTreeNode::forkBox(OBB* box)
{
  /* get the longest axis of the box */
  float               aLength = -1.0f;                     //!< the length of the longest axis
  int                 axisIndex = 0;                       //!< this is the nr of the longest axis

  for(int i = 0; i < 3; ++i)
  {
    if( aLength < box->halfLength[i])
    {
      aLength = box->halfLength[i];
      axisIndex = i;
    }
  }

  printf("\nlongest axis is: nr %i with a half-length of: %f\n", axisIndex, aLength);


  /* get the closest vertex near the center */
  float               dist = 999999.0f;                    //!< the smallest distance to each vertex
  float               tmpDist;                             //!< temporary distance
  int                 vertexIndex;
  Plane               middlePlane(*box->axis[axisIndex], *box->center); //!< the middle plane

  for(int i = 0; i < box->numOfVertices; ++i)
  {
    tmpDist = fabs(middlePlane.distancePoint(box->vertices[i]));
    if( tmpDist < dist)
    {
      dist = tmpDist;
      vertexIndex = i;
    }
  }

  printf("\nthe clostest vertex is nr: %i, with a dist of: %f\n", vertexIndex ,dist);


  /* now definin the separation plane through this specified nearest point and partition
  the points depending on which side they are located
  */
  Plane              separationPlane(*box->axis[axisIndex], box->vertices[vertexIndex]);  //!< separation plane
  tList<sVec3D>      partition1;                           //!< the vertex partition 1
  tList<sVec3D>      partition2;                           //!< the vertex partition 2

  for(int i = 0; i < box->numOfVertices; ++i)
  {
    if( separationPlane.distancePoint(box->vertices[i]) > 0.0f)
      partition1.add(&box->vertices[i]);
    else
      partition2.add(&box->vertices[i]);
  }
  partition1.add(&box->vertices[vertexIndex]);

  printf("\npartition1: got %i vertices/ partition 2: got %i vertices\n", partition1.getSize(), partition2.getSize());


  /* now comes the separation into two different sVec3D arrays */
  tIterator<sVec3D>* iterator;                             //!< the iterator to go through the lists
  sVec3D*            element;                              //!< the elements
  int                index;                                //!< index storage place
  sVec3D*            vertList1;                            //!< the vertex list 1
  sVec3D*            vertList2;                            //!< the vertex list 2

  vertList1 = new sVec3D[partition1.getSize()];
  vertList2 = new sVec3D[partition2.getSize()];

  iterator = partition1.getIterator();
  element = iterator->nextElement();
  index = 0;
  while( element != NULL)
  {
    vertList1[index][0] = element[0][0];
    vertList1[index][1] = element[0][1];
    vertList1[index][2] = element[0][2];
    ++index;
    element = iterator->nextElement();
  }

//   printf("\npartition 1:\n");
//   for(int i = 0; i < partition1.getSize(); ++i)
//   {
//     printf("v[%i][0] = %f\n", i, vertList1[i][0]);
//     printf("v[%i][1] = %f\n", i, vertList1[i][1]);
//     printf("v[%i][2] = %f\n", i, vertList1[i][2]);
//   }

  iterator = partition2.getIterator();
  element = iterator->nextElement();
  index = 0;
  while( element != NULL)
  {
    vertList2[index][0] = element[0][0];
    vertList2[index][1] = element[0][1];
    vertList2[index][2] = element[0][2];
    ++index;
    element = iterator->nextElement();
  }

  //delete iterator;
//   printf("\npartition 2:\n");
//   for(int i = 0; i < partition2.getSize(); ++i)
//   {
//     printf("v[%i][0] = %f\n", i, vertList2[i][0]);
//     printf("v[%i][1] = %f\n", i, vertList2[i][1]);
//     printf("v[%i][2] = %f\n", i, vertList2[i][2]);
//   }

  /* now spawn the obb tree: create the nodes and descent */
  OBBTreeNode*       node1 = new OBBTreeNode();
  OBBTreeNode*       node2 = new OBBTreeNode();

  this->nodeLeft = node1;
  this->nodeRight = node2;

  this->nodeLeft->spawnBVTree(depth - 1, vertList1, partition1.getSize());
  this->nodeRight->spawnBVTree(depth - 1, vertList2, partition2.getSize());
}




void OBBTreeNode::collideWith(const BVTree &tree)
{}




void OBBTreeNode::drawBV(int depth) const
{
  glBegin(GL_TRIANGLES);
  glColor3f(1.0, 1.0, 1.0);
  for(int i = 0; i < this->bvElement->numOfVertices; ++i)
    {
      glVertex3f(this->bvElement->vertices[i][0], this->bvElement->vertices[i][1], this->bvElement->vertices[i][2]);
      //printf("v(%f, %f, %f)\n", this->vertices[i][0], this->vertices[i][1], this->vertices[i][2]);
    }
  glEnd();
  //this->drawBVPolygon(depth);
}


void OBBTreeNode::drawBVPolygon(int depth) const
{
  //OBBTree::material->select();

  this->obbTree->getMaterial(depth)->select();

  /* draw world axes */
//   glBegin(GL_LINES);
//   glColor3f(0.0, 0.4, 0.3);
//   glVertex3f(0.0, 0.0, 0.0);
//   glVertex3f(3.0, 0.0, 0.0);
//
//   glVertex3f(0.0, 0.0, 0.0);
//   glVertex3f(0.0, 3.0, 0.0);
//
//   glVertex3f(0.0, 0.0, 0.0);
//   glVertex3f(0.0, 0.0, 3.0);
//   glEnd();



  /* draw the obb axes */
//   glBegin(GL_LINES);
//   glColor3f(0.0, 0.4, 0.3);
//   glVertex3f(this->bvElement->center->x, this->bvElement->center->y, this->bvElement->center->z);
//   glVertex3f(this->bvElement->center->x + this->bvElement->axis[0]->x * this->bvElement->halfLength[0],
//              this->bvElement->center->y + this->bvElement->axis[0]->y * this->bvElement->halfLength[0],
//              this->bvElement->center->z + this->bvElement->axis[0]->z * this->bvElement->halfLength[0]);
//
//   glVertex3f(this->bvElement->center->x, this->bvElement->center->y, this->bvElement->center->z);
//   glVertex3f(this->bvElement->center->x + this->bvElement->axis[1]->x * this->bvElement->halfLength[1],
//              this->bvElement->center->y + this->bvElement->axis[1]->y * this->bvElement->halfLength[1],
//              this->bvElement->center->z + this->bvElement->axis[1]->z * this->bvElement->halfLength[1]);
//
//   glVertex3f(this->bvElement->center->x, this->bvElement->center->y, this->bvElement->center->z);
//   glVertex3f(this->bvElement->center->x + this->bvElement->axis[2]->x * this->bvElement->halfLength[2],
//              this->bvElement->center->y + this->bvElement->axis[2]->y * this->bvElement->halfLength[2],
//              this->bvElement->center->z + this->bvElement->axis[2]->z * this->bvElement->halfLength[2]);
//   glEnd();


  Vector cen = *this->bvElement->center;
  Vector** axis = this->bvElement->axis;
  float* len = this->bvElement->halfLength;

  /* draw bounding box */
  glBegin(GL_LINE_LOOP);
  glColor3f(0.3, 0.4, 0.7);
  glVertex3f(cen.x + axis[0]->x * len[0] + axis[1]->x * len[1] + axis[2]->x * len[2],
             cen.y + axis[0]->y * len[0] + axis[1]->y * len[1] + axis[2]->y * len[2],
             cen.z + axis[0]->z * len[0] + axis[1]->z * len[1] + axis[2]->z * len[2]);
  glVertex3f(cen.x + axis[0]->x * len[0] + axis[1]->x * len[1] - axis[2]->x * len[2],
             cen.y + axis[0]->y * len[0] + axis[1]->y * len[1] - axis[2]->y * len[2],
             cen.z + axis[0]->z * len[0] + axis[1]->z * len[1] - axis[2]->z * len[2]);
  glVertex3f(cen.x + axis[0]->x * len[0] - axis[1]->x * len[1] - axis[2]->x * len[2],
             cen.y + axis[0]->y * len[0] - axis[1]->y * len[1] - axis[2]->y * len[2],
             cen.z + axis[0]->z * len[0] - axis[1]->z * len[1] - axis[2]->z * len[2]);
  glVertex3f(cen.x + axis[0]->x * len[0] - axis[1]->x * len[1] + axis[2]->x * len[2],
             cen.y + axis[0]->y * len[0] - axis[1]->y * len[1] + axis[2]->y * len[2],
             cen.z + axis[0]->z * len[0] - axis[1]->z * len[1] + axis[2]->z * len[2]);
  glEnd();

  glBegin(GL_LINE_LOOP);
  glVertex3f(cen.x + axis[0]->x * len[0] - axis[1]->x * len[1] + axis[2]->x * len[2],
             cen.y + axis[0]->y * len[0] - axis[1]->y * len[1] + axis[2]->y * len[2],
             cen.z + axis[0]->z * len[0] - axis[1]->z * len[1] + axis[2]->z * len[2]);
  glVertex3f(cen.x + axis[0]->x * len[0] - axis[1]->x * len[1] - axis[2]->x * len[2],
             cen.y + axis[0]->y * len[0] - axis[1]->y * len[1] - axis[2]->y * len[2],
             cen.z + axis[0]->z * len[0] - axis[1]->z * len[1] - axis[2]->z * len[2]);
  glVertex3f(cen.x - axis[0]->x * len[0] - axis[1]->x * len[1] - axis[2]->x * len[2],
             cen.y - axis[0]->y * len[0] - axis[1]->y * len[1] - axis[2]->y * len[2],
             cen.z - axis[0]->z * len[0] - axis[1]->z * len[1] - axis[2]->z * len[2]);
  glVertex3f(cen.x - axis[0]->x * len[0] - axis[1]->x * len[1] + axis[2]->x * len[2],
             cen.y - axis[0]->y * len[0] - axis[1]->y * len[1] + axis[2]->y * len[2],
             cen.z - axis[0]->z * len[0] - axis[1]->z * len[1] + axis[2]->z * len[2]);
  glEnd();

  glBegin(GL_LINE_LOOP);
  glVertex3f(cen.x - axis[0]->x * len[0] - axis[1]->x * len[1] + axis[2]->x * len[2],
             cen.y - axis[0]->y * len[0] - axis[1]->y * len[1] + axis[2]->y * len[2],
             cen.z - axis[0]->z * len[0] - axis[1]->z * len[1] + axis[2]->z * len[2]);
  glVertex3f(cen.x - axis[0]->x * len[0] - axis[1]->x * len[1] - axis[2]->x * len[2],
             cen.y - axis[0]->y * len[0] - axis[1]->y * len[1] - axis[2]->y * len[2],
             cen.z - axis[0]->z * len[0] - axis[1]->z * len[1] - axis[2]->z * len[2]);
  glVertex3f(cen.x - axis[0]->x * len[0] + axis[1]->x * len[1] - axis[2]->x * len[2],
             cen.y - axis[0]->y * len[0] + axis[1]->y * len[1] - axis[2]->y * len[2],
             cen.z - axis[0]->z * len[0] + axis[1]->z * len[1] - axis[2]->z * len[2]);
  glVertex3f(cen.x - axis[0]->x * len[0] + axis[1]->x * len[1] + axis[2]->x * len[2],
             cen.y - axis[0]->y * len[0] + axis[1]->y * len[1] + axis[2]->y * len[2],
             cen.z - axis[0]->z * len[0] + axis[1]->z * len[1] + axis[2]->z * len[2]);
  glEnd();

  glBegin(GL_LINE_LOOP);
  glVertex3f(cen.x - axis[0]->x * len[0] + axis[1]->x * len[1] - axis[2]->x * len[2],
             cen.y - axis[0]->y * len[0] + axis[1]->y * len[1] - axis[2]->y * len[2],
             cen.z - axis[0]->z * len[0] + axis[1]->z * len[1] - axis[2]->z * len[2]);
  glVertex3f(cen.x - axis[0]->x * len[0] + axis[1]->x * len[1] + axis[2]->x * len[2],
             cen.y - axis[0]->y * len[0] + axis[1]->y * len[1] + axis[2]->y * len[2],
             cen.z - axis[0]->z * len[0] + axis[1]->z * len[1] + axis[2]->z * len[2]);
  glVertex3f(cen.x + axis[0]->x * len[0] + axis[1]->x * len[1] + axis[2]->x * len[2],
             cen.y + axis[0]->y * len[0] + axis[1]->y * len[1] + axis[2]->y * len[2],
             cen.z + axis[0]->z * len[0] + axis[1]->z * len[1] + axis[2]->z * len[2]);
  glVertex3f(cen.x + axis[0]->x * len[0] + axis[1]->x * len[1] - axis[2]->x * len[2],
             cen.y + axis[0]->y * len[0] + axis[1]->y * len[1] - axis[2]->y * len[2],
             cen.z + axis[0]->z * len[0] + axis[1]->z * len[1] - axis[2]->z * len[2]);
  glEnd();

/*
  glVertex3f(cen.x - axis[0]->x * len[0] + axis[1]->x * len[1] - axis[2]->x * len[2],
  cen.y - axis[0]->y * len[0] + axis[1]->y * len[1] - axis[2]->y * len[2],
  cen.z - axis[0]->z * len[0] + axis[1]->z * len[1] - axis[2]->z * len[2]);
  glVertex3f(cen.x - axis[0]->x * len[0] + axis[1]->x * len[1] + axis[2]->x * len[2],
  cen.y - axis[0]->y * len[0] + axis[1]->y * len[1] + axis[2]->y * len[2],
  cen.z - axis[0]->z * len[0] + axis[1]->z * len[1] + axis[2]->z * len[2]);*/


  glEnd();

  if( this->nodeLeft != NULL && depth != 0 )
    this->nodeLeft->drawBVPolygon(depth - 1);
  if( this->nodeRight != NULL && depth != 0)
    this->nodeRight->drawBVPolygon(depth - 1);

}


void OBBTreeNode::drawBVBlended(int depth) const
{}


void OBBTreeNode::debug()
{

  /*
  for(int i = 0; i < length; i++)
  {
  printf("vertex %i: %f, %f, %f\n", i, verticesList[i][0], verticesList[i][1], verticesList[i][2]);
}
  */
}