svn merge -r 14721:14810 https://svn.blender.org/svnroot/bf-blender/trunk/blender
authorDaniel Genrich <daniel.genrich@gmx.net>
Mon, 12 May 2008 12:24:52 +0000 (12:24 +0000)
committerDaniel Genrich <daniel.genrich@gmx.net>
Mon, 12 May 2008 12:24:52 +0000 (12:24 +0000)
1  2 
source/Makefile
source/blender/blenkernel/intern/collision.c
source/blender/blenkernel/intern/modifier.c
source/blender/blenkernel/intern/particle_system.c
source/blender/blenlib/BLI_kdopbvh.h
source/blender/blenlib/intern/BLI_kdopbvh.c
source/blender/blenloader/intern/readfile.c
source/blender/src/buttons_editing.c
source/blender/src/buttons_object.c
source/blender/src/drawobject.c

diff --cc source/Makefile
Simple merge
index f3637b4dda2e32d64379873f8a64f6ff962933f1,e244ccca306304b80feb18842c8a7bd5603e41c2..9ba47874d3cfbc86491cb6d6cebc9c45bd7ed068
@@@ -183,42 -130,42 +183,43 @@@ Collision modifier code en
   * copied from SOLVE_CUBIC.C --> GSL
   */
  
--/* DG: debug hint! don't forget that all functions were "fabs", "sinf", etc before */
--#define mySWAP(a,b) { float tmp = b ; b = a ; a = tmp ; }
++#define mySWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0)
  
--int gsl_poly_solve_cubic ( float a, float b, float c, float *x0, float *x1, float *x2 )
++int 
++              gsl_poly_solve_cubic (double a, double b, double c, 
++                                                        double *x0, double *x1, double *x2)
  {
--      float q = ( a * a - 3 * b );
--      float r = ( 2 * a * a * a - 9 * a * b + 27 * c );
++      double q = (a * a - 3 * b);
++      double r = (2 * a * a * a - 9 * a * b + 27 * c);
  
--      float Q = q / 9;
--      float R = r / 54;
++      double Q = q / 9;
++      double R = r / 54;
  
--      float Q3 = Q * Q * Q;
--      float R2 = R * R;
++      double Q3 = Q * Q * Q;
++      double R2 = R * R;
  
--      float CR2 = 729 * r * r;
--      float CQ3 = 2916 * q * q * q;
++      double CR2 = 729 * r * r;
++      double CQ3 = 2916 * q * q * q;
  
--      if ( R == 0 && Q == 0 )
++      if (R == 0 && Q == 0)
        {
                *x0 = - a / 3 ;
                *x1 = - a / 3 ;
                *x2 = - a / 3 ;
                return 3 ;
        }
--      else if ( CR2 == CQ3 )
++      else if (CR2 == CQ3) 
        {
--              /* this test is actually R2 == Q3, written in a form suitable
 -                for exact computation with integers */
++      /* this test is actually R2 == Q3, written in a form suitable
 +              for exact computation with integers */
  
--              /* Due to finite precision some float roots may be missed, and
 -                considered to be a pair of complex roots z = x +/- epsilon i
 -                close to the real axis. */
++      /* Due to finite precision some double roots may be missed, and
 +              considered to be a pair of complex roots z = x +/- epsilon i
 +              close to the real axis. */
  
--              float sqrtQ = sqrt ( Q );
++              double sqrtQ = sqrt (Q);
  
--              if ( R > 0 )
++              if (R > 0)
                {
                        *x0 = -2 * sqrtQ  - a / 3;
                        *x1 = sqrtQ - a / 3;
                }
                return 3 ;
        }
--      else if ( CR2 < CQ3 ) /* equivalent to R2 < Q3 */
++      else if (CR2 < CQ3) /* equivalent to R2 < Q3 */
        {
--              float sqrtQ = sqrt ( Q );
--              float sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
--              float theta = acos ( R / sqrtQ3 );
--              float norm = -2 * sqrtQ;
--              *x0 = norm * cos ( theta / 3 ) - a / 3;
--              *x1 = norm * cos ( ( theta + 2.0 * M_PI ) / 3 ) - a / 3;
--              *x2 = norm * cos ( ( theta - 2.0 * M_PI ) / 3 ) - a / 3;
--
++              double sqrtQ = sqrt (Q);
++              double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
++              double theta = acos (R / sqrtQ3);
++              double norm = -2 * sqrtQ;
++              *x0 = norm * cos (theta / 3) - a / 3;
++              *x1 = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3;
++              *x2 = norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3;
++      
                /* Sort *x0, *x1, *x2 into increasing order */
  
--              if ( *x0 > *x1 )
--                      mySWAP ( *x0, *x1 ) ;
--
--              if ( *x1 > *x2 )
++              if (*x0 > *x1)
++                      mySWAP(*x0, *x1) ;
++      
++              if (*x1 > *x2)
                {
--                      mySWAP ( *x1, *x2 ) ;
--
--                      if ( *x0 > *x1 )
--                              mySWAP ( *x0, *x1 ) ;
++                      mySWAP(*x1, *x2) ;
++          
++                      if (*x0 > *x1)
++                              mySWAP(*x0, *x1) ;
                }
--
++      
                return 3;
        }
        else
        {
--              float sgnR = ( R >= 0 ? 1 : -1 );
--              float A = -sgnR * pow ( ABS ( R ) + sqrt ( R2 - Q3 ), 1.0/3.0 );
--              float B = Q / A ;
++              double sgnR = (R >= 0 ? 1 : -1);
++              double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0/3.0);
++              double B = Q / A ;
                *x0 = A + B - a / 3;
                return 1;
        }
  }
  
  
++
  /**
   * gsl_poly_solve_quadratic
   *
   * copied from GSL
   */
