1 /**************************************************************************** 2 * VCGLib o o * 3 * Visual and Computer Graphics Library o o * 4 * _ O _ * 5 * Copyright(C) 2004-2016 \/)\/ * 6 * Visual Computing Lab /\/| * 7 * ISTI - Italian National Research Council | * 8 * \ * 9 * All rights reserved. * 10 * * 11 * This program is free software; you can redistribute it and/or modify * 12 * it under the terms of the GNU General Public License as published by * 13 * the Free Software Foundation; either version 2 of the License, or * 14 * (at your option) any later version. * 15 * * 16 * This program is distributed in the hope that it will be useful, * 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * 20 * for more details. * 21 * * 22 ****************************************************************************/ 23 /**************************************************************************** 24 History 25 26 $Log: not supported by cvs2svn $ 27 Revision 1.15 2007/07/31 12:35:42 ganovelli 28 added ScalarType to tetra3 29 30 Revision 1.14 2006/07/06 12:39:51 ganovelli 31 adde barycenter() 32 33 Revision 1.13 2006/06/06 14:35:31 zifnab1974 34 Changes for compilation on linux AMD64. Some remarks: Linux filenames are case-sensitive. _fileno and _filelength do not exist on linux 35 36 Revision 1.12 2006/03/01 15:59:34 pietroni 37 added InterpolationParameters function 38 39 Revision 1.11 2005/12/12 11:15:26 ganovelli 40 modifications to compile with gcc 41 42 Revision 1.10 2005/11/29 16:20:33 pietroni 43 added IsInside() function 44 45 Revision 1.9 2004/10/13 12:45:51 cignoni 46 Better Doxygen documentation 47 48 Revision 1.8 2004/09/01 12:21:11 pietroni 49 minor changes to comply gcc compiler (typename's ) 50 51 Revision 1.7 2004/07/09 10:08:21 ganovelli 52 ComputeVOlume moved outside the class and other 53 minor corrections 54 55 Revision 1.6 2004/06/25 18:17:03 ganovelli 56 minor changes 57 58 Revision 1.5 2004/05/13 12:51:00 turini 59 Changed SolidAngle : table EV in table EofV 60 Changed DiedralAngle : tables FE and FV in tables FofE and FofV 61 62 Revision 1.4 2004/05/13 08:42:36 pietroni 63 the relation between entities functions are in tetra class (don't neeed template argoument) 64 65 Revision 1.3 2004/04/28 16:31:17 turini 66 Changed : 67 in SolidAngle(vind) : 68 double da0=DiedralAngle(EV(vind,0)); 69 double da1=DiedralAngle(EV(vind,1)); 70 double da2=DiedralAngle(EV(vind,2)); 71 in 72 double da0=DiedralAngle(EofV(vind,0)); 73 double da1=DiedralAngle(EofV(vind,1)); 74 double da2=DiedralAngle(EofV(vind,2)); 75 76 Changed : 77 in DiedralAngle(edgeind) : 78 int f1=FE(edgeind,0); 79 int f2=FE(edgeind,1); 80 in 81 int f1=FofE(edgeind,0); 82 int f2=FofE(edgeind,1); 83 84 Changed : 85 in DiedralAngle(edgeind) : 86 Point3d p0=FV(f1,0)->P(); 87 Point3d p1=FV(f1,1)->P(); 88 Point3d p2=FV(f1,2)->P(); 89 in 90 Point3d p0=_v[FofV(f1,0)]; 91 Point3d p1=_v[FofV(f1,1)]; 92 Point3d p2=_v[FofV(f1,2)]; 93 94 Changed : 95 in DiedralAngle(edgeind) : 96 p0=FV(f2,0)->P(); 97 p1=FV(f2,1)->P(); 98 p2=FV(f2,2)->P(); 99 in 100 p0=_v[FofV(f2,0)]; 101 p1=_v[FofV(f2,1)]; 102 p2=_v[FofV(f2,2)]; 103 104 Revision 1.2 2004/04/28 11:37:15 pietroni 105 *** empty log message *** 106 107 Revision 1.1 2004/04/22 13:19:12 ganovelli 108 first version 109 110 Revision 1.2 2004/04/20 16:26:48 pietroni 111 *** empty log message *** 112 113 Revision 1.1 2004/04/15 08:54:20 pietroni 114 *** empty log message *** 115 116 Revision 1.1 2004/04/08 01:13:31 pietroni 117 Initial commit 118 119 120 ***************************************************************************/ 121 #ifndef __VCG_TETRA3 122 #define __VCG_TETRA3 123 124 #include <vcg/space/point3.h> 125 #include <vcg/math/matrix44.h> 126 #include <vcg/math/matrix33.h> 127 128 #include <algorithm> 129 130 namespace vcg 131 { 132 /** \addtogroup space */ 133 /*@{*/ 134 /** 135 Templated class for storing a generic tetrahedron 136 137 138 */ 139 class Tetra 140 { 141 public: 142 //Tatrahedron Functions to retrieve information about relation between faces of tetrahedron(faces,adges,vertices). 143 VofE(const int & indexE,const int & indexV)144 static int VofE(const int &indexE, const int &indexV) 145 { 146 assert((indexE < 6) && (indexV < 2)); 147 static int edgevert[6][2] = {{0, 1}, 148 {0, 2}, 149 {0, 3}, 150 {1, 2}, 151 {1, 3}, 152 {2, 3}}; 153 return (edgevert[indexE][indexV]); 154 } 155 VofF(const int & indexF,const int & indexV)156 static int VofF(const int &indexF, const int &indexV) 157 { 158 assert((indexF < 4) && (indexV < 3)); 159 static int facevert[4][3] = {{0, 1, 2}, 160 {0, 3, 1}, 161 {0, 2, 3}, 162 {1, 3, 2}}; 163 return (facevert[indexF][indexV]); 164 } 165 EofV(const int & indexV,const int & indexE)166 static int EofV(const int &indexV, const int &indexE) 167 { 168 assert((indexE < 3) && (indexV < 4)); 169 static int vertedge[4][3] = {{0, 1, 2}, 170 {0, 3, 4}, 171 {5, 1, 3}, 172 {4, 5, 2}}; 173 return vertedge[indexV][indexE]; 174 } 175 EofF(const int & indexF,const int & indexE)176 static int EofF(const int &indexF, const int &indexE) 177 { 178 assert((indexF < 4) && (indexE < 3)); 179 static int faceedge[4][3] = {{0, 3, 1}, 180 {2, 4, 0}, 181 {1, 5, 2}, 182 {4, 5, 3}}; 183 return faceedge[indexF][indexE]; 184 } 185 FofV(const int & indexV,const int & indexF)186 static int FofV(const int &indexV, const int &indexF) 187 { 188 assert((indexV < 4) && (indexF < 3)); 189 static int vertface[4][3] = {{0, 1, 2}, 190 {0, 3, 1}, 191 {0, 2, 3}, 192 {1, 3, 2}}; 193 return vertface[indexV][indexF]; 194 } 195 FofE(const int & indexE,const int & indexSide)196 static int FofE(const int &indexE, const int &indexSide) 197 { 198 assert((indexE < 6) && (indexSide < 2)); 199 static int edgeface[6][2] = {{0, 1}, 200 {0, 2}, 201 {1, 2}, 202 {0, 3}, 203 {1, 3}, 204 {2, 3}}; 205 return edgeface[indexE][indexSide]; 206 } 207 VofEE(const int & indexE0,const int & indexE1)208 static int VofEE(const int &indexE0, const int &indexE1) 209 { 210 assert((indexE0 < 6) && (indexE0 >= 0)); 211 assert((indexE1 < 6) && (indexE1 >= 0)); 212 static int edgesvert[6][6] = {{-1, 0, 0, 1, 1, -1}, 213 {0, -1, 0, 2, -1, 2}, 214 {0, 0, -1, -1, 3, 3}, 215 {1, 2, -1, -1, 1, 2}, 216 {1, -1, 3, 1, -1, 3}, 217 {-1, 2, 3, 2, 3, -1}}; 218 return (edgesvert[indexE0][indexE1]); 219 } 220 VofFFF(const int & indexF0,const int & indexF1,const int & indexF2)221 static int VofFFF(const int &indexF0, const int &indexF1, const int &indexF2) 222 { 223 assert((indexF0 < 4) && (indexF0 >= 0)); 224 assert((indexF1 < 4) && (indexF1 >= 0)); 225 assert((indexF2 < 4) && (indexF2 >= 0)); 226 static int facesvert[4][4][4] = { 227 {//0 228 {-1, -1, -1, -1}, 229 {-1, -1, 0, 1}, 230 {-1, 0, -1, 2}, 231 {-1, 1, 2, -1}}, 232 {//1 233 {-1, -1, 0, 1}, 234 {-1, -1, -1, -1}, 235 {0, -1, -1, 3}, 236 {1, -1, 3, -1}}, 237 {//2 238 {-1, 0, -1, 2}, 239 {0, -1, -1, 3}, 240 {-1, -1, -1, -1}, 241 {2, 3, -1, -1}}, 242 {//3 243 {-1, 1, 2, -1}, 244 {1, -1, 3, -1}, 245 {2, 3, -1, -1}, 246 {-1, -1, -1, -1}}}; 247 return facesvert[indexF0][indexF1][indexF2]; 248 } 249 EofFF(const int & indexF0,const int & indexF1)250 static int EofFF(const int &indexF0, const int &indexF1) 251 { 252 assert((indexF0 < 4) && (indexF0 >= 0)); 253 assert((indexF1 < 4) && (indexF1 >= 0)); 254 static int facesedge[4][4] = {{-1, 0, 1, 3}, 255 {0, -1, 2, 4}, 256 {1, 2, -1, 5}, 257 {3, 4, 5, -1}}; 258 return (facesedge[indexF0][indexF1]); 259 } 260 EofVV(const int & indexV0,const int & indexV1)261 static int EofVV(const int &indexV0, const int &indexV1) 262 { 263 assert((indexV0 < 4) && (indexV0 >= 0)); 264 assert((indexV1 < 4) && (indexV1 >= 0)); 265 static int verticesedge[4][4] = {{-1, 0, 1, 2}, 266 {0, -1, 3, 4}, 267 {1, 3, -1, 5}, 268 {2, 4, 5, -1}}; 269 270 return verticesedge[indexV0][indexV1]; 271 } 272 FofVVV(const int & indexV0,const int & indexV1,const int & indexV2)273 static int FofVVV(const int &indexV0, const int &indexV1, const int &indexV2) 274 { 275 276 assert((indexV0 < 4) && (indexV0 >= 0)); 277 assert((indexV1 < 4) && (indexV1 >= 0)); 278 assert((indexV2 < 4) && (indexV2 >= 0)); 279 280 static int verticesface[4][4][4] = { 281 {//0 282 {-1, -1, -1, -1}, 283 {-1, -1, 0, 1}, 284 {-1, 0, -1, 2}, 285 {-1, 1, 2, -1}}, 286 {//1 287 {-1, -1, 0, 1}, 288 {-1, -1, -1, -1}, 289 {0, -1, -1, 3}, 290 {1, -1, 3, -1}}, 291 {//2 292 {-1, 0, -1, 2}, 293 {0, -1, -1, 3}, 294 {-1, -1, -1, -1}, 295 {2, 3, -1, -1}}, 296 {//3 297 {-1, 1, 2, -1}, 298 {1, -1, 3, -1}, 299 {2, 3, -1, -1}, 300 {-1, -1, -1, -1}}}; 301 return verticesface[indexV0][indexV1][indexV2]; 302 } 303 FofEE(const int & indexE0,const int & indexE1)304 static int FofEE(const int &indexE0, const int &indexE1) 305 { 306 assert((indexE0 < 6) && (indexE0 >= 0)); 307 assert((indexE1 < 6) && (indexE1 >= 0)); 308 static int edgesface[6][6] = {{-1, 0, 1, 0, 1, -1}, 309 {0, -1, 2, 0, -1, 2}, 310 {1, 2, -1, -1, 1, 2}, 311 {0, 0, -1, -1, 3, 3}, 312 {1, -1, 1, 3, -1, 3}, 313 {-1, 2, 2, 3, 3, -1}}; 314 315 return edgesface[indexE0][indexE1]; 316 } 317 FoppositeV(const int & indexV)318 static int FoppositeV (const int & indexV) 319 { 320 assert(indexV < 4 && indexV >= 0); 321 static int oppFaces[4] = { 3, 2, 1, 0 }; 322 323 return oppFaces[indexV]; 324 } 325 VoppositeF(const int & indexF)326 static int VoppositeF (const int & indexF) 327 { 328 assert(indexF < 4 && indexF >= 0); 329 static int oppVerts[4] = { 3, 2, 1, 0 }; 330 331 return oppVerts[indexF]; 332 } 333 EoppositeE(const int & indexE)334 static int EoppositeE (const int & indexE) 335 { 336 assert(indexE < 6 && indexE >= 0); 337 return 5 - indexE; 338 } 339 /** @brief Computes the tetrahedron barycenter 340 */ 341 template <class TetraType> Barycenter(const TetraType & t)342 static Point3<typename TetraType::ScalarType> Barycenter(const TetraType &t) 343 { 344 return ((t.cP(0) + t.cP(1) + t.cP(2) + t.cP(3)) / (typename TetraType::ScalarType)4.0); 345 } 346 347 // compute and return the volume of a tetrahedron 348 template <class TetraType> ComputeVolume(const TetraType & t)349 static typename TetraType::ScalarType ComputeVolume(const TetraType &t) 350 { 351 return (typename TetraType::ScalarType)(((t.cP(2) - t.cP(0)) ^ (t.cP(1) - t.cP(0))) * (t.cP(3) - t.cP(0)) / 6.0); 352 } 353 354 /// Returns the normal to the face face of the tetrahedron t 355 template <class TetraType> Normal(const TetraType & t,const int & face)356 static Point3<typename TetraType::ScalarType> Normal(const TetraType &t, const int &face) 357 { 358 return (((t.cP(Tetra::VofF(face, 1)) - t.cP(Tetra::VofF(face, 0))) ^ (t.cP(Tetra::VofF(face, 2)) - t.cP(Tetra::VofF(face, 0)))).Normalize()); 359 } 360 361 template <class TetraType> DihedralAngle(const TetraType & t,const size_t eidx)362 static typename TetraType::ScalarType DihedralAngle(const TetraType &t, const size_t eidx) 363 { 364 typedef typename TetraType::CoordType CoordType; 365 //get two faces incident on eidx 366 int f0 = Tetra::FofE(eidx, 0); 367 int f1 = Tetra::FofE(eidx, 1); 368 369 CoordType p0 = t.cP(Tetra::VofF(f0, 0)); 370 CoordType p1 = t.cP(Tetra::VofF(f0, 1)); 371 CoordType p2 = t.cP(Tetra::VofF(f0, 2)); 372 373 CoordType n0 = ((p2 - p0) ^ (p1 - p0)).normalized(); 374 375 p0 = t.cP(Tetra::VofF(f1, 0)); 376 p1 = t.cP(Tetra::VofF(f1, 1)); 377 p2 = t.cP(Tetra::VofF(f1, 2)); 378 379 CoordType n1 = ((p2 - p0) ^ (p1 - p0)).normalized(); 380 381 return M_PI - double(acos(n0 * n1)); 382 }; 383 384 template <class TetraType> SolidAngle(const TetraType & t,const size_t vidx)385 static typename TetraType::ScalarType SolidAngle(const TetraType &t, const size_t vidx) 386 { 387 typedef typename TetraType::ScalarType ScalarType; 388 ScalarType a0 = DihedralAngle(t, Tetra::EofV(vidx, 0)); 389 ScalarType a1 = DihedralAngle(t, Tetra::EofV(vidx, 1)); 390 ScalarType a2 = DihedralAngle(t, Tetra::EofV(vidx, 2)); 391 392 return (a0 + a1 + a2) - M_PI; 393 }; 394 395 template <class TetraType> AspectRatio(const TetraType & t)396 static typename TetraType::ScalarType AspectRatio(const TetraType &t) 397 { 398 typedef typename TetraType::ScalarType ScalarType; 399 ScalarType a0 = SolidAngle(t, 0); 400 ScalarType a1 = SolidAngle(t, 1); 401 ScalarType a2 = SolidAngle(t, 2); 402 ScalarType a3 = SolidAngle(t, 3); 403 404 return std::min(a0, std::min(a1, std::min(a2, a3))); 405 } 406 }; 407 408 /** 409 Templated class for storing a generic tetrahedron in a 3D space. 410 Note the relation with the Face class of TetraMesh complex, both classes provide the P(i) access functions to their points and therefore they share the algorithms on it (e.g. area, normal etc...) 411 */ 412 template <class ScalarType> 413 class Tetra3 : public Tetra 414 { 415 public: 416 typedef Point3<ScalarType> CoordType; 417 //typedef typename ScalarType ScalarType; 418 419 /********************************************* 420 421 **/ 422 423 private: 424 /// Vector of the 4 points that defines the tetrahedron 425 CoordType _v[4]; 426 427 public: 428 /// Shortcut per accedere ai punti delle facce P(const int j)429 inline CoordType &P(const int j) { return _v[j]; } cP(const int j)430 inline CoordType const &cP(const int j) const { return _v[j]; } 431 P0(const int j)432 inline CoordType &P0(const int j) { return _v[j]; } P1(const int j)433 inline CoordType &P1(const int j) { return _v[(j + 1) % 4]; } P2(const int j)434 inline CoordType &P2(const int j) { return _v[(j + 2) % 4]; } P3(const int j)435 inline CoordType &P3(const int j) { return _v[(j + 3) % 4]; } 436 P0(const int j)437 inline const CoordType &P0(const int j) const { return _v[j]; } P1(const int j)438 inline const CoordType &P1(const int j) const { return _v[(j + 1) % 4]; } P2(const int j)439 inline const CoordType &P2(const int j) const { return _v[(j + 2) % 4]; } P3(const int j)440 inline const CoordType &P3(const int j) const { return _v[(j + 3) % 4]; } 441 cP0(const int j)442 inline const CoordType &cP0(const int j) const { return _v[j]; } cP1(const int j)443 inline const CoordType &cP1(const int j) const { return _v[(j + 1) % 4]; } cP2(const int j)444 inline const CoordType &cP2(const int j) const { return _v[(j + 2) % 4]; } cP3(const int j)445 inline const CoordType &cP3(const int j) const { return _v[(j + 3) % 4]; } 446 447 /// compute and return the barycenter of a tetrahedron ComputeBarycenter()448 CoordType ComputeBarycenter() 449 { 450 return ((_v[0] + _v[1] + _v[2] + _v[3]) / 4); 451 } 452 453 /// compute and return the solid angle on a vertex SolidAngle(int vind)454 double SolidAngle(int vind) 455 { 456 double da0 = DiedralAngle(EofV(vind, 0)); 457 double da1 = DiedralAngle(EofV(vind, 1)); 458 double da2 = DiedralAngle(EofV(vind, 2)); 459 460 return ((da0 + da1 + da2) - M_PI); 461 } 462 463 /// compute and return the diadedral angle on an edge DiedralAngle(int edgeind)464 double DiedralAngle(int edgeind) 465 { 466 int f1 = FofE(edgeind, 0); 467 int f2 = FofE(edgeind, 1); 468 CoordType p0 = _v[FofV(f1, 0)]; 469 CoordType p1 = _v[FofV(f1, 1)]; 470 CoordType p2 = _v[FofV(f1, 2)]; 471 CoordType norm1 = ((p1 - p0) ^ (p2 - p0)); 472 p0 = _v[FofV(f2, 0)]; 473 p1 = _v[FofV(f2, 1)]; 474 p2 = _v[FofV(f2, 2)]; 475 CoordType norm2 = ((p1 - p0) ^ (p2 - p0)); 476 norm1.Normalize(); 477 norm2.Normalize(); 478 return (M_PI - acos(double(norm1 * norm2))); 479 } 480 481 /// compute and return the aspect ratio of the tetrahedron ComputeAspectRatio()482 ScalarType ComputeAspectRatio() 483 { 484 double a0 = SolidAngle(0); 485 double a1 = SolidAngle(1); 486 double a2 = SolidAngle(2); 487 double a3 = SolidAngle(3); 488 return (ScalarType)std::min(a0, std::min(a1, std::min(a2, a3))); 489 } 490 491 ///return true of p is inside tetrahedron's volume IsInside(const CoordType & p)492 bool IsInside(const CoordType &p) 493 { 494 //bb control 495 vcg::Box3<typename CoordType::ScalarType> bb; 496 for (int i = 0; i < 4; i++) 497 bb.Add(_v[i]); 498 499 if (!bb.IsIn(p)) 500 return false; 501 502 vcg::Matrix44<ScalarType> M0; 503 vcg::Matrix44<ScalarType> M1; 504 vcg::Matrix44<ScalarType> M2; 505 vcg::Matrix44<ScalarType> M3; 506 vcg::Matrix44<ScalarType> M4; 507 508 CoordType p1 = _v[0]; 509 CoordType p2 = _v[1]; 510 CoordType p3 = _v[2]; 511 CoordType p4 = _v[3]; 512 513 M0[0][0] = p1.V(0); 514 M0[0][1] = p1.V(1); 515 M0[0][2] = p1.V(2); 516 M0[1][0] = p2.V(0); 517 M0[1][1] = p2.V(1); 518 M0[1][2] = p2.V(2); 519 M0[2][0] = p3.V(0); 520 M0[2][1] = p3.V(1); 521 M0[2][2] = p3.V(2); 522 M0[3][0] = p4.V(0); 523 M0[3][1] = p4.V(1); 524 M0[3][2] = p4.V(2); 525 M0[0][3] = 1; 526 M0[1][3] = 1; 527 M0[2][3] = 1; 528 M0[3][3] = 1; 529 530 M1 = M0; 531 M1[0][0] = p.V(0); 532 M1[0][1] = p.V(1); 533 M1[0][2] = p.V(2); 534 535 M2 = M0; 536 M2[1][0] = p.V(0); 537 M2[1][1] = p.V(1); 538 M2[1][2] = p.V(2); 539 540 M3 = M0; 541 M3[2][0] = p.V(0); 542 M3[2][1] = p.V(1); 543 M3[2][2] = p.V(2); 544 545 M4 = M0; 546 M4[3][0] = p.V(0); 547 M4[3][1] = p.V(1); 548 M4[3][2] = p.V(2); 549 550 ScalarType d0 = M0.Determinant(); 551 ScalarType d1 = M1.Determinant(); 552 ScalarType d2 = M2.Determinant(); 553 ScalarType d3 = M3.Determinant(); 554 ScalarType d4 = M4.Determinant(); 555 556 // all determinant must have same sign 557 return (((d0 > 0) && (d1 > 0) && (d2 > 0) && (d3 > 0) && (d4 > 0)) || ((d0 < 0) && (d1 < 0) && (d2 < 0) && (d3 < 0) && (d4 < 0))); 558 } 559 InterpolationParameters(const CoordType & bq,ScalarType & a,ScalarType & b,ScalarType & c,ScalarType & d)560 void InterpolationParameters(const CoordType &bq, ScalarType &a, ScalarType &b, ScalarType &c, ScalarType &d) 561 { 562 const ScalarType EPSILON = ScalarType(0.000001); 563 Matrix33<ScalarType> M; 564 565 CoordType v0 = P(0) - P(2); 566 CoordType v1 = P(1) - P(2); 567 CoordType v2 = P(3) - P(2); 568 CoordType v3 = bq - P(2); 569 570 M[0][0] = v0.X(); 571 M[1][0] = v0.Y(); 572 M[2][0] = v0.Z(); 573 574 M[0][1] = v1.X(); 575 M[1][1] = v1.Y(); 576 M[2][1] = v1.Z(); 577 578 M[0][2] = v2.X(); 579 M[1][2] = v2.Y(); 580 M[2][2] = v2.Z(); 581 582 Matrix33<ScalarType> inv_M = vcg::Inverse<ScalarType>(M); 583 584 CoordType Barycentric = inv_M * v3; 585 586 a = Barycentric.V(0); 587 b = Barycentric.V(1); 588 d = Barycentric.V(2); 589 c = 1 - (a + b + d); 590 } 591 592 }; //end Class 593 594 /*@}*/ 595 } // namespace vcg 596 597 #endif 598