1 /************************************************************************* 2 * * 3 * Vega FEM Simulation Library Version 3.1 * 4 * * 5 * "distance field" library , Copyright (C) 2007 CMU, 2016 USC * 6 * All rights reserved. * 7 * * 8 * Code author: Jernej Barbic * 9 * http://www.jernejbarbic.com/code * 10 * * 11 * Research: Jernej Barbic, Hongyi Xu, Doug L. James * 12 * * 13 * Funding: National Science Foundation, Link Foundation, * 14 * Zumberge Research and Innovation Fund at USC * 15 * * 16 * This library is free software; you can redistribute it and/or * 17 * modify it under the terms of the BSD-style license that is * 18 * included with this library in the file LICENSE.txt * 19 * * 20 * This library is distributed in the hope that it will be useful, * 21 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the file * 23 * LICENSE.TXT for more details. * 24 * * 25 *************************************************************************/ 26 27 /* 28 These are the key computational routines to compute the distance field. 29 */ 30 31 #include <float.h> 32 #include "objMeshGraph.h" 33 #include "vegalong.h" 34 35 #ifdef COMPUTE_SIGNED_FIELD 36 #ifdef COMPUTE_FLOOD_FIELD computeFloodFillSignedField(ObjMesh * objMeshIn,int resolutionX_,int resolutionY_,int resolutionZ_,int maxTriCount_,int maxDepth_,int zMin_,int zMax_)37 int DistanceField::computeFloodFillSignedField(ObjMesh* objMeshIn, int resolutionX_, int resolutionY_, int resolutionZ_, int maxTriCount_, int maxDepth_,int zMin_, int zMax_ ) 38 #else 39 #ifdef COMPUTE_CLOSEST_POINT 40 int ClosestPointField::computeSignedField(ObjMesh * objMeshIn, int resolutionX_, int resolutionY_, int resolutionZ_, int maxTriCount_, int maxDepth_,int zMin_, int zMax_ ) 41 #else 42 int DistanceField::computeSignedField(ObjMesh * objMeshIn, int resolutionX_, int resolutionY_, int resolutionZ_, int maxTriCount_, int maxDepth_, int zMin_, int zMax_) 43 #endif 44 #endif 45 #else 46 #ifdef COMPUTE_CLOSEST_POINT 47 int ClosestPointField::computeUnsignedField(ObjMesh * objMeshIn, int resolutionX_, int resolutionY_, int resolutionZ_, int maxTriCount_, int maxDepth_, int zMin_, int zMax_) 48 #else 49 int DistanceField::computeUnsignedField(ObjMesh * objMeshIn, int resolutionX_, int resolutionY_, int resolutionZ_, int maxTriCount_, int maxDepth_, int zMin_, int zMax_) 50 #endif 51 #endif 52 { 53 ObjMesh objMesh(*objMeshIn); 54 55 maxTriCount = maxTriCount_; 56 maxDepth = maxDepth_; 57 resolutionX = resolutionX_; 58 resolutionY = resolutionY_; 59 resolutionZ = resolutionZ_; 60 61 #ifdef COMPUTE_SIGNED_FIELD 62 // === check if closed mesh 63 if (!(objMesh.isTriangularMesh())) 64 { 65 printf("Mesh was not triangular. Triangulating..."); fflush(NULL); 66 objMesh.triangulate(); 67 printf(" done\n"); 68 } 69 ObjMeshOrientable * objMeshOrientable = new ObjMeshOrientable(&objMesh, 0); 70 71 int numOrientationFlips = objMeshOrientable->GenerateHalfEdgeDataStructure(); 72 if (numOrientationFlips < 0) 73 { 74 printf("Error: cannot orient the mesh. Geometry is non-manifold.\n"); 75 exit(1); 76 } 77 78 cout << "Number of distinct connected components: " << objMeshOrientable->numConnectedComponents() << endl; 79 80 cout << "Checking if mesh has no boundary..." << endl; 81 82 // check if mesh has no boundary 83 if (objMeshOrientable->hasBoundary()) 84 { 85 cout << "Error: mesh has boundary. Signed distance field is ill-defined." << endl; 86 cout << " Num boundary edges: " << objMeshOrientable->numBoundaryEdges() << endl; 87 int edge = objMeshOrientable->boundaryEdge(0); 88 cout << " A boundary edge: " << objMesh.getPosition(objMeshOrientable->halfEdge(edge).startVertex()) << " " << 89 objMesh.getPosition(objMeshOrientable->halfEdge(edge).endVertex()) << endl; 90 91 return 1; 92 } 93 94 cout << "Mesh has no boundary (i.e. is closed surface)." << endl; 95 96 cout << "Checking if input mesh is oriented consistently..."; 97 if (numOrientationFlips != 0) 98 { 99 cout << " no." << endl; 100 cout << "Error: triangles in the input mesh are not oriented consistently." << endl; 101 return 2; 102 } 103 else 104 cout << " yes." << endl; 105 106 delete(objMeshOrientable); 107 #endif 108 109 ObjMeshGraph * meshGraph = NULL; 110 if (computeVoronoiDiagram) 111 { 112 meshGraph = new ObjMeshGraph(&objMesh); 113 printf("Original graph:\n"); 114 meshGraph->PrintInfo(); 115 //meshGraph->ExpandNeighbors(); 116 //printf("Expanded graph:\n"); 117 //meshGraph->PrintInfo(); 118 if (voronoiDiagram != NULL) 119 free(voronoiDiagram); 120 voronoiDiagram = (int*) malloc (sizeof(int) * (resolutionX+1) * (resolutionY+1) * (resolutionZ+1)); 121 } 122 123 // === build octree 124 125 printf("Preparing to build the octree. Max triangle count per cell: %d . Max depth: %d .\n", maxTriCount, maxDepth);fflush(NULL); 126 127 #ifdef COMPUTE_SIGNED_FIELD 128 ObjMeshOctree<TriangleWithCollisionInfoAndPseudoNormals> * objMeshOctree = 129 new ObjMeshOctree<TriangleWithCollisionInfoAndPseudoNormals>( &objMesh, maxTriCount, maxDepth ); 130 #else 131 ObjMeshOctree<TriangleWithCollisionInfo> * objMeshOctree = 132 new ObjMeshOctree<TriangleWithCollisionInfo>( &objMesh, maxTriCount, maxDepth ); 133 #endif 134 135 136 #ifndef COMPUTE_FLOOD_FIELD 137 if (useAutomaticBox && !bboxComputed) 138 { 139 computeBoundingBox(&objMesh, resolutionX, resolutionY, resolutionZ); 140 } 141 #endif 142 143 setGridParameters(); 144 zMin = (zMin_==-1) ? 0 : zMin_; 145 zMax = (zMax_==-1) ? resolutionZ : zMax_; 146 147 #ifndef COMPUTE_FLOOD_FIELD 148 //distanceData = (float*) malloc (sizeof(float)*(resolutionX+1)*(resolutionY+1)*(resolutionZ+1)); 149 vegalong size = sizeof(float)*(resolutionX+1)*(resolutionY+1)*(zMax - zMin +1); 150 distanceData = (float*) realloc (distanceData, size); 151 152 // initialize to FLT_MAX 153 for(vegalong ii=0; ii<=resolutionX; ii++) 154 for(vegalong jj=0; jj<=resolutionY; jj++) 155 for(vegalong kk=zMin; kk<=zMax; kk++) 156 { 157 vegalong index = (kk-zMin) * (resolutionX+1) * (resolutionY+1) + jj * (resolutionX+1) + ii; 158 distanceData[index] = FLT_MAX; 159 } 160 #else 161 if(zMin == 0) 162 { 163 cout << "computing the flood fill tag..." << endl; 164 computeFloodFillTag(&objMesh); 165 } 166 #endif 167 168 // debug 169 #ifdef GENERATE_DEBUG_DATA 170 pseudoData = (float*) malloc (6*sizeof(float)*(resolutionX+1)*(resolutionY+1)*(resolutionZ+1)); 171 #endif 172 173 #ifdef COMPUTE_CLOSEST_POINT 174 closestPointData = (float*) malloc (3*sizeof(float)*(resolutionX+1)*(resolutionY+1)*(resolutionZ+1)); 175 #endif 176 177 cout << "Computing the distance field..." << endl; 178 cout << "Corners of the bounding box are:" << endl; 179 cout << " " << bmin_ << endl; 180 cout << " " << bmax_ << endl; 181 cout << "Grid sizes are: " << side[0] / resolutionX << " " << side[1] / resolutionY << " " << side[2] / resolutionZ << endl; 182 183 #ifdef COMPUTE_SIGNED_FIELD 184 #ifdef COMPUTE_FLOOD_FIELD 185 ZigZagFloodFillSigned((void*)objMeshOctree, (void*)meshGraph); 186 #else 187 ZigZagSigned((void*)objMeshOctree, (void*)meshGraph); 188 #endif 189 #else 190 ZigZagUnsigned((void*)objMeshOctree, (void*)meshGraph); 191 #endif 192 193 delete(objMeshOctree); 194 delete(meshGraph); 195 196 if(zMax == resolutionZ) 197 { 198 //free(floodFillTag); 199 #ifdef COMPUTE_SIGNED_FIELD 200 cout << endl << "Signed distance field successfully computed..." << endl; 201 #else 202 cout << endl << "Unsigned distance field successfully computed..." << endl; 203 #endif 204 } 205 206 return 0; 207 } 208 209 #ifdef COMPUTE_SIGNED_FIELD 210 #ifdef COMPUTE_FLOOD_FIELD ZigZagFloodFillSigned(void * objMeshOctree_,void * meshGraph_)211 int DistanceField::ZigZagFloodFillSigned(void * objMeshOctree_, void * meshGraph_) 212 #else 213 #ifdef COMPUTE_CLOSEST_POINT 214 int ClosestPointField::ZigZagSigned(void * objMeshOctree_, void * meshGraph_) 215 #else 216 int DistanceField::ZigZagSigned(void * objMeshOctree_, void * meshGraph_) 217 #endif 218 #endif 219 #ifdef COMPUTE_FLOOD_FIELD 220 #define ZIGZAGROUTINE ZigZagFloodFillSigned 221 #else 222 #define ZIGZAGROUTINE ZigZagSigned 223 #endif 224 #else 225 #ifdef COMPUTE_CLOSEST_POINT 226 int ClosestPointField::ZigZagUnsigned(void * objMeshOctree_, void * meshGraph_) 227 #else 228 int DistanceField::ZigZagUnsigned(void * objMeshOctree_, void * meshGraph_) 229 #endif 230 #define ZIGZAGROUTINE ZigZagUnsigned 231 #endif 232 { 233 return ZIGZAGROUTINE (objMeshOctree_, meshGraph_, zMin, zMax); 234 } 235 #undef ZIGZAGROUTINE 236 237 #ifdef COMPUTE_SIGNED_FIELD 238 #ifdef COMPUTE_FLOOD_FIELD ZigZagFloodFillSigned(void * objMeshOctree_,void * meshGraph_,int zLo,int zHi,int asterisk)239 int DistanceField::ZigZagFloodFillSigned(void * objMeshOctree_, void * meshGraph_, int zLo, int zHi, int asterisk) 240 #else 241 #ifdef COMPUTE_CLOSEST_POINT 242 int ClosestPointField::ZigZagSigned(void * objMeshOctree_, void * meshGraph_, int zLo, int zHi, int asterisk) 243 #else 244 int DistanceField::ZigZagSigned(void * objMeshOctree_, void * meshGraph_, int zLo, int zHi, int asterisk) 245 #endif 246 #endif 247 #else 248 #ifdef COMPUTE_CLOSEST_POINT 249 int ClosestPointField::ZigZagUnsigned(void * objMeshOctree_, void * meshGraph_, int zLo, int zHi, int asterisk) 250 #else 251 int DistanceField::ZigZagUnsigned(void * objMeshOctree_, void * meshGraph_, int zLo, int zHi, int asterisk) 252 #endif 253 #endif 254 { 255 ObjMeshGraph * meshGraph = (ObjMeshGraph*) meshGraph_; 256 257 // triangleList will hold the triangles retrieved by range query 258 #ifdef COMPUTE_SIGNED_FIELD 259 vector< TriangleWithCollisionInfoAndPseudoNormals* > triangleList; 260 ObjMeshOctree<TriangleWithCollisionInfoAndPseudoNormals> * objMeshOctree = (ObjMeshOctree<TriangleWithCollisionInfoAndPseudoNormals> *) objMeshOctree_; 261 #else 262 vector< TriangleWithCollisionInfo* > triangleList; 263 ObjMeshOctree<TriangleWithCollisionInfo> * objMeshOctree = (ObjMeshOctree<TriangleWithCollisionInfo> *) objMeshOctree_; 264 #endif 265 266 // do the zig-zag 267 vegalong i=0, j=0, k=0; 268 k = zLo; 269 270 int diri=1; 271 int dirj=1; 272 273 double radius = len(side); 274 275 int needtoCheckFlag = 0; 276 int skipFlag = 0; 277 needtoCheckFlag = needtoCheckFlag; // just for removing warnings 278 do 279 { 280 vegalong index = (k-zMin) * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i; 281 skipFlag = 0; 282 #ifdef COMPUTE_FLOOD_FIELD 283 if(!needtoCheckFlag && floodFillTag[index] == 0) 284 skipFlag = 1; 285 #endif 286 287 float closestDistance = 0; 288 if(!skipFlag) 289 { 290 Vec3d currentPosition 291 (bmin_[0] + 1.0 * i * gridX, 292 bmin_[1] + 1.0 * j * gridY, 293 bmin_[2] + 1.0 * k * gridZ); 294 295 //printf("%d %d %d\n",i,j,k); 296 // process grid point i,j,k 297 triangleList.clear(); 298 299 while (triangleList.size() <= 0) // we should only go through this loop exactly once 300 { 301 //printf("radius: %f pos: %f %f %f\n", radius, currentPosition[0], currentPosition[1], currentPosition[2]); 302 objMeshOctree->rangeQuery(triangleList, SimpleSphere(currentPosition, radius)); 303 304 if (triangleList.size() <= 0) 305 { // should never happen... but to stay robust 306 cout << "Warning: range query didn't find any triangles. Increasing radius by a factor of 2 and re-trying." << endl; 307 radius *= 2; 308 } 309 310 //if (triangleList.size() > 5) 311 //makeUniqueList(triangleList); 312 } 313 314 // find closest triangle among the retrieved ones 315 // initialization: 316 double closestDistance2 = 1.25 * radius * radius; // there will be somebody within that radius (actually even without the factor "1.25") 317 // (factor "1.25" added to account for numerical round-off) 318 319 int closestFeature = -1; 320 int closestTriangle = -1; 321 closestTriangle = closestTriangle; // to avoid compiler warnings 322 int indexClosestTriangle = -1; 323 324 #ifdef COMPUTE_CLOSEST_POINT 325 Vec3d closestPosition(DBL_MAX, DBL_MAX, DBL_MAX); 326 #endif 327 328 for (unsigned int l=0; l<triangleList.size(); l++) 329 { 330 int closestLocalFeature = -1; 331 332 #ifdef COMPUTE_CLOSEST_POINT 333 double alpha, beta, gamma; 334 double d2 = triangleList[l]->distanceToPoint2(currentPosition, &closestLocalFeature, &alpha, &beta, &gamma); 335 if (d2 < closestDistance2) 336 { 337 closestDistance2 = d2; 338 closestPosition = triangleList[l]->getBarycentricLocation(alpha, beta, gamma); 339 closestFeature = closestLocalFeature; 340 closestTriangle = l; 341 indexClosestTriangle = triangleList[l]->index(); 342 } 343 #else 344 double d2 = triangleList[l]->distanceToPoint2(currentPosition,&closestLocalFeature); 345 if (d2 < closestDistance2) 346 { 347 closestDistance2 = d2; 348 closestFeature = closestLocalFeature; 349 closestTriangle = l; 350 indexClosestTriangle = triangleList[l]->index(); 351 } 352 #endif 353 } 354 355 // now, indexClosestTriangle and closestFeature are set to the closest triangle and the closest feature on the triangle 356 357 if (computeVoronoiDiagram) 358 { 359 int closestSite = meshGraph->graphID(indexClosestTriangle, closestFeature); 360 voronoiDiagram[k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i] = closestSite; 361 //printf("cs: %d cf: %d\n", closestSite, closestFeature); 362 } 363 364 if (closestFeature < 0) // should never happen 365 { 366 cout << "Internal error: did not find any triangle within the guaranteed radius:" << radius << endl; 367 return 3; 368 } 369 370 #ifdef GENERATE_DEBUG_DATA 371 if ((i==34) && (j==17) && ((k==44) || (k==45))) 372 { 373 printf("Closest index: %d Feature: %d Min dist2: %.15f\n\n", indexClosestTriangle, closestFeature, closestDistance2); 374 } 375 if ((i==34) && (j==17) && ((k==44) || (k==45))) 376 printf("Grid location: %.15f %.15f %.15f\n", currentPosition[0], currentPosition[1], currentPosition[2]); 377 #endif 378 379 // square root... 380 closestDistance = (float) sqrt(closestDistance2); 381 382 #ifdef COMPUTE_SIGNED_FIELD 383 // determine sign, as in [Baerentzen 2002] 384 Vec3d pseudoNormal = triangleList[closestTriangle]->pseudoNormal(closestFeature); 385 Vec3d pseudoClosestPosition = triangleList[closestTriangle]->pseudoClosestPosition(closestFeature); 386 387 #ifdef COMPUTE_FLOOD_FIELD 388 if (dot(pseudoNormal,currentPosition-pseudoClosestPosition) <= 0) // interior, must flip sign 389 { 390 closestDistance *= -1; 391 needtoCheckFlag = 1; 392 } 393 else 394 needtoCheckFlag = 0; 395 #else 396 if (dot(pseudoNormal,currentPosition-pseudoClosestPosition) < 0) // interior, must flip sign 397 closestDistance *= -1; 398 #endif 399 //printf("%G ", closestDistance); 400 #endif 401 402 // register result 403 distanceData[index] = closestDistance; 404 405 #ifdef COMPUTE_CLOSEST_POINT 406 vegalong dataIndex = 3 * (k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i); 407 closestPointData[dataIndex+0] = closestPosition[0]; 408 closestPointData[dataIndex+1] = closestPosition[1]; 409 closestPointData[dataIndex+2] = closestPosition[2]; 410 #endif 411 412 // store debug info to disk file 413 #ifdef GENERATE_DEBUG_DATA 414 //pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+0] = pseudoNormal[0]; 415 //pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+1] = pseudoNormal[1]; 416 //pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+2] = pseudoNormal[2]; 417 //pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+3] = pseudoClosestPosition[0]; 418 //pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+4] = pseudoClosestPosition[1]; 419 //pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+5] = pseudoClosestPosition[2]; 420 pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+0] = 0; 421 pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+1] = 0; 422 pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+2] = 0; 423 pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+3] = closestPosition[0]; 424 pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+4] = closestPosition[1]; 425 pseudoData[6*(k * (resolutionX+1) * (resolutionY+1) + j * (resolutionX+1) + i)+5] = closestPosition[2]; 426 #endif 427 } 428 429 closestDistance = distanceData[index]; 430 431 vegalong oldi = i; 432 vegalong oldj = j; 433 vegalong oldk = k; 434 435 // increment i,j,k 436 i += diri; 437 438 if (i>resolutionX) 439 { 440 diri = -1; // start going to the left 441 i = resolutionX; 442 j += dirj; 443 } 444 445 if (i<0) 446 { 447 diri = +1; // start going to the right 448 i = 0; 449 j += dirj; 450 } 451 452 if (j>resolutionY) 453 { 454 dirj = -1; // start going down 455 j = resolutionY; 456 if (asterisk) 457 cout << "*"; 458 cout << k << " " << flush; 459 k += 1; 460 } 461 462 if (j<0) 463 { 464 dirj = +1; // start going up 465 j = 0; 466 if (asterisk) 467 cout << "*"; 468 cout << k << " " << flush; 469 k += 1; 470 } 471 472 // spatial coherence: update radius for next iteration 473 double delta=0.0; // this will be modified below 474 if (i != oldi) 475 delta = 1.01 * gridX; 476 477 if (j != oldj) 478 delta = 1.01 * gridY; 479 480 if (k != oldk) 481 delta = 1.01 * gridZ; 482 483 radius = fabs(closestDistance) + delta; 484 } 485 while (k <= zHi); 486 487 return 0; 488 } 489 490