Extras/ConvexDecomposition/concavity.cpp

   1 #include "float_math.h"
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <string.h>
   5 #include <assert.h>
   6
   7 #include <vector>
   8
   9 /*----------------------------------------------------------------------
  10                 Copyright (c) 2004 Open Dynamics Framework Group
  11                                         www.physicstools.org
  12                 All rights reserved.
  13
  14                 Redistribution and use in source and binary forms, with or without modification, are permitted provided
  15                 that the following conditions are met:
  16
  17                 Redistributions of source code must retain the above copyright notice, this list of conditions
  18                 and the following disclaimer.
  19
  20                 Redistributions in binary form must reproduce the above copyright notice,
  21                 this list of conditions and the following disclaimer in the documentation
  22                 and/or other materials provided with the distribution.
  23
  24                 Neither the name of the Open Dynamics Framework Group nor the names of its contributors may
  25                 be used to endorse or promote products derived from this software without specific prior written permission.
  26
  27                 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 'AS IS' AND ANY EXPRESS OR IMPLIED WARRANTIES,
  28                 INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  29                 DISCLAIMED. IN NO EVENT SHALL THE INTEL OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  30                 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  31                 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
  32                 IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  33                 THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  34 -----------------------------------------------------------------------*/
  35
  36 // http://codesuppository.blogspot.com
  37 //
  38 // mailto: jratcliff@infiniplex.net
  39 //
  40 // http://www.amillionpixels.us
  41 //
  42
  43 #include "concavity.h"
  44 #include "raytri.h"
  45 #include "bestfit.h"
  46 #include "cd_hull.h"
  47 #include "meshvolume.h"
  48 #include "cd_vector.h"
  49 #include "splitplane.h"
  50 #include "ConvexDecomposition.h"
  51
  52
  53 #define WSCALE 4
  54 #define CONCAVE_THRESH 0.05f
  55
  56 namespace ConvexDecomposition
  57 {
  58
  59 unsigned int getDebugColor(void)
  60 {
  61         static unsigned int colors[8] =
  62         {
  63                 0xFF0000,
  64           0x00FF00,
  65                 0x0000FF,
  66                 0xFFFF00,
  67                 0x00FFFF,
  68                 0xFF00FF,
  69                 0xFFFFFF,
  70                 0xFF8040
  71         };
  72
  73         static int count = 0;
  74
  75         count++;
  76
  77         if ( count == 8 ) count = 0;
  78
  79         assert( count >= 0 && count < 8 );
  80
  81         unsigned int color = colors[count];
  82
  83   return color;
  84
  85 }
  86
  87 class Wpoint
  88 {
  89 public:
  90   Wpoint(const Vector3d &p,float w)
  91   {
  92     mPoint = p;
  93     mWeight = w;
  94   }
  95
  96   Vector3d mPoint;
  97   float           mWeight;
  98 };
  99
 100 typedef std::vector< Wpoint > WpointVector;
 101
 102
 103 static inline float DistToPt(const float *p,const float *plane)
 104 {
 105         float x = p[0];
 106         float y = p[1];
 107         float z = p[2];
 108         float d = x*plane[0] + y*plane[1] + z*plane[2] + plane[3];
 109         return d;
 110 }
 111
 112
 113 static void intersect(const float *p1,const float *p2,float *split,const float *plane)
 114 {
 115
 116   float dp1 = DistToPt(p1,plane);
 117
 118   float dir[3];
 119
 120   dir[0] = p2[0] - p1[0];
 121   dir[1] = p2[1] - p1[1];
 122   dir[2] = p2[2] - p1[2];
 123
 124   float dot1 = dir[0]*plane[0] + dir[1]*plane[1] + dir[2]*plane[2];
 125   float dot2 = dp1 - plane[3];
 126
 127   float    t = -(plane[3] + dot2 ) / dot1;
 128
 129   split[0] = (dir[0]*t)+p1[0];
 130   split[1] = (dir[1]*t)+p1[1];
 131   split[2] = (dir[2]*t)+p1[2];
 132 }
 133
 134
 135 class CTri
 136 {
 137 public:
 138         CTri(void) { };
 139
 140   CTri(const float *p1,const float *p2,const float *p3,unsigned int i1,unsigned int i2,unsigned int i3)
 141   {
 142     mProcessed = 0;
 143     mI1 = i1;
 144     mI2 = i2;
 145     mI3 = i3;
 146
 147         mP1.Set(p1);
 148         mP2.Set(p2);
 149         mP3.Set(p3);
 150
 151         mPlaneD = mNormal.ComputePlane(mP1,mP2,mP3);
 152         }
 153
 154   float Facing(const CTri &t)
 155   {
 156                 float d = mNormal.Dot(t.mNormal);
 157                 return d;
 158   }
 159
 160   // clip this line segment against this triangle.
 161   bool clip(const Vector3d &start,Vector3d &end) const
 162   {
 163     Vector3d sect;
 164
 165     bool hit = lineIntersectsTriangle(start.Ptr(), end.Ptr(), mP1.Ptr(), mP2.Ptr(), mP3.Ptr(), sect.Ptr() );
 166
 167     if ( hit )
 168     {
 169       end = sect;
 170     }
 171     return hit;
 172   }
 173
 174         bool Concave(const Vector3d &p,float &distance,Vector3d &n) const
 175         {
 176                 n.NearestPointInTriangle(p,mP1,mP2,mP3);
 177                 distance = p.Distance(n);
 178                 return true;
 179         }
 180
 181         void addTri(unsigned int *indices,unsigned int i1,unsigned int i2,unsigned int i3,unsigned int &tcount) const
 182         {
 183                 indices[tcount*3+0] = i1;
 184                 indices[tcount*3+1] = i2;
 185                 indices[tcount*3+2] = i3;
 186                 tcount++;
 187         }
 188
 189         float getVolume(ConvexDecompInterface *callback) const
 190         {
 191                 unsigned int indices[8*3];
 192
 193
 194     unsigned int tcount = 0;
 195
 196     addTri(indices,0,1,2,tcount);
 197     addTri(indices,3,4,5,tcount);
 198
 199     addTri(indices,0,3,4,tcount);
 200     addTri(indices,0,4,1,tcount);
 201
 202     addTri(indices,1,4,5,tcount);
 203     addTri(indices,1,5,2,tcount);
 204
 205     addTri(indices,0,3,5,tcount);
 206     addTri(indices,0,5,2,tcount);
 207
 208     const float *vertices = mP1.Ptr();
 209
 210                 if ( callback )
 211                 {
 212                         unsigned int color = getDebugColor();
 213
 214 #if 0
 215                 Vector3d d1 = mNear1;
 216                 Vector3d d2 = mNear2;
 217                 Vector3d d3 = mNear3;
 218
 219                 callback->ConvexDebugPoint(mP1.Ptr(),0.01f,0x00FF00);
 220                 callback->ConvexDebugPoint(mP2.Ptr(),0.01f,0x00FF00);
 221                 callback->ConvexDebugPoint(mP3.Ptr(),0.01f,0x00FF00);
 222                 callback->ConvexDebugPoint(d1.Ptr(),0.01f,0xFF0000);
 223                 callback->ConvexDebugPoint(d2.Ptr(),0.01f,0xFF0000);
 224                 callback->ConvexDebugPoint(d3.Ptr(),0.01f,0xFF0000);
 225
 226                 callback->ConvexDebugTri(mP1.Ptr(), d1.Ptr(),  d1.Ptr(),0x00FF00);
 227                 callback->ConvexDebugTri(mP2.Ptr(), d2.Ptr(),  d2.Ptr(),0x00FF00);
 228                 callback->ConvexDebugTri(mP3.Ptr(), d3.Ptr(),  d3.Ptr(),0x00FF00);
 229
 230 #else
 231                         for (unsigned int i=0; i<tcount; i++)
 232                         {
 233                                 unsigned int i1 = indices[i*3+0];
 234                                 unsigned int i2 = indices[i*3+1];
 235                                 unsigned int i3 = indices[i*3+2];
 236
 237                                 const float *p1 = &vertices[ i1*3 ];
 238                                 const float *p2 = &vertices[ i2*3 ];
 239                                 const float *p3 = &vertices[ i3*3 ];
 240
 241                                 callback->ConvexDebugTri(p1,p2,p3,color);
 242
 243                         }
 244 #endif
 245                 }
 246
 247                 float v = computeMeshVolume(mP1.Ptr(), tcount, indices );
 248
 249                 return v;
 250
 251         }
 252
 253         float raySect(const Vector3d &p,const Vector3d &dir,Vector3d &sect) const
 254         {
 255                 float plane[4];
 256
 257     plane[0] = mNormal.x;
 258     plane[1] = mNormal.y;
 259     plane[2] = mNormal.z;
 260     plane[3] = mPlaneD;
 261
 262                 Vector3d dest = p+dir*100000;
 263
 264     intersect( p.Ptr(), dest.Ptr(), sect.Ptr(), plane );
 265
 266     return sect.Distance(p); // return the intersection distance.
 267
 268         }
 269
 270   float planeDistance(const Vector3d &p) const
 271   {
 272                 float plane[4];
 273
 274     plane[0] = mNormal.x;
 275     plane[1] = mNormal.y;
 276     plane[2] = mNormal.z;
 277     plane[3] = mPlaneD;
 278
 279                 return DistToPt( p.Ptr(), plane );
 280
 281   }
 282
 283         bool samePlane(const CTri &t) const
 284         {
 285                 const float THRESH = 0.001f;
 286     float dd = fabsf( t.mPlaneD - mPlaneD );
 287     if ( dd > THRESH ) return false;
 288     dd = fabsf( t.mNormal.x - mNormal.x );
 289     if ( dd > THRESH ) return false;
 290     dd = fabsf( t.mNormal.y - mNormal.y );
 291     if ( dd > THRESH ) return false;
 292     dd = fabsf( t.mNormal.z - mNormal.z );
 293     if ( dd > THRESH ) return false;
 294     return true;
 295         }
 296
 297         bool hasIndex(unsigned int i) const
 298         {
 299                 if ( i == mI1 || i == mI2 || i == mI3 ) return true;
 300                 return false;
 301         }
 302
 303   bool sharesEdge(const CTri &t) const
 304   {
 305     bool ret = false;
 306     unsigned int count = 0;
 307
 308                 if ( t.hasIndex(mI1) ) count++;
 309           if ( t.hasIndex(mI2) ) count++;
 310                 if ( t.hasIndex(mI3) ) count++;
 311
 312     if ( count >= 2 ) ret = true;
 313
 314     return ret;
 315   }
 316
 317   void debug(unsigned int color,ConvexDecompInterface *callback)
 318   {
 319     callback->ConvexDebugTri( mP1.Ptr(), mP2.Ptr(), mP3.Ptr(), color );
 320     callback->ConvexDebugTri( mP1.Ptr(), mP1.Ptr(), mNear1.Ptr(), 0xFF0000 );
 321     callback->ConvexDebugTri( mP2.Ptr(), mP2.Ptr(), mNear2.Ptr(), 0xFF0000 );
 322     callback->ConvexDebugTri( mP2.Ptr(), mP3.Ptr(), mNear3.Ptr(), 0xFF0000 );
 323     callback->ConvexDebugPoint( mNear1.Ptr(), 0.01f, 0xFF0000 );
 324     callback->ConvexDebugPoint( mNear2.Ptr(), 0.01f, 0xFF0000 );
 325     callback->ConvexDebugPoint( mNear3.Ptr(), 0.01f, 0xFF0000 );
 326   }
 327
 328   float area(void)
 329   {
 330                 float a = mConcavity*mP1.Area(mP2,mP3);
 331     return a;
 332   }
 333
 334   void addWeighted(WpointVector &list,ConvexDecompInterface *callback)
 335   {
 336
 337     Wpoint p1(mP1,mC1);
 338     Wpoint p2(mP2,mC2);
 339     Wpoint p3(mP3,mC3);
 340
 341                 Vector3d d1 = mNear1 - mP1;
 342                 Vector3d d2 = mNear2 - mP2;
 343                 Vector3d d3 = mNear3 - mP3;
 344
 345                 d1*=WSCALE;
 346                 d2*=WSCALE;
 347                 d3*=WSCALE;
 348
 349                 d1 = d1 + mP1;
 350                 d2 = d2 + mP2;
 351           d3 = d3 + mP3;
 352
 353     Wpoint p4(d1,mC1);
 354     Wpoint p5(d2,mC2);
 355     Wpoint p6(d3,mC3);
 356
 357     list.push_back(p1);
 358     list.push_back(p2);
 359     list.push_back(p3);
 360
 361     list.push_back(p4);
 362     list.push_back(p5);
 363     list.push_back(p6);
 364
 365 #if 0
 366                 callback->ConvexDebugPoint(mP1.Ptr(),0.01f,0x00FF00);
 367                 callback->ConvexDebugPoint(mP2.Ptr(),0.01f,0x00FF00);
 368                 callback->ConvexDebugPoint(mP3.Ptr(),0.01f,0x00FF00);
 369                 callback->ConvexDebugPoint(d1.Ptr(),0.01f,0xFF0000);
 370                 callback->ConvexDebugPoint(d2.Ptr(),0.01f,0xFF0000);
 371                 callback->ConvexDebugPoint(d3.Ptr(),0.01f,0xFF0000);
 372
 373                 callback->ConvexDebugTri(mP1.Ptr(), d1.Ptr(),  d1.Ptr(),0x00FF00);
 374                 callback->ConvexDebugTri(mP2.Ptr(), d2.Ptr(),  d2.Ptr(),0x00FF00);
 375                 callback->ConvexDebugTri(mP3.Ptr(), d3.Ptr(),  d3.Ptr(),0x00FF00);
 376
 377                 Vector3d np1 = mP1 + mNormal*0.05f;
 378                 Vector3d np2 = mP2 + mNormal*0.05f;
 379                 Vector3d np3 = mP3 + mNormal*0.05f;
 380
 381                 callback->ConvexDebugTri(mP1.Ptr(), np1.Ptr(), np1.Ptr(), 0xFF00FF );
 382                 callback->ConvexDebugTri(mP2.Ptr(), np2.Ptr(), np2.Ptr(), 0xFF00FF );
 383                 callback->ConvexDebugTri(mP3.Ptr(), np3.Ptr(), np3.Ptr(), 0xFF00FF );
 384
 385                 callback->ConvexDebugPoint( np1.Ptr(), 0.01F, 0XFF00FF );
 386                 callback->ConvexDebugPoint( np2.Ptr(), 0.01F, 0XFF00FF );
 387                 callback->ConvexDebugPoint( np3.Ptr(), 0.01F, 0XFF00FF );
 388
 389 #endif
 390
 391
 392
 393   }
 394
 395   Vector3d      mP1;
 396   Vector3d      mP2;
 397   Vector3d      mP3;
 398   Vector3d mNear1;
 399   Vector3d mNear2;
 400   Vector3d mNear3;
 401   Vector3d mNormal;
 402   float           mPlaneD;
 403   float           mConcavity;
 404   float           mC1;
 405   float           mC2;
 406   float           mC3;
 407   unsigned int    mI1;
 408   unsigned int    mI2;
 409   unsigned int    mI3;
 410   int             mProcessed; // already been added...
 411 };
 412
 413 typedef std::vector< CTri > CTriVector;
 414
 415 bool featureMatch(CTri &m,const CTriVector &tris,ConvexDecompInterface *callback,const CTriVector &input_mesh)
 416 {
 417
 418   bool ret = false;
 419
 420   float neardot = 0.707f;
 421
 422   m.mConcavity = 0;
 423
 424         //gLog->Display("*********** FEATURE MATCH *************\r\n");
 425         //gLog->Display("Plane: %0.4f,%0.4f,%0.4f   %0.4f\r\n", m.mNormal.x, m.mNormal.y, m.mNormal.z, m.mPlaneD );
 426         //gLog->Display("*********************************************\r\n");
 427
 428         CTriVector::const_iterator i;
 429
 430         CTri nearest;
 431
 432
 433         for (i=tris.begin(); i!=tris.end(); ++i)
 434         {
 435                 const CTri &t = (*i);
 436
 437
 438         //gLog->Display("   HullPlane: %0.4f,%0.4f,%0.4f   %0.4f\r\n", t.mNormal.x, t.mNormal.y, t.mNormal.z, t.mPlaneD );
 439
 440                 if ( t.samePlane(m) )
 441                 {
 442                         //gLog->Display("*** PLANE MATCH!!!\r\n");
 443                         ret = false;
 444                         break;
 445                 }
 446
 447           float dot = t.mNormal.Dot(m.mNormal);
 448
 449           if ( dot > neardot )
 450           {
 451
 452       float d1 = t.planeDistance( m.mP1 );
 453       float d2 = t.planeDistance( m.mP2 );
 454       float d3 = t.planeDistance( m.mP3 );
 455
 456       if ( d1 > 0.001f || d2 > 0.001f || d3 > 0.001f ) // can't be near coplaner!
 457       {
 458
 459                 neardot = dot;
 460
 461         Vector3d n1,n2,n3;
 462
 463         t.raySect( m.mP1, m.mNormal, m.mNear1 );
 464         t.raySect( m.mP2, m.mNormal, m.mNear2 );
 465         t.raySect( m.mP3, m.mNormal, m.mNear3 );
 466
 467                                 nearest = t;
 468
 469                         ret = true;
 470                   }
 471
 472           }
 473         }
 474
 475         if ( ret )
 476         {
 477     if ( 0 )
 478     {
 479       CTriVector::const_iterator i;
 480       for (i=input_mesh.begin(); i!=input_mesh.end(); ++i)
 481       {
 482         const CTri &c = (*i);
 483         if ( c.mI1 != m.mI1 && c.mI2 != m.mI2 && c.mI3 != m.mI3 )
 484         {
 485                           c.clip( m.mP1, m.mNear1 );
 486                                   c.clip( m.mP2, m.mNear2 );
 487                                   c.clip( m.mP3, m.mNear3 );
 488         }
 489       }
 490     }
 491
 492           //gLog->Display("*********************************************\r\n");
 493         //gLog->Display("   HullPlaneNearest: %0.4f,%0.4f,%0.4f   %0.4f\r\n", nearest.mNormal.x, nearest.mNormal.y, nearest.mNormal.z, nearest.mPlaneD );
 494
 495                 m.mC1 = m.mP1.Distance( m.mNear1 );
 496                 m.mC2 = m.mP2.Distance( m.mNear2 );
 497                 m.mC3 = m.mP3.Distance( m.mNear3 );
 498
 499                 m.mConcavity = m.mC1;
 500
 501                 if ( m.mC2 > m.mConcavity ) m.mConcavity = m.mC2;
 502                 if ( m.mC3 > m.mConcavity ) m.mConcavity = m.mC3;
 503
 504     #if 0
 505                 callback->ConvexDebugTri( m.mP1.Ptr(), m.mP2.Ptr(), m.mP3.Ptr(), 0x00FF00 );
 506                 callback->ConvexDebugTri( m.mNear1.Ptr(), m.mNear2.Ptr(), m.mNear3.Ptr(), 0xFF0000 );
 507
 508                 callback->ConvexDebugTri( m.mP1.Ptr(), m.mP1.Ptr(), m.mNear1.Ptr(), 0xFFFF00 );
 509                 callback->ConvexDebugTri( m.mP2.Ptr(), m.mP2.Ptr(), m.mNear2.Ptr(), 0xFFFF00 );
 510                 callback->ConvexDebugTri( m.mP3.Ptr(), m.mP3.Ptr(), m.mNear3.Ptr(), 0xFFFF00 );
 511     #endif
 512
 513         }
 514         else
 515         {
 516                 //gLog->Display("No match\r\n");
 517         }
 518
 519         //gLog->Display("*********************************************\r\n");
 520         return ret;
 521 }
 522
 523 bool isFeatureTri(CTri &t,CTriVector &flist,float fc,ConvexDecompInterface *callback,unsigned int color)
 524 {
 525   bool ret = false;
 526
 527   if ( t.mProcessed == 0 ) // if not already processed
 528   {
 529
 530     float c = t.mConcavity / fc; // must be within 80% of the concavity of the parent.
 531
 532     if ( c > 0.85f )
 533     {
 534       // see if this triangle is a 'feature' triangle.  Meaning it shares an
 535       // edge with any existing feature triangle and is within roughly the same
 536       // concavity of the parent.
 537                         if ( flist.size() )
 538                         {
 539                           CTriVector::iterator i;
 540                           for (i=flist.begin(); i!=flist.end(); ++i)
 541                           {
 542                                   CTri &ftri = (*i);
 543                                   if ( ftri.sharesEdge(t) )
 544                                   {
 545                                           t.mProcessed = 2; // it is now part of a feature.
 546                                           flist.push_back(t); // add it to the feature list.
 547 //                                        callback->ConvexDebugTri( t.mP1.Ptr(), t.mP2.Ptr(),t.mP3.Ptr(), color );
 548                                           ret = true;
 549                                           break;
 550           }
 551                                 }
 552                         }
 553                         else
 554                         {
 555                                 t.mProcessed = 2;
 556                                 flist.push_back(t); // add it to the feature list.
 557 //                              callback->ConvexDebugTri( t.mP1.Ptr(), t.mP2.Ptr(),t.mP3.Ptr(), color );
 558                                 ret = true;
 559                         }
 560     }
 561     else
 562     {
 563       t.mProcessed = 1; // eliminated for this feature, but might be valid for the next one..
 564     }
 565
 566   }
 567   return ret;
 568 }
 569
 570 float computeConcavity(unsigned int vcount,
 571                        const float *vertices,
 572                        unsigned int tcount,
 573                        const unsigned int *indices,
 574                        ConvexDecompInterface *callback,
 575                        float *plane,      // plane equation to split on
 576                        float &volume)
 577 {
 578
 579
 580         float cret = 0;
 581         volume = 1;
 582
 583         HullResult  result;
 584   HullLibrary hl;
 585   HullDesc    desc;
 586
 587         desc.mMaxFaces = 256;
 588         desc.mMaxVertices = 256;
 589         desc.SetHullFlag(QF_TRIANGLES);
 590
 591
 592   desc.mVcount       = vcount;
 593   desc.mVertices     = vertices;
 594   desc.mVertexStride = sizeof(float)*3;
 595
 596   HullError ret = hl.CreateConvexHull(desc,result);
 597
 598   if ( ret == QE_OK )
 599   {
 600 #if 0
 601                 float bmin[3];
 602                 float bmax[3];
 603
 604                 float dx = bmax[0] - bmin[0];
 605                 float dy = bmax[1] - bmin[1];
 606                 float dz = bmax[2] - bmin[2];
 607
 608                 Vector3d center;
 609
 610                 center.x = bmin[0] + dx*0.5f;
 611                 center.y = bmin[1] + dy*0.5f;
 612                 center.z = bmin[2] + dz*0.5f;
 613 #endif
 614
 615                 volume = computeMeshVolume2( result.mOutputVertices, result.mNumFaces, result.mIndices );
 616
 617 #if 1
 618                 // ok..now..for each triangle on the original mesh..
 619                 // we extrude the points to the nearest point on the hull.
 620                 const unsigned int *source = result.mIndices;
 621
 622                 CTriVector tris;
 623
 624     for (unsigned int i=0; i<result.mNumFaces; i++)
 625     {
 626         unsigned int i1 = *source++;
 627         unsigned int i2 = *source++;
 628         unsigned int i3 = *source++;
 629
 630         const float *p1 = &result.mOutputVertices[i1*3];
 631         const float *p2 = &result.mOutputVertices[i2*3];
 632         const float *p3 = &result.mOutputVertices[i3*3];
 633
 634 //                      callback->ConvexDebugTri(p1,p2,p3,0xFFFFFF);
 635
 636                         CTri t(p1,p2,p3,i1,i2,i3); //
 637                         tris.push_back(t);
 638                 }
 639
 640     // we have not pre-computed the plane equation for each triangle in the convex hull..
 641
 642                 float totalVolume = 0;
 643
 644                 CTriVector ftris; // 'feature' triangles.
 645
 646                 const unsigned int *src = indices;
 647
 648
 649     float maxc=0;
 650
 651
 652                 if ( 1 )
 653                 {
 654       CTriVector input_mesh;
 655       if ( 1 )
 656       {
 657                     const unsigned int *src = indices;
 658                         for (unsigned int i=0; i<tcount; i++)
 659                         {
 660
 661                 unsigned int i1 = *src++;
 662                 unsigned int i2 = *src++;
 663                 unsigned int i3 = *src++;
 664
 665                 const float *p1 = &vertices[i1*3];
 666                 const float *p2 = &vertices[i2*3];
 667                 const float *p3 = &vertices[i3*3];
 668
 669                                 CTri t(p1,p2,p3,i1,i2,i3);
 670           input_mesh.push_back(t);
 671         }
 672       }
 673
 674       CTri  maxctri;
 675
 676                         for (unsigned int i=0; i<tcount; i++)
 677                         {
 678
 679                 unsigned int i1 = *src++;
 680                 unsigned int i2 = *src++;
 681                 unsigned int i3 = *src++;
 682
 683                 const float *p1 = &vertices[i1*3];
 684                 const float *p2 = &vertices[i2*3];
 685                 const float *p3 = &vertices[i3*3];
 686
 687                                 CTri t(p1,p2,p3,i1,i2,i3);
 688
 689                                 featureMatch(t, tris, callback, input_mesh );
 690
 691                                 if ( t.mConcavity > CONCAVE_THRESH )
 692                                 {
 693
 694           if ( t.mConcavity > maxc )
 695           {
 696             maxc = t.mConcavity;
 697             maxctri = t;
 698           }
 699
 700                                 float v = t.getVolume(0);
 701                                 totalVolume+=v;
 702                                 ftris.push_back(t);
 703                         }
 704
 705                         }
 706                 }
 707
 708           if ( ftris.size()  && 0 )
 709           {
 710
 711       // ok..now we extract the triangles which form the maximum concavity.
 712       CTriVector major_feature;
 713       float maxarea = 0;
 714
 715       while ( maxc > CONCAVE_THRESH  )
 716       {
 717
 718         unsigned int color = getDebugColor();  //
 719
 720         CTriVector flist;
 721
 722         bool found;
 723
 724         float totalarea = 0;
 725
 726         do
 727         {
 728           found = false;
 729           CTriVector::iterator i;
 730           for (i=ftris.begin(); i!=ftris.end(); ++i)
 731           {
 732             CTri &t = (*i);
 733             if ( isFeatureTri(t,flist,maxc,callback,color) )
 734             {
 735               found = true;
 736               totalarea+=t.area();
 737             }
 738                                 }
 739         } while ( found );
 740
 741
 742         if ( totalarea > maxarea )
 743         {
 744           major_feature = flist;
 745           maxarea = totalarea;
 746         }
 747
 748         maxc = 0;
 749
 750         for (unsigned int i=0; i<ftris.size(); i++)
 751         {
 752           CTri &t = ftris[i];
 753           if ( t.mProcessed != 2 )
 754           {
 755             t.mProcessed = 0;
 756             if ( t.mConcavity > maxc )
 757             {
 758               maxc = t.mConcavity;
 759             }
 760           }
 761         }
 762       }
 763
 764       unsigned int color = getDebugColor();
 765
 766       WpointVector list;
 767       for (unsigned int i=0; i<major_feature.size(); ++i)
 768       {
 769         major_feature[i].addWeighted(list,callback);
 770         major_feature[i].debug(color,callback);
 771       }
 772
 773       getBestFitPlane( list.size(), &list[0].mPoint.x, sizeof(Wpoint), &list[0].mWeight, sizeof(Wpoint), plane );
 774
 775                 computeSplitPlane( vcount, vertices, tcount, indices, callback, plane );
 776
 777
 778                 }
 779           else
 780           {
 781                 computeSplitPlane( vcount, vertices, tcount, indices, callback, plane );
 782           }
 783 #endif
 784
 785                 cret = totalVolume;
 786
 787           hl.ReleaseResult(result);
 788   }
 789
 790
 791         return cret;
 792 }
 793
 794
 795 }