开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > OpenGL > 查看源码
Utils3d.cpp
资源名称：display_list.zip　[点击查看]
上传用户：hcfgz168
上传日期：2011-09-11
资源大小：116k
文件大小：26k
源码类别：

OpenGL

开发平台：

WINDOWS

Utils3d.cpp：源码内容
  1. //********************************************
  2. // Utils3d.cpp
  3. //********************************************
  4. // pierre.alliez@cnet.francetelecom.fr
  5. // Created : 29/01/97
  6. // Modified : 14/09/98
  7. //********************************************
  9. #include "stdafx.h"
  10. #include "math.h"
  11. #include "Utils3d.h"
  14. int round(double x)
  15. {
  16. int _ceil = (int)ceil(x);
  17. int _floor = (int)floor(x);
  18. if((_ceil-x) < (x-_floor))
  19. return _ceil;
  20. else
  21. return _floor;
  22. }
  24. //**************************************
  25. // SinAngle between two faces
  26. // u  ^ v  = w
  27. // |w| = |u||v| |sin(u,v)|
  28. //**************************************
  29. double SinAngle(CFace3d *pFace1,
  30.   CFace3d *pFace2)
  31. {
  32. ASSERT(pFace1 != NULL);
  33. ASSERT(pFace2 != NULL);
  35. pFace1->CalculateNormal();
  36. pFace2->CalculateNormal();
  38. if(!pFace1->HasNormal() || !pFace2->HasNormal())
  39. return 0.0f;
  41. // Avoid illegal call of non-static function
  42. CVector3d w = CVector3d::Inner(pFace1->GetNormal(),pFace2->GetNormal());
  43. double Norm1 = pFace1->GetNormal()->GetNormL2();
  44. double Norm2 = pFace2->GetNormal()->GetNormL2();
  45. double Norm_w = w.GetNormL2();
  47. if(Norm1 * Norm2)
  48. return (Norm_w / (Norm1 * Norm2));
  49. else
  50. {
  51. TRACE("** Norms : %g %gn",Norm1,Norm2);
  52. return 0.0;
  53. }
  54. }
  57. //**************************************
  58. // CosAngle between two faces
  59. // uv = |u||v| cos(u,v)
  60. //**************************************
  61. double CosAngle(CFace3d *pFace1,
  62.   CFace3d *pFace2)
  63. {
  64. ASSERT(pFace1 != NULL);
  65. ASSERT(pFace2 != NULL);
  67. pFace1->CalculateNormal();
  68. pFace2->CalculateNormal();
  70. if(!pFace1->HasNormal() || !pFace2->HasNormal())
  71. return 0.0f;
  73. CVector3d *pU = pFace1->GetNormal();
  74. CVector3d *pV = pFace2->GetNormal();
  75. double NormU = pU->GetNormL2();
  76. double NormV = pV->GetNormL2();
  77. double product = NormU*NormV;
  79. // Check
  80. if(product == 0)
  81. {
  82. TRACE("** Norms : %g %gn",NormU,NormV);
  83. return 0.0;
  84. }
  86. double scalar = ::Scalar(pU,pV);
  88. return (scalar / product);
  89. }
  91. //**************************************
  92. // Angle between two faces (in radians)
  93. // we use this formula
  94. // uv = |u||v| cos(u,v)
  95. // u  ^ v  = w
  96. // |w| = |u||v| |sin(u,v)|
  97. //**************************************
  98. double Angle(CFace3d *pFace1,
  99.  CFace3d *pFace2)
  100. {
  101. ASSERT(pFace1 != NULL);
  102. ASSERT(pFace2 != NULL);
  104. pFace1->CalculateNormal();
  105. pFace2->CalculateNormal();
  107. if(!pFace1->HasNormal() || !pFace2->HasNormal())
  108. return 0.0f;
  110. CVector3d *pU = pFace1->GetNormal();
  111. CVector3d *pV = pFace2->GetNormal();
  113. return Angle(pU,pV);
  114. }
  116. //**************************************
  117. // Angle between two vectors (in radians)
  118. // we use this formula
  119. // uv = |u||v| cos(u,v)
  120. // u  ^ v  = w
  121. // |w| = |u||v| |sin(u,v)|
  122. //**************************************
  123. double Angle(CVector3d *pU,
  124.  CVector3d *pV)
  125. {
  126. double NormU = pU->GetNormL2();
  127. double NormV = pV->GetNormL2();
  128. double product = NormU*NormV;
  130. // Check
  131. if(product == 0)
  132. {
  133. TRACE("** Norms : %g %gn",NormU,NormV);
  134. return 0.0;
  135. }
  137. double scalar = ::Scalar(pU,pV);
  138. ASSERT(product != 0);
  139. double cosinus = scalar / product;
  141. // Sinus
  142. CVector3d w = CVector3d::Inner(pU,pV);
  143. double Norm_w = w.GetNormL2();
  145. double AbsSinus = Norm_w / product;
  147. // Remove degeneracy
  148. AbsSinus = (AbsSinus > 1) ? 1 : AbsSinus;
  149. AbsSinus = (AbsSinus < -1) ? -1 : AbsSinus;
  150. if(AbsSinus > 1)
  151. {
  152. TRACE("Remove degeneracy (AbsSinus : %lg)n",AbsSinus);
  153. }
  155. if(cosinus >= 0)
  156. return asin(AbsSinus);
  157. else
  158. return (PI-asin(AbsSinus));
  159. }
  162. //**************************************
  163. // AddFaceNoDuplicate
  164. //**************************************
  165. int AddFaceNoDuplicate(CArray3d<CFace3d> &array,
  166.  CFace3d *pFace)
  167. {
  168. int size = array.GetSize();
  169. for(int i=0;i<size;i++)
  170. {
  171. if(array[i] == pFace)
  172. return 0;
  173. }
  174. array.Add(pFace);
  175. //TRACE("** AddFaceNoDuplicate : %xn",pFace);
  177. return 1;
  178. }
  180. //**************************************
  181. // AddFaceNoDuplicate
  182. //**************************************
  183. int AddFaceRecursive(CArray3d<CFace3d> &array,
  184.  CFace3d *pFace,
  185.  CVertex3d *pVertex)
  186. {
  187. // Add vertex, if needed
  188. if(pFace->HasVertex(pVertex))
  189. {
  190. if(!AddFaceNoDuplicate(array,pFace))
  191. return 0;
  192. }
  193. else
  194. return 0;
  196. // Search in neighbors
  197. for(int i=0;i<3;i++)
  198. {
  199. CFace3d *pFaceNeighbor = pFace->f(i);
  200. if(pFaceNeighbor != NULL)
  201. if(pFaceNeighbor->HasVertex(pVertex))
  202. AddFaceRecursive(array,pFaceNeighbor,pVertex);
  203. }
  205. return 1;
  206. }
  209. //**************************************
  210. // NbSharpEdge
  211. //**************************************
  212. int NbSharpEdge(CArray3d<CFace3d> &array,
  213. double threshold)
  214. {
  215. int size = array.GetSize();
  217. // Check 
  218. if(size == 0)
  219. return 0;
  221. // Only one face, 2 sharp edges
  222. if(size == 1)
  223. return 2;
  225. int NbSharp = 0;
  227. for(int i=0;i<size-1;i++)
  228. {
  229. CFace3d *pFace = array[i];
  230. for(int j=i+1;j<size;j++)
  231. {
  232. CFace3d *pFaceTest = array[j];
  233. if(pFace->HasNeighbor(pFaceTest))
  234. NbSharp += (SinAngle(pFace,pFaceTest) >= threshold);
  235. }
  236. }
  238. return NbSharp;
  239. }
  242. //**************************************
  243. // NormalSum
  244. //**************************************
  245. int NormalSum(CArray3d<CFace3d> &array,
  246. double *pSum)
  247. {
  248. int size = array.GetSize();
  250. // Check 
  251. if(size == 0)
  252. return 0;
  254. // Only one face
  255. if(size == 1)
  256. return 0;
  258. double sum = 0;
  259. for(int i=0;i<size-1;i++)
  260. {
  261. CFace3d *pFace = array[i];
  262. for(int j=i+1;j<size;j++)
  263. {
  264. CFace3d *pFaceTest = array[j];
  265. if(pFace->HasNeighbor(pFaceTest))
  266. sum += SinAngle(pFace,pFaceTest);
  267. }
  268. }
  270. *pSum = sum;
  271. //TRACE("Value : %gn",sum);
  272. return 1;
  273. }
  275. //**************************************
  276. // NormalMax
  277. // (in radians)
  278. //**************************************
  279. int NormalMax(CArray3d<CFace3d> &array,
  280. double *pMax)
  281. {
  282. int size = array.GetSize();
  284. // Check 
  285. if(size == 0)
  286. return 0;
  288. // Only one face
  289. if(size == 1)
  290. return 0;
  292. double max = 0;
  293. for(int i=0;i<size-1;i++)
  294. {
  295. CFace3d *pFace = array[i];
  296. for(int j=i+1;j<size;j++)
  297. {
  298. CFace3d *pFaceTest = array[j];
  299. if(pFace->HasNeighbor(pFaceTest))
  300. {
  301. double sin = SinAngle(pFace,pFaceTest);
  302. max = (sin > max) ? sin : max;
  303. }
  304. }
  305. }
  307. *pMax = max;
  308. max = (max > 1) ? 1 : max; // precision work-around
  309. //TRACE("Value : %gn",max);
  311. return 1;
  312. }
  315. int MaxAngleBetweenFaces(CArray3d<CFace3d> &array,
  316.  double *pMax)
  317. {
  318. int size = array.GetSize();
  319. // Check 
  320. if(size == 0)
  321. return 0;
  322. // Only one face
  323. if(size == 1)
  324. return 0;
  325. double max = 0;
  326. for(int i=0;i<size-1;i++)
  327. {
  328. CFace3d *pFace = array[i];
  329. for(int j=i+1;j<size;j++)
  330. {
  331. CFace3d *pFaceTest = array[j];
  332. if(pFace->HasNeighbor(pFaceTest))
  333. {
  334. double angle = ::Angle(pFace,pFaceTest);
  335. max = (angle > max) ? angle : max;
  336. }
  337. }
  338. }
  340. *pMax = max;
  341. return 1;
  342. }
  345. //**************************************
  346. // DistanceSquare
  347. //**************************************
  348. double DistanceSquare(CVertex3d *pVertex1,
  349.  CVertex3d *pVertex2)
  350. {
  351. return ( (double)(pVertex1->x() - pVertex2->x()) * 
  352.        (double)(pVertex1->x() - pVertex2->x()) +
  353.    (double)(pVertex1->y() - pVertex2->y()) * 
  354.  (double)(pVertex1->y() - pVertex2->y()) +
  355.    (double)(pVertex1->z() - pVertex2->z()) * 
  356.  (double)(pVertex1->z() - pVertex2->z()) );
  357. }
  359. //**************************************
  360. // Distance
  361. //**************************************
  362. double Distance(CVertex3d *pVertex1,
  363.   CVertex3d *pVertex2)
  364. {
  365. return sqrt(DistanceSquare(pVertex1,pVertex2));
  366. }
  368. //**************************************
  369. // DistanceVertexToLine
  370. // Line go through Va and Vb
  371. // Projected vertex is stored in 
  372. // pVertexProjected
  373. //**************************************
  374. double DistanceVertexToLine(CVertex3d *pVertex,
  375.  CVertex3d *pVa,
  376.  CVertex3d *pVb,
  377.  CVertex3d *pVertexProjected)
  378. {
  379. CVector3d VectorAC(pVa,pVertex);
  380. CVector3d VectorAB(pVa,pVb);
  381. double DistanceSquareAB = DistanceSquare(pVa,pVb);
  382. double temp = Scalar(&VectorAC,&VectorAB) / DistanceSquareAB;
  383. VectorAB *= (float)temp;
  384. pVertexProjected->Set(pVa);
  385. pVertexProjected->Move(&VectorAB);
  386. return Distance(pVertexProjected,pVertex);
  387. }
  389. //**************************************
  390. // DistanceVertexToLine
  391. // Line go through Va and Vb
  392. // Projected vertex is stored in 
  393. // pVertexProjected
  394. //**************************************
  395. int ProjectVertexOnLine(CVertex3d *pVertex,
  396.    CVertex3d *pVa,
  397.  CVertex3d *pVb,
  398.  CVertex3d *pVertexProjected)
  399. {
  400. CVector3d VectorAC(pVa,pVertex);
  401. CVector3d VectorAB(pVa,pVb);
  402. double DistanceSquareAB = DistanceSquare(pVa,pVb);
  403. if(DistanceSquareAB == 0.0f)
  404. return 0; // no line
  405. double temp = Scalar(&VectorAC,&VectorAB) / DistanceSquareAB;
  406. VectorAB *= (float)temp;
  407. pVertexProjected->Set(pVa);
  408. pVertexProjected->Move(&VectorAB);
  409. return 1;
  410. }
  413. //**************************************
  414. // DistanceVertexToFace
  415. // Projected vertex is stored in 
  416. // pVertexProjected
  417. //**************************************
  418. double SquareDistanceVertexToFaceIfInf(CVertex3d *pVertex,
  419. CFace3d *pFace,
  420. CVertex3d *pVertexProjected,
  421. double OldDistance)
  422. {
  423. // Project point on triangle's plane
  424. ProjectPointOnFace(pVertex,pFace,pVertexProjected);
  426. if(DistanceSquare(pVertexProjected,pVertex) > OldDistance && 
  427.  !VertexInFace(pVertexProjected,pFace))
  428. return -1.0; // flag
  429. else
  430. return SquareDistanceVertexToFace(pVertex,pFace,pVertexProjected);
  431. }
  433. //**************************************
  434. // DistanceVertexToFace
  435. // Projected vertex is stored in 
  436. // pVertexProjected
  437. //**************************************
  438. double SquareDistanceVertexToFace(CVertex3d *pVertex,
  439.  CFace3d *pFace,
  440.  CVertex3d *pVertexProjected)
  441. {
  442. // ** DEBUG **
  444. // Project point on triangle's plane
  445. ProjectPointOnFace(pVertex,pFace,pVertexProjected);
  447. // Is point in triangle  ?
  448. if(VertexInFace(pVertexProjected,pFace))
  449. {
  450. // TRACE(" vertex in face... (projected : %g %g %g)n",pVertexProjected->x(),pVertexProjected->y(),pVertexProjected->z());
  451. return DistanceSquare(pVertexProjected,pVertex);
  452. }
  454. // Out, also project "projected vertex" on 3 lines
  455. CVertex3d VertexProjectedOnLine[3];
  456. for(int i=0;i<3;i++)
  457. {
  458. ProjectVertexOnLine(pVertexProjected,pFace->v(i),pFace->v((i+1)%3),&VertexProjectedOnLine[i]);
  459. //// TRACE("  ProjectVertexOnLine : (%g %g %g)...n",VertexProjectedOnLine[i].x(),VertexProjectedOnLine[i].y(),VertexProjectedOnLine[i].z());
  460. }
  462. // Is point on segment of face ?
  463. int success = 0;
  464. for(i=0;i<3;i++)
  465. if(VertexOnSegment(&VertexProjectedOnLine[i],pFace->v(i),pFace->v((i+1)%3)))
  466. {
  467. success = 1;
  468. VertexProjectedOnLine[i].SetFlag(1);
  469. //// TRACE("  vertex on segment of face (%g %g %g)...n",VertexProjectedOnLine[i].x(),VertexProjectedOnLine[i].y(),VertexProjectedOnLine[i].z());
  470. }
  472. // At least one vertex on segment
  473. if(success)
  474. {
  475. // Compare with projected points in segment and vertices of triangle
  477. // Get nearest vertex
  478. CVertex3d *pVertexTmp = pFace->FindNearestVertex(pVertexProjected);
  479. CVertex3d VertexProjectedFinal(pVertexTmp);
  480. double MinDistance = DistanceSquare(pVertexProjected,pVertexTmp);
  482. for(i=0;i<3;i++)
  483. if(VertexProjectedOnLine[i].GetFlag())
  484. {
  485. VertexProjectedOnLine[i].SetFlag(0); // Restore
  486. double tmp = DistanceSquare(pVertexProjected,&VertexProjectedOnLine[i]);
  487. if(tmp <= MinDistance)
  488. {
  489. VertexProjectedFinal.Set(VertexProjectedOnLine[i]);
  490. MinDistance = tmp;
  491. }
  492. }
  493. pVertexProjected->Set(VertexProjectedFinal);
  494. // TRACE("  nearest vertex on face (%g %g %g)...n",VertexProjectedFinal.x(),VertexProjectedFinal.y(),VertexProjectedFinal.z());
  495. return DistanceSquare(&VertexProjectedFinal,pVertex);
  496. }
  498. // Out, we must search nearest vertex
  499. CVertex3d *pVertexProjectedFinal = pFace->FindNearestVertex(pVertexProjected);
  501. pVertexProjected->Set(pVertexProjectedFinal);
  502. // TRACE("  nearest vertex on face (%g %g %g)...n",pVertexProjectedFinal->x(),pVertexProjectedFinal->y(),pVertexProjectedFinal->z());
  503. return DistanceSquare(pVertexProjectedFinal,pVertex);
  504. }
  507. //********************************************
  508. // SquareDistanceToFaceFast
  509. //********************************************
  510. // Search min distance between pVertex and 
  511. // mesh's faces around pFace
  512. //********************************************
  513. double SquareDistanceToSetOfFace(CVertex3d *pVertex,
  514. CVertex3d *pVertexRef,
  515. CFace3d **ppFace,
  516. CArray3d<CFace3d> &ArrayFace)
  517. {
  518. ASSERT(pVertex != NULL);
  519. ASSERT(pVertexRef != NULL);
  521. int size = ArrayFace.GetSize();
  522. ASSERT(size != 0);
  524. // For each face in set
  525. CVertex3d vertex;
  526. double distance = MAX_DOUBLE;
  528. for(int i=0;i<size;i++)
  529. {
  530. double tmp = SquareDistanceVertexToFaceIfInf(pVertex,ArrayFace[i],&vertex,distance);
  531. if(tmp < distance && tmp != -1.0f)
  532. {
  533. distance = tmp;
  534. pVertexRef->Set(vertex);
  535. *ppFace = ArrayFace[i];
  536. }
  537. }
  539. //TRACE("Dist : %g vertex(%g,%g,%g)",distance,pVertex->x(),pVertex->y(),pVertex->z());
  541. return distance;
  542. }
  544. //********************************************
  545. // MaxSquareDistanceVertexFace
  546. // For each vertex, search smallest distance 
  547. // to set of faces (pArrayFace)
  548. // Return max founded
  549. //********************************************
  550. double MaxSquareDistanceVertexFace(CArray3d<CVertex3d> *pArrayVertex,
  551.      CArray3d<CFace3d> *pArrayFace,
  552. CVertex3d **ppVertex,
  553. CFace3d **ppFace)
  554. {
  555. int NbVertex = pArrayVertex->GetSize();
  556. int NbFace = pArrayFace->GetSize();
  558. // Check
  559. ASSERT(NbVertex != 0 && NbFace != 0);
  561. // Non-optimized
  562. CVertex3d v;
  563. double max = 0.0;
  564. double distance;
  565. for(int i=0;i<NbVertex;i++)
  566. {
  567. // Get distance to set of faces
  568. distance = MAX_FLOAT;
  569. for(int j=0;j<NbFace;j++)
  570. {
  571. double d = ::SquareDistanceVertexToFace(pArrayVertex->GetAt(i),
  572.                                       pArrayFace->GetAt(j),&v);
  573. if(d < distance)
  574. distance = d;
  575. }
  576. // Get max for each vertex
  577.   if(distance > max)
  578. {
  579. max = distance;
  580. if(ppFace != NULL && ppVertex != NULL)
  581. {
  582. *ppFace = pArrayFace->GetAt(j);
  583. *ppVertex = pArrayVertex->GetAt(i);
  584. }
  585. }
  586. }
  587. return max;
  588. }
  591. //**************************************
  592. // Scalar product
  593. //**************************************
  594. double Scalar(CVector3d *pV1,
  595. CVector3d *pV2)
  596. {
  597. return ((double)pV1->x()*(double)pV2->x() +  
  598.       (double)pV1->y()*(double)pV2->y() + 
  599. (double)pV1->z()*(double)pV2->z());
  600. }
  602. //**************************************
  603. // Opposite per coord
  604. //**************************************
  605. int OppositePerCoord(CVector3d *pV1,
  606.  CVector3d *pV2)
  607. {
  608. return (pV1->x()*pV2->x() < 0.0f &&  
  609.       pV1->y()*pV2->y() < 0.0f &&  
  610. pV1->z()*pV2->z() < 0.0f);
  611. }
  613. //**************************************
  614. // Scalar product
  615. //**************************************
  616. double Scalar(CVector3d &v1,
  617. CVector3d &v2)
  618. {
  619. return ((double)v1.x() * (double)v2.x() +  
  620.       (double)v1.y() * (double)v2.y() + 
  621. (double)v1.z() * (double)v2.z());
  622. }
  624. //********************************************
  625. // Inner
  626. //********************************************
  627. CVector3d Inner(CVector3d* u,
  628. CVector3d* v)
  629. {
  630. // w = u ^ v
  631. CVector3d w;
  632. w.Set((float)((double)u->y()*(double)v->z()-(double)u->z()*(double)v->y()),
  633.     (float)((double)u->z()*(double)v->x()-(double)u->x()*(double)v->z()),
  634.     (float)((double)u->x()*(double)v->y()-(double)u->y()*(double)v->x()));
  635. return w;
  636. }
  640. //**************************************
  641. // PlaneEquation
  642. // Plane : ax + by + cz + d = 0
  643. // Normal : a,b,c
  644. // d : distance from origin along normal 
  645. //     to the plane
  646. //**************************************
  647. int PlaneEquationFromFace(CFace3d *pFace,
  648. float *a,
  649. float *b,
  650. float *c,
  651. float *d)
  652. {
  653. ASSERT(a != NULL);
  654. ASSERT(b != NULL);
  655. ASSERT(c != NULL);
  656. ASSERT(d != NULL);
  658. pFace->CalculateNormal();
  659. if(!pFace->HasNormal())
  660. return 0;
  662. CVector3d *pNormal = pFace->GetNormal();
  663. CVertex3d *pVertex = pFace->v1();
  665. ASSERT(pNormal != NULL);
  666. ASSERT(pVertex != NULL);
  668. // Store coeff
  669. *a = pNormal->x();
  670. *b = pNormal->y();
  671. *c = pNormal->z();
  672. *d = -(pNormal->x()*pVertex->x()+
  673.      pNormal->y()*pVertex->y()+
  674.  pNormal->z()*pVertex->z());
  676. return 1;
  677. }
  679. //**************************************
  680. // VertexOnSegment
  681. //**************************************
  682. int VertexOnSegment(CVertex3d *pVertex,
  683. CVertex3d *pV1,
  684. CVertex3d *pV2)
  685. {
  686. CVector3d v1(pV1,pVertex);
  687. CVector3d v2(pV2,pVertex);
  688. return (Scalar(v1,v2) <= 0);
  689. }
  691. //**************************************
  692. // ProjectPointOnPlane
  693. // Plane : ax + by + cz + d = 0
  694. // We store projected vertex coordinates 
  695. // in pVertexProjected
  696. //**************************************
  697. int ProjectPointOnPlane(CVertex3d *pVertex,
  698. CVertex3d *pVertexProjected,
  699. float a,
  700. float b,
  701. float c,
  702. float d)
  703. {
  704. // We need valid plane
  705. float div = a*a+b*b+c*c;
  706. ASSERT(div != 0.0f);
  707. if(div == 0.0f)
  708. return 0;
  710. float x = pVertex->x();
  711. float y = pVertex->y();
  712. float z = pVertex->z();
  713. float t = -(a*x + b*y + c*z + d) / div;
  714. pVertexProjected->Set(a*t+x,b*t+y,c*t+z);
  716. return 1;
  717. }
  719. //**************************************
  720. // ProjectPointOnFace
  721. // We store projected vertex coordinates 
  722. // in pVertexProjected
  723. //**************************************
  724. int ProjectPointOnFace(CVertex3d *pVertex,
  725. CFace3d *pFace,
  726. CVertex3d *pVertexProjected)
  727. {
  728. // Get plane coeff
  729. float a,b,c,d;
  730. if(PlaneEquationFromFace(pFace,&a,&b,&c,&d))
  731. return ProjectPointOnPlane(pVertex,pVertexProjected,a,b,c,d);
  732. else
  733. return 0;
  734. //TRACE("** Plane equation : %f %f %f %fn",a,b,c,d);
  735. // Project point on triangle's plane
  736. }
  739. //**************************************
  740. // IsVertexInFace
  741. // We know vertex is in plane of
  742. // triangle, and just test whenever
  743. // vertex V is in triangular face F
  744. //**************************************
  745. int VertexInFace(CVertex3d *pVertex,
  746.  CFace3d *pFace)
  747. {
  748. CVector3d vv1(pVertex,pFace->v1());
  749. CVector3d vv2(pVertex,pFace->v2());
  750. CVector3d vv3(pVertex,pFace->v3());
  752. CVector3d w1 = CVector3d::Inner(&vv1,&vv2);
  753. CVector3d w2 = CVector3d::Inner(&vv2,&vv3);
  754. CVector3d w3 = CVector3d::Inner(&vv3,&vv1);
  756. return (Scalar(&w1,&w2) >= 0 && 
  757.       Scalar(&w2,&w3) >= 0 && 
  758. Scalar(&w1,&w3) >= 0);
  759. }
  762. //**************************************
  763. // Spring term 
  764. // SIGMA(j,k)(k|Vj - Vk|^2);
  765. //**************************************
  766. double Spring(CVertex3d *pVertex,
  767.   double k)
  768. {
  769. // Check
  770. int NbVertexNeighbor = pVertex->NbVertexNeighbor();
  771. if(NbVertexNeighbor == 0)
  772. return 0.0f;
  774. double sum = 0.0f;
  775. for(int i=0;i<NbVertexNeighbor;i++)
  776. sum += DistanceSquare(pVertex,pVertex->GetVertexNeighbor(i));
  778. return (sum*k);
  779. }
  781. //**************************************
  782. // Area
  783. //**************************************
  784. double Area(CVertex3d *pV0,
  785. CVertex3d *pV1,
  786. CVertex3d *pV2)
  787. {
  788. ASSERT(pV0 != NULL);
  789. ASSERT(pV1 != NULL);
  790. ASSERT(pV2 != NULL);
  791. CVector3d bc(pV1,pV2);
  792. CVector3d ba(pV1,pV0);
  793. //bc.Trace();
  794. //ba.Trace();
  795. CVector3d v = ::Inner(&bc,&ba);
  796. //v.Trace();
  797. return v.GetNormL2()/2.0;
  798. }
  800. //**************************************
  801. // GetMeanArea
  802. //**************************************
  803. double GetMeanArea(CArray3d<CFace3d> *pArrayFace)
  804. {
  805. int size = pArrayFace->GetSize();
  806. ASSERT(size > 0);
  807. double sum = 0.0f;
  808. for(int i=0;i<size;i++)
  809. sum += pArrayFace->GetAt(i)->Area();
  810. return (sum/(double)size);
  811. }
  813. //**************************************
  814. // GetMinArea
  815. //**************************************
  816. double GetMinArea(CArray3d<CFace3d> *pArrayFace)
  817. {
  818. int size = pArrayFace->GetSize();
  819. ASSERT(size > 0);
  820. double min = MAX_DOUBLE;
  821. for(int i=0;i<size;i++)
  822. {
  823. double area = pArrayFace->GetAt(i)->Area();
  824. //TRACE("Area : %gn",area);
  825. min = (area < min) ? area : min;
  826. }
  827. return min;
  828. }
  831. //**************************************
  832. // FormFunction
  833. //**************************************
  834. // Ratio (area(013)/area(012))
  835. // Triangle : 012
  836. // Vertex on which we want to evaluate 
  837. // form function : 3
  838. //**************************************
  839. double FormFunction(CVertex3d *pV0,
  840. CVertex3d *pV1,
  841. CVertex3d *pV2,
  842. CVertex3d *pV3)
  843. {
  844. double area = Area(pV0,pV1,pV2);
  845. ASSERT(area != 0.0);
  846. return Area(pV0,pV1,pV3)/area;
  847. }
  850. //**************************************
  851. // IntersectionLineFace
  852. //**************************************
  853. int IntersectionLineFace(CVertex3d *pV0,
  854.  CVertex3d *pV1,
  855.  CFace3d *pFace,
  856.  CVertex3d *pVertexResult)
  857. {
  858. float a,b,c,d;
  859. PlaneEquationFromFace(pFace,&a,&b,&c,&d);
  860. float i = pV1->x()-pV0->x();
  861. float j = pV1->y()-pV0->y();
  862. float k = pV1->z()-pV0->z();
  864. float den = a*i + b*j + c*k;
  865. float num = -(a * pV0->x() + b * pV0->y() + c * pV0->z() + d);
  867. // parallel or face contain line
  868. if(den == 0.0f)
  869. return 0;
  871. float t = num / den;
  872. pVertexResult->Set(pV0->x()+t*i,pV0->y()+t*j,pV0->z()+t*k);
  874. return VertexInFace(pVertexResult,pFace);
  875. }
  877. //********************************************
  878. // IntersectionLineFaceRecursive
  879. // (propagation from pFaceStart)
  880. //********************************************
  881. int IntersectionLineFaceRecursive(CVertex3d *pV0,
  882. CVertex3d *pV1,
  883. CFace3d *pFaceStart,
  884. CVertex3d *pVertexResult,
  885. CFace3d **ppFaceResult,
  886. CArray3d<CFace3d> *pArrayFace,
  887. int MaxNbFace)
  888. {
  889. if(pArrayFace->GetSize() >= MaxNbFace)
  890. return 0;
  892. pArrayFace->Add(pFaceStart);
  893. if(IntersectionLineFace(pV0,pV1,pFaceStart,pVertexResult))
  894. {
  895. *ppFaceResult = pFaceStart;
  896. return 1;
  897. }
  899. //TRACE("Size : %d (max : %d)n",pArrayFace->GetSize(),MaxNbFace);
  901. pFaceStart->SetFlag(1);
  903. // neighboring
  904. for(int i=0;i<3;i++)
  905. {
  906. CFace3d *pFace = pFaceStart->f(i);
  907. if(pFace != NULL)  // neighbor exists
  908. if(pFace->GetFlag() != 1) // not yet
  909. if(IntersectionLineFaceRecursive(pV0,pV1,pFace,
  910.              pVertexResult,ppFaceResult,pArrayFace,MaxNbFace))
  911. return 1;
  912. }
  914. return 0;
  915. }
  917. //********************************************
  918. // NearestIntersectionWithLine
  919. // First intersection founded (propagation
  920. // from pFaceStart)
  921. //********************************************
  922. int NearestIntersectionWithLine(CVertex3d *pV0,
  923. CVertex3d *pV1,
  924. CFace3d *pFaceStart,
  925. CVertex3d *pVertexResult,
  926. CFace3d **ppFaceResult,
  927. int MaxNbFace)
  928. {
  929. // Check
  930. ASSERT(pV0 != NULL);
  931. ASSERT(pV1 != NULL);
  932. ASSERT(pFaceStart != NULL);
  933. ASSERT(pVertexResult != NULL);
  935. // Use flag
  936. CArray3d<CFace3d> ArrayFace;
  937. int success = IntersectionLineFaceRecursive(pV0,pV1,pFaceStart,
  938.                                           pVertexResult,ppFaceResult,&ArrayFace,MaxNbFace);
  939. for(int i=0;i<ArrayFace.GetSize();i++)
  940. ArrayFace[i]->SetFlag(0);
  942. //TRACE("NbFace : %dn",ArrayFace.GetSize());
  943. ArrayFace.SetSize(0);
  945. return success;
  946. }
  949. //********************************************
  950. // NearestIntersectionWithLine
  951. // Only search in pArrayFace
  952. //********************************************
  953. int NearestIntersectionWithLine(CArray3d<CFace3d> *pArrayFace,
  954. CVertex3d *pV0,
  955. CVertex3d *pV1,
  956. CVertex3d *pVertexResult,
  957. CFace3d **ppFaceResult,
  958. int *pNbVisitedFace)
  959. {
  960. int success = 0;
  961. *ppFaceResult = NULL;
  962. int size = pArrayFace->GetSize();
  963. *pNbVisitedFace = size;
  964. double SqDistance = MAX_DOUBLE;
  965. // For each face
  966. for(int i=0;i<size;i++)
  967. {
  968. CVertex3d v;
  969. if(IntersectionLineFace(pV0,pV1,pArrayFace->GetAt(i),&v))
  970. {
  971. success++;
  972. double d = ::DistanceSquare(pV0,&v);
  973. if(d < SqDistance)
  974. {
  975. SqDistance = d;
  976. pVertexResult->Set(v);
  977. // this face contains intersection
  978. *ppFaceResult = pArrayFace->GetAt(i);
  979. }
  980. }
  981. }
  982. //TRACE("NbIntersection : %dn",success);
  983. if(success)
  984. {
  985. ASSERT(ppFaceResult != NULL);
  986. }
  987. return success; // NbIntersection
  988. }
  991. //********************************************
  992. // NearestIntersectionWithLine
  993. // Only search in pFace and its neighbors
  994. //********************************************
  995. int NearestIntersectionWithLineFromFaceAndNeighbors(CFace3d *pFace,
  996. CVertex3d *pV0,
  997. CVertex3d *pV1,
  998. CVertex3d *pVertexResult,
  999. CFace3d **ppFaceResult)
  1000. {
  1001. int success = 0;
  1002. *ppFaceResult = NULL;
  1003. double SqDistance = MAX_DOUBLE;
  1004. CVertex3d v;
  1006. // pFace alone
  1007. if(IntersectionLineFace(pV0,pV1,pFace,&v))
  1008. {
  1009. success++;
  1010. double d = ::DistanceSquare(pV0,&v);
  1011. if(d < SqDistance)
  1012. {
  1013. SqDistance = d;
  1014. pVertexResult->Set(v);
  1015. *ppFaceResult = pFace;
  1016. }
  1017. }
  1019. // Its neighbors
  1020. for(int i=0;i<3;i++)
  1021. {
  1022. CFace3d *pFaceNeighbor = pFace->f(i);
  1023. if(pFaceNeighbor != NULL) // neighbor
  1024. if(IntersectionLineFace(pV0,pV1,pFaceNeighbor,&v))
  1025. {
  1026. success++;
  1027. double d = ::DistanceSquare(pV0,&v);
  1028. if(d < SqDistance)
  1029. {
  1030. SqDistance = d;
  1031. pVertexResult->Set(v);
  1032. *ppFaceResult = pFaceNeighbor;
  1033. }
  1034. }
  1035. }
  1037. if(success) { ASSERT((*ppFaceResult) != NULL); }
  1038. return success;
  1039. }
  1043. int GravityCenter(CArray3d<CVertex3d> *pArray,
  1044. CVertex3d *pVertexCenter)
  1045. {
  1046. int size = pArray->GetSize();
  1047. for(int j=0;j<3;j++)
  1048. {
  1049. float v = 0.0f;
  1050. for(int i=0;i<size;i++)
  1051. v += pArray->GetAt(i)->Get(j);
  1052. v /= size;
  1053. pVertexCenter->Set(j,v);
  1054. }
  1055. return 1;
  1056. }
  1058. // ** EOF **