--int gsl_poly_solve_quadratic ( float a, float b, float c,  float *x0, float *x1 )
++int 
++              gsl_poly_solve_quadratic (double a, double b, double c, 
++                                                                double *x0, double *x1)
  {
--      float disc = b * b - 4 * a * c;
++      double disc = b * b - 4 * a * c;
++
++      if (a == 0) /* Handle linear case */
++      {
++              if (b == 0)
++              {
++                      return 0;
++              }
++              else
++              {
++                      *x0 = -c / b;
++                      return 1;
++              };
++      }
  
--      if ( disc > 0 )
++      if (disc > 0)
        {
--              if ( b == 0 )
++              if (b == 0)
                {
--                      float r = ABS ( 0.5 * sqrt ( disc ) / a );
++                      double r = fabs (0.5 * sqrt (disc) / a);
                        *x0 = -r;
                        *x1 =  r;
                }
                else
                {
--                      float sgnb = ( b > 0 ? 1 : -1 );
--                      float temp = -0.5 * ( b + sgnb * sqrt ( disc ) );
--                      float r1 = temp / a ;
--                      float r2 = c / temp ;
++                      double sgnb = (b > 0 ? 1 : -1);
++                      double temp = -0.5 * (b + sgnb * sqrt (disc));
++                      double r1 = temp / a ;
++                      double r2 = c / temp ;
  
--                      if ( r1 < r2 )
++                      if (r1 < r2) 
                        {
                                *x0 = r1 ;
                                *x1 = r2 ;
--                      }
--                      else
++                      } 
++                      else 
                        {
                                *x0 = r2 ;
                                *x1 = r1 ;
                }
                return 2;
        }
--      else if ( disc == 0 )
++      else if (disc == 0) 
        {
                *x0 = -0.5 * b / a ;
                *x1 = -0.5 * b / a ;
  
  
  
++
  /*
   * See Bridson et al. "Robust Treatment of Collision, Contact and Friction for Cloth Animation"
   *     page 4, left column
   */
--
--int cloth_get_collision_time ( float a[3], float b[3], float c[3], float d[3], float e[3], float f[3], float solution[3] )
++int cloth_get_collision_time ( double a[3], double b[3], double c[3], double d[3], double e[3], double f[3], double solution[3] )
  {
        int num_sols = 0;
  
--      float g = -a[2] * c[1] * e[0] + a[1] * c[2] * e[0] +
--                a[2] * c[0] * e[1] - a[0] * c[2] * e[1] -
--                a[1] * c[0] * e[2] + a[0] * c[1] * e[2];
--
--      float h = -b[2] * c[1] * e[0] + b[1] * c[2] * e[0] - a[2] * d[1] * e[0] +
--                a[1] * d[2] * e[0] + b[2] * c[0] * e[1] - b[0] * c[2] * e[1] +
--                a[2] * d[0] * e[1] - a[0] * d[2] * e[1] - b[1] * c[0] * e[2] +
--                b[0] * c[1] * e[2] - a[1] * d[0] * e[2] + a[0] * d[1] * e[2] -
--                a[2] * c[1] * f[0] + a[1] * c[2] * f[0] + a[2] * c[0] * f[1] -
--                a[0] * c[2] * f[1] - a[1] * c[0] * f[2] + a[0] * c[1] * f[2];
--
--      float i = -b[2] * d[1] * e[0] + b[1] * d[2] * e[0] +
--                b[2] * d[0] * e[1] - b[0] * d[2] * e[1] -
--                b[1] * d[0] * e[2] + b[0] * d[1] * e[2] -
--                b[2] * c[1] * f[0] + b[1] * c[2] * f[0] -
--                a[2] * d[1] * f[0] + a[1] * d[2] * f[0] +
--                b[2] * c[0] * f[1] - b[0] * c[2] * f[1] +
--                a[2] * d[0] * f[1] - a[0] * d[2] * f[1] -
--                b[1] * c[0] * f[2] + b[0] * c[1] * f[2] -
--                a[1] * d[0] * f[2] + a[0] * d[1] * f[2];
--
--      float j = -b[2] * d[1] * f[0] + b[1] * d[2] * f[0] +
--                b[2] * d[0] * f[1] - b[0] * d[2] * f[1] -
--                b[1] * d[0] * f[2] + b[0] * d[1] * f[2];
++      // x^0 - checked 
++      double g =      a[0] * c[1] * e[2] - a[0] * c[2] * e[1] +
++                              a[1] * c[2] * e[0] - a[1] * c[0] * e[2] + 
++                              a[2] * c[0] * e[1] - a[2] * c[1] * e[0];
++      
++      // x^1
++      double h = -b[2] * c[1] * e[0] + b[1] * c[2] * e[0] - a[2] * d[1] * e[0] +
++                      a[1] * d[2] * e[0] + b[2] * c[0] * e[1] - b[0] * c[2] * e[1] +
++                      a[2] * d[0] * e[1] - a[0] * d[2] * e[1] - b[1] * c[0] * e[2] +
++                      b[0] * c[1] * e[2] - a[1] * d[0] * e[2] + a[0] * d[1] * e[2] -
++                      a[2] * c[1] * f[0] + a[1] * c[2] * f[0] + a[2] * c[0] * f[1] -
++                      a[0] * c[2] * f[1] - a[1] * c[0] * f[2] + a[0] * c[1] * f[2];
++
++      // x^2
++      double i = -b[2] * d[1] * e[0] + b[1] * d[2] * e[0] +
++                      b[2] * d[0] * e[1] - b[0] * d[2] * e[1] -
++                      b[1] * d[0] * e[2] + b[0] * d[1] * e[2] -
++                      b[2] * c[1] * f[0] + b[1] * c[2] * f[0] -
++                      a[2] * d[1] * f[0] + a[1] * d[2] * f[0] +
++                      b[2] * c[0] * f[1] - b[0] * c[2] * f[1] + 
++                      a[2] * d[0] * f[1] - a[0] * d[2] * f[1] -
++                      b[1] * c[0] * f[2] + b[0] * c[1] * f[2] -
++                      a[1] * d[0] * f[2] + a[0] * d[1] * f[2];
++      
++      // x^3 - checked
++      double j = -b[2] * d[1] * f[0] + b[1] * d[2] * f[0] +
++                      b[2] * d[0] * f[1] - b[0] * d[2] * f[1] -
++                      b[1] * d[0] * f[2] + b[0] * d[1] * f[2];
++      
++      /*
++      printf("r1: %lf\n", a[0] * c[1] * e[2] - a[0] * c[2] * e[1]);
++      printf("r2: %lf\n", a[1] * c[2] * e[0] - a[1] * c[0] * e[2]);
++      printf("r3: %lf\n", a[2] * c[0] * e[1] - a[2] * c[1] * e[0]);
++      
++      printf("x1 x: %f, y: %f, z: %f\n", a[0], a[1], a[2]);
++      printf("x2 x: %f, y: %f, z: %f\n", c[0], c[1], c[2]);
++      printf("x3 x: %f, y: %f, z: %f\n", e[0], e[1], e[2]);
++      
++      printf("v1 x: %f, y: %f, z: %f\n", b[0], b[1], b[2]);
++      printf("v2 x: %f, y: %f, z: %f\n", d[0], d[1], d[2]);
++      printf("v3 x: %f, y: %f, z: %f\n", f[0], f[1], f[2]);
++      
++      printf("t^3: %lf, t^2: %lf, t^1: %lf, t^0: %lf\n", j, i, h, g);
++      */
  
        // Solve cubic equation to determine times t1, t2, t3, when the collision will occur.
--      if ( ABS ( j ) > ALMOST_ZERO )
++      if ( ABS ( j ) > DBL_EPSILON )
        {
                i /= j;
                h /= j;
                g /= j;
--
                num_sols = gsl_poly_solve_cubic ( i, h, g, &solution[0], &solution[1], &solution[2] );
        }
--      else if ( ABS ( i ) > ALMOST_ZERO )
++      else
        {
                num_sols = gsl_poly_solve_quadratic ( i, h, g, &solution[0], &solution[1] );
                solution[2] = -1.0;
        }
--      else if ( ABS ( h ) > ALMOST_ZERO )
--      {
--              solution[0] = -g / h;
--              solution[1] = solution[2] = -1.0;
--              num_sols = 1;
--      }
--      else if ( ABS ( g ) > ALMOST_ZERO )
--      {
--              solution[0] = 0;
--              solution[1] = solution[2] = -1.0;
--              num_sols = 1;
--      }
++
++      // printf("num_sols: %d, sol1: %lf, sol2: %lf, sol3: %lf\n", num_sols, solution[0],  solution[1],  solution[2]);
  
        // Discard negative solutions
--      if ( ( num_sols >= 1 ) && ( solution[0] < 0 ) )
++      if ( ( num_sols >= 1 ) && ( solution[0] < DBL_EPSILON ) )
        {
                --num_sols;
                solution[0] = solution[num_sols];
        }
--      if ( ( num_sols >= 2 ) && ( solution[1] < 0 ) )
++      if ( ( num_sols >= 2 ) && ( solution[1] < DBL_EPSILON ) )
        {
                --num_sols;
                solution[1] = solution[num_sols];
        }
--      if ( ( num_sols == 3 ) && ( solution[2] < 0 ) )
++      if ( ( num_sols == 3 ) && ( solution[2] < DBL_EPSILON ) )
        {
                --num_sols;
        }
@@@ -736,128 -703,137 +762,210 @@@ int cloth_are_edges_adjacent ( ClothMod
        return 0;
  }
  
- void cloth_collision_moving_edges ( ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair )
 -void cloth_collision_moving_edges ( ClothModifierData *clmd, ClothModifierData *coll_clmd, CollisionTree *tree1, CollisionTree *tree2 )
++int cloth_collision_moving_edges ( ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair )
  {
        EdgeCollPair edgecollpair;
 -      Cloth *cloth1=NULL, *cloth2=NULL;
 -      MFace *face1=NULL, *face2=NULL;
 -      ClothVertex *verts1=NULL, *verts2=NULL;
 +      Cloth *cloth1=NULL;
 +      ClothVertex *verts1=NULL;
        unsigned int i = 0, j = 0, k = 0;
        int numsolutions = 0;
--      float a[3], b[3], c[3], d[3], e[3], f[3], solution[3];
++      double x1[3], v1[3], x2[3], v2[3], x3[3], v3[3];
++      double solution[3];
 +      MVert *verts2 = collmd->current_x; // old x
 +      MVert *velocity2 = collmd->current_v; // velocity
-       float mintime = 0;
++      float mintime = FLT_MAX;
++      float distance;
++      float triA[3][3], triB[3][3];
++      int result = 0;
  
        cloth1 = clmd->clothObject;
 -      cloth2 = coll_clmd->clothObject;
 -
        verts1 = cloth1->verts;
 -      verts2 = cloth2->verts;
 -
 -      face1 = & ( cloth1->mfaces[tree1->tri_index] );
 -      face2 = & ( cloth2->mfaces[tree2->tri_index] );
 -
 -      for ( i = 0; i < 5; i++ )
++      
++      /*
++      double p[4][3] = {{0,0,0},{0,2,0},{1,1,-1},{1,1,1}};
++      double v[4][3] = {{0,0,0},{1,0,0},{-2,0,0},{-2,0,0}};
++      
++      double pp[2][3] = {{-1,-1,-1}, {2,2,2}};
++      
++      
++      VECSUB ( x1, p[1], p[0] );
++      VECSUB ( v1, v[1], v[0] );
++                      
++      VECSUB ( x2, p[2], p[0] );
++      VECSUB ( v2, v[2], v[0] );
++                      
++      VECSUB ( x3, p[3], p[0] );
++      VECSUB ( v3, v[3], v[0] );
++
++      printf("x1 x: %f, y: %f, z: %f\n", x1[0], x1[1], x1[2]);
++      printf("x2 x: %f, y: %f, z: %f\n", x2[0], x2[1], x2[2]);
++      printf("x3 x: %f, y: %f, z: %f\n", x3[0], x3[1], x3[2]);
++      
++      printf("v1 x: %f, y: %f, z: %f\n", v1[0], v1[1], v1[2]);
++      printf("v2 x: %f, y: %f, z: %f\n", v2[0], v2[1], v2[2]);
++      printf("v3 x: %f, y: %f, z: %f\n", v3[0], v3[1], v3[2]);
 +
++      numsolutions = cloth_get_collision_time ( x1, v1, x2, v2, x3, v3, solution );
++      
++      for ( k = 0; k < numsolutions; k++ )
++              printf("mintime: %f\n", solution[k]);
++      
++      mintime = solution[0];
++      
++      // move triangles to collision point in time
++      VECADDS(triA[0], pp[0], v[0], solution[0]);
++      VECADDS(triA[1], p[0], v[0], solution[0]);
++      VECADDS(triA[2], p[1], v[1], solution[0]);
++              
++      VECADDS(triB[0], pp[1], v[0], solution[0]);
++      VECADDS(triB[1], p[2], v[2], solution[0]);
++      VECADDS(triB[2], p[3], v[3], solution[0]);
++              
++              // check distance there
++      distance = plNearestPoints (triA[0], triA[1], triA[2], triB[0], triB[1], triB[2], collpair->pa,collpair->pb,collpair->vector );
++      
++      printf("mintime: %f, dist: %f\n", mintime, distance);
++      
++      exit(0);
++      */
 +      for(i = 0; i < 9; i++)
        {
 -              if ( i == 0 )
 +              // 9 edge - edge possibilities
 +              
 +              if(i == 0) // cloth edge: 1-2; coll edge: 1-2
                {
 -                      edgecollpair.p11 = face1->v1;
 -                      edgecollpair.p12 = face1->v2;
 +                      edgecollpair.p11 = collpair->ap1;
 +                      edgecollpair.p12 = collpair->ap2;
 +                      
 +                      edgecollpair.p21 = collpair->bp1;
 +                      edgecollpair.p22 = collpair->bp2;
                }
 -              else if ( i == 1 )
 +              else if(i == 1) // cloth edge: 1-2; coll edge: 2-3
                {
 -                      edgecollpair.p11 = face1->v2;
 -                      edgecollpair.p12 = face1->v3;
 +                      edgecollpair.p11 = collpair->ap1;
 +                      edgecollpair.p12 = collpair->ap2;
 +                      
 +                      edgecollpair.p21 = collpair->bp2;
 +                      edgecollpair.p22 = collpair->bp3;
                }
 -              else if ( i == 2 )
 +              else if(i == 2) // cloth edge: 1-2; coll edge: 1-3
                {
 -                      if ( face1->v4 )
 -                      {
 -                              edgecollpair.p11 = face1->v3;
 -                              edgecollpair.p12 = face1->v4;
 -                      }
 -                      else
 -                      {
 -                              edgecollpair.p11 = face1->v3;
 -                              edgecollpair.p12 = face1->v1;
 -                              i+=5; // get out of here after this edge pair is handled
 -                      }
 +                      edgecollpair.p11 = collpair->ap1;
 +                      edgecollpair.p12 = collpair->ap2;
 +                      
 +                      edgecollpair.p21 = collpair->bp1;
 +                      edgecollpair.p22 = collpair->bp3;
                }
 -              else if ( i == 3 )
 +              else if(i == 3) // cloth edge: 2-3; coll edge: 1-2
                {
 -                      if ( face1->v4 )
 -                      {
 -                              edgecollpair.p11 = face1->v4;
 -                              edgecollpair.p12 = face1->v1;
 -                      }
 -                      else
 -                              continue;
 +                      edgecollpair.p11 = collpair->ap2;
 +                      edgecollpair.p12 = collpair->ap3;
 +                      
 +                      edgecollpair.p21 = collpair->bp1;
 +                      edgecollpair.p22 = collpair->bp2;
                }
 -              else
 +              else if(i == 4) // cloth edge: 2-3; coll edge: 2-3
                {
 -                      edgecollpair.p11 = face1->v3;
 -                      edgecollpair.p12 = face1->v1;
 +                      edgecollpair.p11 = collpair->ap2;
 +                      edgecollpair.p12 = collpair->ap3;
 +                      
 +                      edgecollpair.p21 = collpair->bp2;
 +                      edgecollpair.p22 = collpair->bp3;
                }
 -
 -
 -              for ( j = 0; j < 5; j++ )
 +              else if(i == 5) // cloth edge: 2-3; coll edge: 1-3
                {
 -                      if ( j == 0 )
 -                      {
 -                              edgecollpair.p21 = face2->v1;
 -                              edgecollpair.p22 = face2->v2;
 -                      }
 -                      else if ( j == 1 )
 -                      {
 -                              edgecollpair.p21 = face2->v2;
 -                              edgecollpair.p22 = face2->v3;
 -                      }
 -                      else if ( j == 2 )
 -                      {
 -                              if ( face2->v4 )
 -                              {
 -                                      edgecollpair.p21 = face2->v3;
 -                                      edgecollpair.p22 = face2->v4;
 -                              }
 -                              else
 -                              {
 -                                      edgecollpair.p21 = face2->v3;
 -                                      edgecollpair.p22 = face2->v1;
 -                              }
 -                      }
 -                      else if ( j == 3 )
 -                      {
 -                              if ( face2->v4 )
 -                              {
 -                                      edgecollpair.p21 = face2->v4;
 -                                      edgecollpair.p22 = face2->v1;
 -                              }
 -                              else
 -                                      continue;
 -                      }
 -                      else
 -                      {
 -                              edgecollpair.p21 = face2->v3;
 -                              edgecollpair.p22 = face2->v1;
 -                      }
 -
 -
 -                      if ( !cloth_are_edges_adjacent ( clmd, coll_clmd, &edgecollpair ) )
 +                      edgecollpair.p11 = collpair->ap2;
 +                      edgecollpair.p12 = collpair->ap3;
 +                      
 +                      edgecollpair.p21 = collpair->bp1;
 +                      edgecollpair.p22 = collpair->bp3;
 +              }
 +              else if(i ==6) // cloth edge: 1-3; coll edge: 1-2
 +              {
 +                      edgecollpair.p11 = collpair->ap1;
 +                      edgecollpair.p12 = collpair->ap3;
 +                      
 +                      edgecollpair.p21 = collpair->bp1;
 +                      edgecollpair.p22 = collpair->bp2;
 +              }
 +              else if(i ==7) // cloth edge: 1-3; coll edge: 2-3
 +              {
 +                      edgecollpair.p11 = collpair->ap1;
 +                      edgecollpair.p12 = collpair->ap3;
 +                      
 +                      edgecollpair.p21 = collpair->bp2;
 +                      edgecollpair.p22 = collpair->bp3;
 +              }
 +              else if(i == 8) // cloth edge: 1-3; coll edge: 1-3
 +              {
 +                      edgecollpair.p11 = collpair->ap1;
 +                      edgecollpair.p12 = collpair->ap3;
 +                      
 +                      edgecollpair.p21 = collpair->bp1;
 +                      edgecollpair.p22 = collpair->bp3;
 +              }
 +              
 +              if ( !cloth_are_edges_adjacent ( clmd, collmd, &edgecollpair ) )
 +              {
 +                      // always put coll points in p21/p22
-                       VECSUB ( a, verts1[edgecollpair.p12].txold, verts1[edgecollpair.p11].txold );
-                       VECSUB ( b, verts1[edgecollpair.p12].tv, verts1[edgecollpair.p11].tv );
-                       VECSUB ( c, verts2[edgecollpair.p21].co, verts1[edgecollpair.p11].txold );
-                       VECSUB ( d, velocity2[edgecollpair.p21].co, verts1[edgecollpair.p11].tv );
-                       VECSUB ( e, verts2[edgecollpair.p22].co, verts1[edgecollpair.p11].txold );
-                       VECSUB ( f, velocity2[edgecollpair.p22].co, verts1[edgecollpair.p11].v );
-       
-                       numsolutions = cloth_get_collision_time ( a, b, c, d, e, f, solution );
++                      VECSUB ( x1, verts1[edgecollpair.p12].txold, verts1[edgecollpair.p11].txold );
++                      VECSUB ( v1, verts1[edgecollpair.p12].tv, verts1[edgecollpair.p11].tv );
++                      
++                      VECSUB ( x2, verts2[edgecollpair.p21].co, verts1[edgecollpair.p11].txold );
++                      VECSUB ( v2, velocity2[edgecollpair.p21].co, verts1[edgecollpair.p11].tv );
++                      
++                      VECSUB ( x3, verts2[edgecollpair.p22].co, verts1[edgecollpair.p11].txold );
++                      VECSUB ( v3, velocity2[edgecollpair.p22].co, verts1[edgecollpair.p11].tv );
++                      /*
++                      printf("A x: %f, y: %f, z: %f\n", a[0], a[1], a[2]);
++                      printf("B x: %f, y: %f, z: %f\n", b[0], b[1], b[2]);
++                      printf("C x: %f, y: %f, z: %f\n", c[0], c[1], c[2]);
++                      printf("D x: %f, y: %f, z: %f\n", d[0], d[1], d[2]);
++                      printf("E x: %f, y: %f, z: %f\n", e[0], e[1], e[2]);
++                      printf("F x: %f, y: %f, z: %f\n", f[0], f[1], f[2]);
++                      exit(0);
++                      */
++                      numsolutions = cloth_get_collision_time ( x1, v1, x2, v2, x3, v3, solution );
 +      
 +                      for ( k = 0; k < numsolutions; k++ )
                        {
-                               if ( ( solution[k] >= 0.0 ) && ( solution[k] <= 1.0 ) )
 -                              VECSUB ( a, verts1[edgecollpair.p12].xold, verts1[edgecollpair.p11].xold );
 -                              VECSUB ( b, verts1[edgecollpair.p12].v, verts1[edgecollpair.p11].v );
 -                              VECSUB ( c, verts1[edgecollpair.p21].xold, verts1[edgecollpair.p11].xold );
 -                              VECSUB ( d, verts1[edgecollpair.p21].v, verts1[edgecollpair.p11].v );
 -                              VECSUB ( e, verts2[edgecollpair.p22].xold, verts1[edgecollpair.p11].xold );
 -                              VECSUB ( f, verts2[edgecollpair.p22].v, verts1[edgecollpair.p11].v );
 -
 -                              numsolutions = cloth_get_collision_time ( a, b, c, d, e, f, solution );
 -
 -                              for ( k = 0; k < numsolutions; k++ )
++                              if ( ( solution[k] >= DBL_EPSILON ) && ( solution[k] <= 1.0 ) )
                                {
 -                                      if ( ( solution[k] >= 0.0 ) && ( solution[k] <= 1.0 ) )
 -                                      {
 -                                              //float out_collisionTime = solution[k];
 -
 -                                              // TODO: check for collisions
 -
 -                                              // TODO: put into (edge) collision list
 -
 -                                              // printf("Moving edge found!\n");
 -                                      }
 +                                      //float out_collisionTime = solution[k];
 +      
 +                                      // TODO: check for collisions
 +      
 +                                      // TODO: put into (edge) collision list
 +                                      
-                                       mintime = MIN2(mintime, solution[k]);
-       
-                                       printf("Moving edge found!, mintime: %f\n", mintime);
++                                      mintime = MIN2(mintime, (float)solution[k]);
++                                      
++//                                    printf("mt: %f, %lf, %f\n", mintime, solution[k], (float)solution[k]);
++                                      
++                                      result = 1;
 +                                      break;
                                }
                        }
                }
        }
++      
++      if(result)
++      {
++              // move triangles to collision point in time
++              VECADDS(triA[0], verts1[collpair->ap1].txold, verts1[collpair->ap1].tv, mintime);
++              VECADDS(triA[1], verts1[collpair->ap2].txold, verts1[collpair->ap2].tv, mintime);
++              VECADDS(triA[2], verts1[collpair->ap3].txold, verts1[collpair->ap3].tv, mintime);
++              
++              VECADDS(triB[0], collmd->current_x[collpair->bp1].co, collmd->current_v[collpair->bp1].co, mintime);
++              VECADDS(triB[1], collmd->current_x[collpair->bp2].co, collmd->current_v[collpair->bp2].co, mintime);
++              VECADDS(triB[2], collmd->current_x[collpair->bp3].co, collmd->current_v[collpair->bp3].co, mintime);
++              
++              // check distance there
++              distance = plNearestPoints (triA[0], triA[1], triA[2], triB[0], triB[1], triB[2], collpair->pa,collpair->pb,collpair->vector );
++              
++              printf("mintime: %f, dist: %f\n", mintime, distance);
++      }
++      
++      return result;
  }
  
  void cloth_collision_moving_tris ( ClothModifierData *clmd, ClothModifierData *coll_clmd, CollisionTree *tree1, CollisionTree *tree2 )
        ClothVertex *verts1=NULL, *verts2=NULL;
        unsigned int i = 0, j = 0, k = 0;
        int numsolutions = 0;
--      float a[3], b[3], c[3], d[3], e[3], f[3], solution[3];
++      float a[3], b[3], c[3], d[3], e[3], f[3];
++      double solution[3];
  
        for ( i = 0; i < 2; i++ )
        {
  
                                for ( k = 0; k < numsolutions; k++ )
                                {
--                                      if ( ( solution[k] >= 0.0 ) && ( solution[k] <= 1.0 ) )
++                                      if ( ( solution[k] >= ALMOST_ZERO ) && ( solution[k] <= 1.0 ) )
                                        {
                                                //float out_collisionTime = solution[k];
  
@@@ -961,27 -936,24 +1070,29 @@@ void cloth_collision_moving ( ClothModi
        cloth_collision_moving_tris ( clmd, coll_clmd, tree1, tree2 );
        cloth_collision_moving_tris ( coll_clmd, clmd, tree2, tree1 );
  }
 +*/
  
 -void cloth_free_collision_list ( ClothModifierData *clmd )
 +int cloth_collision_moving ( ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair, CollPair *collision_end )
  {
 -      // free collision list
 -      if ( clmd->coll_parms->collision_list )
 -      {
 -              LinkNode *search = clmd->coll_parms->collision_list;
 -              while ( search )
 -              {
 -                      CollPair *coll_pair = search->link;
 +      int result = 0;
 +      Cloth *cloth1;
 +      float w1, w2, w3, u1, u2, u3;
 +      float v1[3], v2[3], relativeVelocity[3];
 +      float magrelVel;
 +      float epsilon2 = BLI_bvhtree_getepsilon ( collmd->bvhtree );
  
 -                      MEM_freeN ( coll_pair );
 -                      search = search->next;
 -              }
 -              BLI_linklist_free ( clmd->coll_parms->collision_list,NULL );
 +      cloth1 = clmd->clothObject;
  
 -              clmd->coll_parms->collision_list = NULL;
 +      for ( ; collpair != collision_end; collpair++ )
 +      {
 +              // only handle moving collisions here
 +              if (!( collpair->flag & COLLISION_IN_FUTURE ))
 +                      continue;
 +              
 +              cloth_collision_moving_edges ( clmd, collmd, collpair);
        }
++      
++      return 1;
  }
  
  int cloth_bvh_objcollisions_do ( ClothModifierData * clmd, CollisionModifierData *collmd, float step, float dt )
index 596c381b8964cdad9b3510ebcc18563af68ec74e,596c381b8964cdad9b3510ebcc18563af68ec74e..8dd3b1d0324c0bc1c100cee20c0ff34b99914693
@@@ -46,6 -46,6 +46,7 @@@
  #include "DNA_curve_types.h"
  #include "DNA_group_types.h"
  #include "DNA_scene_types.h"
++#include "DNA_sph_types.h"
  #include "DNA_texture_types.h"
  
  #include "BLI_rand.h"
@@@ -74,6 -74,6 +75,7 @@@
  #include "BKE_mesh.h"
  #include "BKE_modifier.h"
  #include "BKE_scene.h"
++#include "BKE_sph.h"
  
  #include "BSE_headerbuttons.h"
  
@@@ -4634,6 -4634,6 +4636,50 @@@ static void particles_fluid_step(Objec
                snprintf(debugStrBuffer,256,"readFsPartData::done - particles:%d, active:%d, file:%d, mask:%d  \n", psys->totpart,activeParts,fileParts,readMask);
                elbeemDebugOut(debugStrBuffer);
        } // fluid sim particles done
++      else
++      {
++              // check for sph modifier
++              SphModifierData *sphmd = (SphModifierData *)modifiers_findByType(ob, eModifierType_Sph);
++              
++              // check for an sph object
++              if(sphmd)
++              {       
++                      // check for existing coordinates
++                      if(sphmd->sim_parms->co)
++                      {
++                              int i = 0;
++                              ParticleSettings *part = psys->part;
++                              ParticleData *pa=0;
++                              float null[3] = {0,0,0};
++                              
++                              if(ob==G.obedit) // off...
++                                      return;
++                              
++                              part->totpart= sphmd->sim_parms->numpart;
++                              part->sta=part->end = 1.0f;
++                              part->lifetime = G.scene->r.efra + 1;
++
++                              /* initialize particles */
++                              realloc_particles(ob, psys, part->totpart);
++                              initialize_all_particles(ob, psys, 0);
++                              
++                              printf("sphmd->sim_parms->numpart: %ld\n", sphmd->sim_parms->numpart);
++                              
++                              for(i = 0, pa=psys->particles; i < sphmd->sim_parms->numpart; i++, pa++)
++                              {
++                                      pa->size = 0.01f;
++                                      VECCOPY(pa->state.co, sphmd->sim_parms->co + i*3);
++                                      VECCOPY(pa->state.vel, null);
++                                      
++                                      pa->state.ave[0] = pa->state.ave[1] = pa->state.ave[2] = 0.0f;
++                                      pa->state.rot[0] = 1.0;
++                                      pa->state.rot[1] = pa->state.rot[2] = pa->state.rot[3] = 0.0;
++
++                                      pa->alive = PARS_ALIVE;
++                              }
++                      }
++              }
++      }
        #endif // DISABLE_ELBEEM
  }
  
index 51f87b26aafd2ef826d11f0af9e64f7f025c944a,0000000000000000000000000000000000000000..3261984da7604f7e97389536e6138a44266d9e54
mode 100644,000000..100644
--- /dev/null
@@@ -1,63 -1,0 +1,63 @@@
-  * Contributor(s): Daniel (Genscher)
 +/**
 + *
 + * ***** BEGIN GPL LICENSE BLOCK *****
 + *
 + * 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
 + * of the License, or (at your option) any later version.
 + *
 + * This program is distributed in the hope that it will be useful,
 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 + * GNU General Public License for more details.
 + *
 + * You should have received a copy of the GNU General Public License
 + * along with this program; if not, write to the Free Software Foundation,
 + * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 + *
 + * The Original Code is Copyright (C) 2006 by NaN Holding BV.
 + * All rights reserved.
 + *
 + * The Original Code is: all of this file.
 + *
++ * Contributor(s): Daniel Genrich, Jose Pinto
 + *
 + * ***** END GPL LICENSE BLOCK *****
 + */
 +
 +
 +#ifndef BLI_KDOPBVH_H
 +#define BLI_KDOPBVH_H
 +
 +#include <float.h>
 +
 +struct BVHTree;
 +typedef struct BVHTree BVHTree;
 +
 +typedef struct BVHTreeOverlap {
 +      int indexA;
 +      int indexB;
 +} BVHTreeOverlap;
 +
 +BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis);
 +void BLI_bvhtree_free(BVHTree *tree);
 +
 +/* construct: first insert points, then call balance */
 +int BLI_bvhtree_insert(BVHTree *tree, int index, float *co, int numpoints);
 +void BLI_bvhtree_balance(BVHTree *tree);
 +
 +/* update: first update points/nodes, then call update_tree to refit the bounding volumes */
 +int BLI_bvhtree_update_node(BVHTree *tree, int index, float *co, float *co_moving, int numpoints);
 +void BLI_bvhtree_update_tree(BVHTree *tree);
 +
 +/* collision/overlap: check two trees if they overlap, alloc's *overlap with length of the int return value */
 +BVHTreeOverlap *BLI_bvhtree_overlap(BVHTree *tree1, BVHTree *tree2, int *result);
 +
 +float BLI_bvhtree_getepsilon(BVHTree *tree);
 +
 +#endif // BLI_KDOPBVH_H
 +
 +
 +
 +
index 8be52854a7bf29e79db38a8284714f622fce9539,0000000000000000000000000000000000000000..6a1abb5d8ad1b6b5d8a896335ad5047dea8677b2
mode 100644,000000..100644
--- /dev/null
@@@ -1,786 -1,0 +1,768 @@@
- * Contributor(s): Daniel Genrich
 +/*  kdop.c      
 +* 
 +*
 +* ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
 +*
 +* 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
 +* of the License, or (at your option) any later version. The Blender
 +* Foundation also sells licenses for use in proprietary software under
 +* the Blender License.  See http://www.blender.org/BL/ for information
 +* about this.
 +*
 +* This program is distributed in the hope that it will be useful,
 +* but WITHOUT ANY WARRANTY; without even the implied warranty of
 +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 +* GNU General Public License for more details.
 +*
 +* You should have received a copy of the GNU General Public License
 +* along with this program; if not, write to the Free Software Foundation,
 +* Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 +*
 +* The Original Code is Copyright (C) Blender Foundation
 +* All rights reserved.
 +*
 +* The Original Code is: all of this file.
 +*
-       int     *orig_index; /* mapping for orig_index to node_index */
++* Contributor(s): Daniel Genrich, Jose Pinto
 +*
 +* ***** END GPL/BL DUAL LICENSE BLOCK *****
 +*/
 +
 +#include "math.h"
 +#include <stdio.h>
 +#include <stdlib.h> 
 +#include <string.h>
 +
 +#include "MEM_guardedalloc.h"
 +
 +#include "BKE_utildefines.h"
 +
 +#include "BLI_kdopbvh.h"
 +
 +#ifdef _OPENMP
 +#include <omp.h>
 +#endif
 +
 +typedef struct BVHNode
 +{
 +      struct BVHNode *children[8]; // max 8 children
 +      struct BVHNode *parent; // needed for bottom - top update
 +      float bv[26]; // Bounding volume of all nodes, max 13 axis
 +      int index; /* face, edge, vertex index */
 +      char totnode; // how many nodes are used, used for speedup
 +      char traversed;  // how many nodes already traversed until this level?
 +      char main_axis;
 +} BVHNode;
 +
 +struct BVHTree
 +{
 +      BVHNode **nodes;
 +      BVHNode *nodearray; /* pre-alloc branch nodes */
-               MEM_freeN(tree->orig_index);
 +      float   epsilon; /* epslion is used for inflation of the k-dop     */
 +      int     totleaf; // leafs
 +      int     totbranch;
 +      char    tree_type; // type of tree (4 => quadtree)
 +      char    axis; // kdop type (6 => OBB, 7 => AABB, ...)
 +      char    start_axis, stop_axis; // KDOP_AXES array indices according to axis
 +};
 +
 +typedef struct BVHOverlapData 
 +{  
 +      BVHTree *tree1, *tree2; 
 +      BVHTreeOverlap *overlap; 
 +      int i, max_overlap; /* i is number of overlaps */
 +} BVHOverlapData;
 +////////////////////////////////////////
 +
 +
 +////////////////////////////////////////////////////////////////////////
 +// Bounding Volume Hierarchy Definition
 +// 
 +// Notes: From OBB until 26-DOP --> all bounding volumes possible, just choose type below
 +// Notes: You have to choose the type at compile time ITM
 +// Notes: You can choose the tree type --> binary, quad, octree, choose below
 +////////////////////////////////////////////////////////////////////////
 +
 +static float KDOP_AXES[13][3] =
 +{ {1.0, 0, 0}, {0, 1.0, 0}, {0, 0, 1.0}, {1.0, 1.0, 1.0}, {1.0, -1.0, 1.0}, {1.0, 1.0, -1.0},
 +{1.0, -1.0, -1.0}, {1.0, 1.0, 0}, {1.0, 0, 1.0}, {0, 1.0, 1.0}, {1.0, -1.0, 0}, {1.0, 0, -1.0},
 +{0, 1.0, -1.0}
 +};
 +
 +//////////////////////////////////////////////////////////////////////////////////////////////////////
 +// Introsort 
 +// with permission deriven from the following Java code:
 +// http://ralphunden.net/content/tutorials/a-guide-to-introsort/
 +// and he derived it from the SUN STL 
 +//////////////////////////////////////////////////////////////////////////////////////////////////////
 +static int size_threshold = 16;
 +/*
 +* Common methods for all algorithms
 +*/
 +static void bvh_exchange(BVHNode **a, int i, int j)
 +{
 +      BVHNode *t=a[i];
 +      a[i]=a[j];
 +      a[j]=t;
 +}
 +static int floor_lg(int a)
 +{
 +      return (int)(floor(log(a)/log(2)));
 +}
 +
 +/*
 +* Insertion sort algorithm
 +*/
 +static void bvh_insertionsort(BVHNode **a, int lo, int hi, int axis)
 +{
 +      int i,j;
 +      BVHNode *t;
 +      for (i=lo; i < hi; i++)
 +      {
 +              j=i;
 +              t = a[i];
 +              while((j!=lo) && (t->bv[axis] < (a[j-1])->bv[axis]))
 +              {
 +                      a[j] = a[j-1];
 +                      j--;
 +              }
 +              a[j] = t;
 +      }
 +}
 +
 +static int bvh_partition(BVHNode **a, int lo, int hi, BVHNode * x, int axis)
 +{
 +      int i=lo, j=hi;
 +      while (1)
 +      {
 +              while ((a[i])->bv[axis] < x->bv[axis]) i++;
 +              j=j-1;
 +              while (x->bv[axis] < (a[j])->bv[axis]) j=j-1;
 +              if(!(i < j))
 +                      return i;
 +              bvh_exchange(a, i,j);
 +              i++;
 +      }
 +}
 +
 +/*
 +* Heapsort algorithm
 +*/
 +static void bvh_downheap(BVHNode **a, int i, int n, int lo, int axis)
 +{
 +      BVHNode * d = a[lo+i-1];
 +      int child;
 +      while (i<=n/2)
 +      {
 +              child = 2*i;
 +              if ((child < n) && ((a[lo+child-1])->bv[axis] < (a[lo+child])->bv[axis]))
 +              {
 +                      child++;
 +              }
 +              if (!(d->bv[axis] < (a[lo+child-1])->bv[axis])) break;
 +              a[lo+i-1] = a[lo+child-1];
 +              i = child;
 +      }
 +      a[lo+i-1] = d;
 +}
 +
 +static void bvh_heapsort(BVHNode **a, int lo, int hi, int axis)
 +{
 +      int n = hi-lo, i;
 +      for (i=n/2; i>=1; i=i-1)
 +      {
 +              bvh_downheap(a, i,n,lo, axis);
 +      }
 +      for (i=n; i>1; i=i-1)
 +      {
 +              bvh_exchange(a, lo,lo+i-1);
 +              bvh_downheap(a, 1,i-1,lo, axis);
 +      }
 +}
 +
 +static BVHNode *bvh_medianof3(BVHNode **a, int lo, int mid, int hi, int axis) // returns Sortable
 +{
 +      if ((a[mid])->bv[axis] < (a[lo])->bv[axis])
 +      {
 +              if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
 +                      return a[mid];
 +              else
 +              {
 +                      if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
 +                              return a[hi];
 +                      else
 +                              return a[lo];
 +              }
 +      }
 +      else
 +      {
 +              if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
 +              {
 +                      if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
 +                              return a[lo];
 +                      else
 +                              return a[hi];
 +              }
 +              else
 +                      return a[mid];
 +      }
 +}
 +/*
 +* Quicksort algorithm modified for Introsort
 +*/
 +static void bvh_introsort_loop (BVHNode **a, int lo, int hi, int depth_limit, int axis)
 +{
 +      int p;
 +
 +      while (hi-lo > size_threshold)
 +      {
 +              if (depth_limit == 0)
 +              {
 +                      bvh_heapsort(a, lo, hi, axis);
 +                      return;
 +              }
 +              depth_limit=depth_limit-1;
 +              p=bvh_partition(a, lo, hi, bvh_medianof3(a, lo, lo+((hi-lo)/2)+1, hi-1, axis), axis);
 +              bvh_introsort_loop(a, p, hi, depth_limit, axis);
 +              hi=p;
 +      }
 +}
 +
 +static void sort(BVHNode **a0, int begin, int end, int axis)
 +{
 +      if (begin < end)
 +      {
 +              BVHNode **a=a0;
 +              bvh_introsort_loop(a, begin, end, 2*floor_lg(end-begin), axis);
 +              bvh_insertionsort(a, begin, end, axis);
 +      }
 +}
 +void sort_along_axis(BVHTree *tree, int start, int end, int axis)
 +{
 +      sort(tree->nodes, start, end, axis);
 +}
 +
 +//////////////////////////////////////////////////////////////////////////////////////////////////////
 +
 +void BLI_bvhtree_free(BVHTree *tree)
 +{     
 +      if(tree)
 +      {
 +              MEM_freeN(tree->nodes);
 +              MEM_freeN(tree->nodearray);
-               tree->orig_index = (int *)MEM_callocN(sizeof(int)*(numbranches+maxsize + tree_type), "BVHIndexArray");
-               
-               if(!tree->orig_index)
-               {
-                       MEM_freeN(tree);
-                       MEM_freeN(tree->nodes);
-                       MEM_freeN(tree->nodearray);
-                       return NULL;
-               }
-               
 +              MEM_freeN(tree);
 +      }
 +}
 +
 +BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis)
 +{
 +      BVHTree *tree;
 +      int numbranches=0, i;
 +      
 +      // only support up to octree
 +      if(tree_type > 8)
 +              return NULL;
 +
 +      tree = (BVHTree *)MEM_callocN(sizeof(BVHTree), "BVHTree");
 +      
 +      if(tree)
 +      {
 +              // calculate max number of branches, our bvh kdop is "almost perfect"
 +              for(i = 1; i <= (int)ceil((float)((float)log(maxsize)/(float)log(tree_type))); i++)
 +                      numbranches += (pow(tree_type, i) / tree_type);
 +              
 +              tree->nodes = (BVHNode **)MEM_callocN(sizeof(BVHNode *)*(numbranches+maxsize + tree_type), "BVHNodes");
 +              
 +              if(!tree->nodes)
 +              {
 +                      MEM_freeN(tree);
 +                      return NULL;
 +              }
 +              
 +              tree->nodearray = (BVHNode *)MEM_callocN(sizeof(BVHNode)*(numbranches+maxsize + tree_type), "BVHNodeArray");
 +              
 +              if(!tree->nodearray)
 +              {
 +                      MEM_freeN(tree);
 +                      MEM_freeN(tree->nodes);
 +                      return NULL;
 +              }
 +              
-       // put indices into array for O(1) access
-       for(i = 0; i < tree->totleaf; i++)
-       {
-               tree->orig_index[tree->nodes[i]->index] = i;
-       }
-       
 +              tree->epsilon = epsilon;
 +              tree->tree_type = tree_type; 
 +              tree->axis = axis;
 +              
 +              if(axis == 26)
 +              {
 +                      tree->start_axis = 0;
 +                      tree->stop_axis = 13;
 +              }
 +              else if(axis == 18)
 +              {
 +                      tree->start_axis = 7;
 +                      tree->stop_axis = 13;
 +              }
 +              else if(axis == 14)
 +              {
 +                      tree->start_axis = 0;
 +                      tree->stop_axis = 7;
 +              }
 +              else if(axis == 8) // AABB
 +              {
 +                      tree->start_axis = 0;
 +                      tree->stop_axis = 4;
 +              }
 +              else if(axis == 6) // OBB
 +              {
 +                      tree->start_axis = 0;
 +                      tree->stop_axis = 3;
 +              }
 +              else
 +              {
 +                      BLI_bvhtree_free(tree);
 +                      return NULL;
 +              }
 +      }
 +
 +      return tree;
 +}
 +
 +static void create_kdop_hull(BVHTree *tree, BVHNode *node, float *co, int numpoints, int moving)
 +{
 +      float newminmax;
 +      int i, k;
 +      
 +      // don't init boudings for the moving case
 +      if(!moving)
 +      {
 +              for (i = tree->start_axis; i < tree->stop_axis; i++)
 +              {
 +                      node->bv[2*i] = FLT_MAX;
 +                      node->bv[2*i + 1] = -FLT_MAX;
 +              }
 +      }
 +      
 +      for(k = 0; k < numpoints; k++)
 +      {
 +              // for all Axes.
 +              for (i = tree->start_axis; i < tree->stop_axis; i++)
 +              {
 +                      newminmax = INPR(&co[k * 3], KDOP_AXES[i]);
 +                      if (newminmax < node->bv[2 * i])
 +                              node->bv[2 * i] = newminmax;
 +                      if (newminmax > node->bv[(2 * i) + 1])
 +                              node->bv[(2 * i) + 1] = newminmax;
 +              }
 +      }
 +}
 +
 +// depends on the fact that the BVH's for each face is already build
 +static void refit_kdop_hull(BVHTree *tree, int start, int end, float *bv)
 +{
 +      float newmin,newmax;
 +      int i, j;
 +      
 +      for (i = tree->start_axis; i < tree->stop_axis; i++)
 +      {
 +              bv[2*i] = FLT_MAX;
 +              bv[2*i + 1] = -FLT_MAX;
 +      }
 +
 +      for (j = start; j < end; j++)
 +      {
 +                // for all Axes.
 +              for (i = tree->start_axis; i < tree->stop_axis; i++)
 +              {
 +                      newmin = tree->nodes[j]->bv[(2 * i)];   
 +                      if ((newmin < bv[(2 * i)]))
 +                              bv[(2 * i)] = newmin;
 + 
 +                      newmax = tree->nodes[j]->bv[(2 * i) + 1];
 +                      if ((newmax > bv[(2 * i) + 1]))
 +                              bv[(2 * i) + 1] = newmax;
 +              }
 +      }
 +}
 +
 +int BLI_bvhtree_insert(BVHTree *tree, int index, float *co, int numpoints)
 +{
 +      BVHNode *node= NULL;
 +      int i;
 +      
 +      // insert should only possible as long as tree->totbranch is 0
 +      if(tree->totbranch > 0)
 +              return 0;
 +      
 +      if(tree->totleaf+1 >= MEM_allocN_len(tree->nodes))
 +              return 0;
 +      
 +      // TODO check if have enough nodes in array
 +      
 +      node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
 +      tree->totleaf++;
 +      
 +      create_kdop_hull(tree, node, co, numpoints, 0);
 +      
 +      // inflate the bv with some epsilon
 +      for (i = tree->start_axis; i < tree->stop_axis; i++)
 +      {
 +              node->bv[(2 * i)] -= tree->epsilon; // minimum 
 +              node->bv[(2 * i) + 1] += tree->epsilon; // maximum 
 +      }
 +
 +      node->index= index;
 +      
 +      return 1;
 +}
 +
 +// only supports x,y,z axis in the moment
 +// but we should use a plain and simple function here for speed sake
 +static char get_largest_axis(float *bv)
 +{
 +      float middle_point[3];
 +
 +      middle_point[0] = (bv[1]) - (bv[0]); // x axis
 +      middle_point[1] = (bv[3]) - (bv[2]); // y axis
 +      middle_point[2] = (bv[5]) - (bv[4]); // z axis
 +      if (middle_point[0] > middle_point[1]) 
 +      {
 +              if (middle_point[0] > middle_point[2])
 +                      return 1; // max x axis
 +              else
 +                      return 5; // max z axis
 +      }
 +      else 
 +      {
 +              if (middle_point[1] > middle_point[2])
 +                      return 3; // max y axis
 +              else
 +                      return 5; // max z axis
 +      }
 +}
 +
 +static void bvh_div_nodes(BVHTree *tree, BVHNode *node, int start, int end, char lastaxis)
 +{
 +      char laxis;
 +      int i, tend;
 +      BVHNode *tnode;
 +      int slice = (end-start+tree->tree_type-1)/tree->tree_type;      //division rounded up
 +      
 +      // Determine which axis to split along
 +      laxis = get_largest_axis(node->bv);
 +
 +      // Sort along longest axis
 +      if(laxis!=lastaxis)
 +              sort_along_axis(tree, start, end, laxis);
 +      
 +      // split nodes along longest axis
 +      for (i=0; start < end; start += slice, i++) //i counts the current child
 +      {       
 +              tend = start + slice;
 +              
 +              if(tend > end) tend = end;
 +              
 +              if(tend-start == 1)     // ok, we have 1 left for this node
 +              {
 +                      node->children[i] = tree->nodes[start];
 +                      node->children[i]->parent = node;
 +              }
 +              else
 +              {
 +                      tnode = node->children[i] = tree->nodes[tree->totleaf  + tree->totbranch] = &(tree->nodearray[tree->totbranch + tree->totleaf]);
 +                      tree->totbranch++;
 +                      tnode->parent = node;
 +
 +                      refit_kdop_hull(tree, start, tend, tnode->bv);
 +                      bvh_div_nodes(tree, tnode, start, tend, laxis);
 +              }
 +              node->totnode++;
 +      }
 +      
 +      return;
 +}
 +
 +static void verify_tree(BVHTree *tree)
 +{
 +      int i, j, check = 0;
 +      
 +      // check the pointer list
 +      for(i = 0; i < tree->totleaf; i++)
 +      {
 +              if(tree->nodes[i]->parent == NULL)
 +                      printf("Leaf has no parent: %d\n", i);
 +              else
 +              {
 +                      for(j = 0; j < tree->tree_type; j++)
 +                      {
 +                              if(tree->nodes[i]->parent->children[j] == tree->nodes[i])
 +                                      check = 1;
 +                      }
 +                      if(!check)
 +                      {
 +                              printf("Parent child relationship doesn't match: %d\n", i);
 +                      }
 +                      check = 0;
 +              }
 +      }
 +      
 +      // check the leaf list
 +      for(i = 0; i < tree->totleaf; i++)
 +      {
 +              if(tree->nodearray[i].parent == NULL)
 +                      printf("Leaf has no parent: %d\n", i);
 +              else
 +              {
 +                      for(j = 0; j < tree->tree_type; j++)
 +                      {
 +                              if(tree->nodearray[i].parent->children[j] == &tree->nodearray[i])
 +                                      check = 1;
 +                      }
 +                      if(!check)
 +                      {
 +                              printf("Parent child relationship doesn't match: %d\n", i);
 +                      }
 +                      check = 0;
 +              }
 +      }
 +      
 +      printf("branches: %d, leafs: %d, total: %d\n", tree->totbranch, tree->totleaf, tree->totbranch + tree->totleaf);
 +}
 +      
 +void BLI_bvhtree_balance(BVHTree *tree)
 +{
 +      BVHNode *node;
 +      int i;
 +      
 +      if(tree->totleaf == 0)
 +              return;
 +      
 +      // create root node
 +      node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
 +      tree->totbranch++;
 +      
 +      // refit root bvh node
 +      refit_kdop_hull(tree, 0, tree->totleaf, tree->nodes[tree->totleaf]->bv);
 +      // create + balance tree
 +      bvh_div_nodes(tree, tree->nodes[tree->totleaf], 0, tree->totleaf, 0);
 +      
-       node = tree->nodes[tree->orig_index[index]];
 +      verify_tree(tree);
 +}
 +
 +// overlap - is it possbile for 2 bv's to collide ?
 +static int tree_overlap(float *bv1, float *bv2, int start_axis, int stop_axis)
 +{
 +      float *bv1_end = bv1 + (stop_axis<<1);
 +              
 +      bv1 += start_axis<<1;
 +      bv2 += start_axis<<1;
 +      
 +      // test all axis if min + max overlap
 +      for (; bv1 != bv1_end; bv1+=2, bv2+=2)
 +      {
 +              if ((*(bv1) > *(bv2 + 1)) || (*(bv2) > *(bv1 + 1))) 
 +                      return 0;
 +      }
 +      
 +      return 1;
 +}
 +
 +static void traverse(BVHOverlapData *data, BVHNode *node1, BVHNode *node2)
 +{
 +      int j;
 +      
 +      if(tree_overlap(node1->bv, node2->bv, MIN2(data->tree1->start_axis, data->tree2->start_axis), MIN2(data->tree1->stop_axis, data->tree2->stop_axis)))
 +      {
 +              // check if node1 is a leaf
 +              if(node1->index)
 +              {
 +                      // check if node2 is a leaf
 +                      if(node2->index)
 +                      {
 +                              
 +                              if(node1 == node2)
 +                              {
 +                                      return;
 +                              }
 +                                      
 +                              if(data->i >= data->max_overlap)
 +                              {       
 +                                      // try to make alloc'ed memory bigger
 +                                      data->overlap = realloc(data->overlap, sizeof(BVHTreeOverlap)*data->max_overlap*2);
 +                                      
 +                                      if(!data->overlap)
 +                                      {
 +                                              printf("Out of Memory in traverse\n");
 +                                              return;
 +                                      }
 +                                      data->max_overlap *= 2;
 +                              }
 +                              
 +                              // both leafs, insert overlap!
 +                              data->overlap[data->i].indexA = node1->index;
 +                              data->overlap[data->i].indexB = node2->index;
 +
 +                              data->i++;
 +                      }
 +                      else
 +                      {
 +                              for(j = 0; j < data->tree2->tree_type; j++)
 +                              {
 +                                      if(node2->children[j])
 +                                              traverse(data, node1, node2->children[j]);
 +                              }
 +                      }
 +              }
 +              else
 +              {
 +                      
 +                      for(j = 0; j < data->tree2->tree_type; j++)
 +                      {
 +                              if(node1->children[j])
 +                                      traverse(data, node1->children[j], node2);
 +                      }
 +              }
 +      }
 +      return;
 +}
 +
 +BVHTreeOverlap *BLI_bvhtree_overlap(BVHTree *tree1, BVHTree *tree2, int *result)
 +{
 +      int j, total = 0;
 +      BVHTreeOverlap *overlap = NULL, *to = NULL;
 +      BVHOverlapData *data[tree1->tree_type];
 +      
 +      // check for compatibility of both trees (can't compare 14-DOP with 18-DOP)
 +      if((tree1->axis != tree2->axis) && ((tree1->axis == 14) || tree2->axis == 14))
 +              return 0;
 +      
 +      // fast check root nodes for collision before doing big splitting + traversal
 +      if(!tree_overlap(tree1->nodes[tree1->totleaf]->bv, tree2->nodes[tree2->totleaf]->bv, MIN2(tree1->start_axis, tree2->start_axis), MIN2(tree1->stop_axis, tree2->stop_axis)))
 +              return 0;
 +      
 +      for(j = 0; j < tree1->tree_type; j++)
 +      {
 +              data[j] = (BVHOverlapData *)MEM_callocN(sizeof(BVHOverlapData), "BVHOverlapData");
 +              
 +              // init BVHOverlapData
 +              data[j]->overlap = (BVHTreeOverlap *)malloc(sizeof(BVHTreeOverlap)*MAX2(tree1->totleaf, tree2->totleaf));
 +              data[j]->tree1 = tree1;
 +              data[j]->tree2 = tree2;
 +              data[j]->max_overlap = MAX2(tree1->totleaf, tree2->totleaf);
 +              data[j]->i = 0;
 +      }
 +
 +#pragma omp parallel for private(j) schedule(static)
 +      for(j = 0; j < tree1->tree_type; j++)
 +      {
 +              traverse(data[j], tree1->nodes[tree1->totleaf]->children[j], tree2->nodes[tree2->totleaf]);
 +      }
 +      
 +      for(j = 0; j < tree1->tree_type; j++)
 +              total += data[j]->i;
 +      
 +      to = overlap = (BVHTreeOverlap *)MEM_callocN(sizeof(BVHTreeOverlap)*total, "BVHTreeOverlap");
 +      
 +      for(j = 0; j < tree1->tree_type; j++)
 +      {
 +              memcpy(to, data[j]->overlap, data[j]->i*sizeof(BVHTreeOverlap));
 +              to+=data[j]->i;
 +      }
 +      
 +      for(j = 0; j < tree1->tree_type; j++)
 +      {
 +              free(data[j]->overlap);
 +              MEM_freeN(data[j]);
 +      }
 +      
 +      (*result) = total;
 +      return overlap;
 +}
 +
 +
 +// bottom up update of bvh tree:
 +// join the 4 children here
 +static void node_join(BVHTree *tree, BVHNode *node)
 +{
 +      int i, j;
 +      
 +      for (i = tree->start_axis; i < tree->stop_axis; i++)
 +      {
 +              node->bv[2*i] = FLT_MAX;
 +              node->bv[2*i + 1] = -FLT_MAX;
 +      }
 +      
 +      for (i = 0; i < tree->tree_type; i++)
 +      {
 +              if (node->children[i]) 
 +              {
 +                      for (j = tree->start_axis; j < tree->stop_axis; j++)
 +                      {
 +                              // update minimum 
 +                              if (node->children[i]->bv[(2 * j)] < node->bv[(2 * j)]) 
 +                                      node->bv[(2 * j)] = node->children[i]->bv[(2 * j)];
 +                              
 +                              // update maximum 
 +                              if (node->children[i]->bv[(2 * j) + 1] > node->bv[(2 * j) + 1])
 +                                      node->bv[(2 * j) + 1] = node->children[i]->bv[(2 * j) + 1];
 +                      }
 +              }
 +              else
 +                      break;
 +      }
 +}
 +
 +// call before BLI_bvhtree_update_tree()
 +int BLI_bvhtree_update_node(BVHTree *tree, int index, float *co, float *co_moving, int numpoints)
 +{
 +      BVHNode *node= NULL;
 +      int i = 0;
 +      
 +      // check if index exists
 +      if(index > tree->totleaf)
 +              return 0;
 +      
++      node = tree->nodearray + index;
 +      
 +      create_kdop_hull(tree, node, co, numpoints, 0);
 +      
 +      if(co_moving)
 +              create_kdop_hull(tree, node, co_moving, numpoints, 1);
 +      
 +      // inflate the bv with some epsilon
 +      for (i = tree->start_axis; i < tree->stop_axis; i++)
 +      {
 +              node->bv[(2 * i)] -= tree->epsilon; // minimum 
 +              node->bv[(2 * i) + 1] += tree->epsilon; // maximum 
 +      }
 +      
 +      return 1;
 +}
 +
 +// call BLI_bvhtree_update_node() first for every node/point/triangle
 +void BLI_bvhtree_update_tree(BVHTree *tree)
 +{
 +      BVHNode *leaf, *parent;
 +      
 +      // reset tree traversing flag
 +      for (leaf = tree->nodearray + tree->totleaf; leaf != tree->nodearray + tree->totleaf + tree->totbranch; leaf++)
 +              leaf->traversed = 0;
 +      
 +      for (leaf = tree->nodearray; leaf != tree->nodearray + tree->totleaf; leaf++)
 +      {
 +              for (parent = leaf->parent; parent; parent = parent->parent)
 +              {
 +                      parent->traversed++;    // we tried to go up in hierarchy 
 +                      if (parent->traversed < parent->totnode) 
 +                              break;  // we do not need to check further 
 +                      else 
 +                              node_join(tree, parent);
 +              }
 +      }
 +}
 +
 +float BLI_bvhtree_getepsilon(BVHTree *tree)
 +{
 +      return tree->epsilon;
 +}
Simple merge
Simple merge
Simple merge