vcglib: triangle3.h Source File

00001 /****************************************************************************
00002 * VCGLib                                                            o o     *
00003 * Visual and Computer Graphics Library                            o     o   *
00004 *                                                                _   O  _   *
00005 * Copyright(C) 2004                                                \/)\/    *
00006 * Visual Computing Lab                                            /\/|      *
00007 * ISTI - Italian National Research Council                           |      *
00008 *                                                                    \      *
00009 * All rights reserved.                                                      *
00010 *                                                                           *
00011 * This program is free software; you can redistribute it and/or modify      *   
00012 * it under the terms of the GNU General Public License as published by      *
00013 * the Free Software Foundation; either version 2 of the License, or         *
00014 * (at your option) any later version.                                       *
00015 *                                                                           *
00016 * This program is distributed in the hope that it will be useful,           *
00017 * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
00018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
00019 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
00020 * for more details.                                                         *
00021 *                                                                           *
00022 ****************************************************************************/
00023 /****************************************************************************
00024   History
00025 
00026 $Log: not supported by cvs2svn $
00027 Revision 1.21  2007/12/02 07:39:19  cignoni
00028 disambiguated sqrt call
00029 
00030 Revision 1.20  2007/11/26 14:11:38  ponchio
00031 Added Mean Ratio metric for triangle quality.
00032 
00033 Revision 1.19  2007/11/19 17:04:05  ponchio
00034 QualityRadii values fixed.
00035 
00036 Revision 1.18  2007/11/18 19:12:54  ponchio
00037 Typo (missing comma).
00038 
00039 Revision 1.17  2007/11/16 14:22:35  ponchio
00040 Added qualityRadii: computes inradius /circumradius.
00041 (ok the name is ugly...)
00042 
00043 Revision 1.16  2007/10/10 15:11:30  ponchio
00044 Added Circumcenter function.
00045 
00046 Revision 1.15  2007/05/10 09:31:15  cignoni
00047 Corrected InterpolationParameters invocation
00048 
00049 Revision 1.14  2007/05/04 16:33:27  ganovelli
00050 moved InterpolationParamaters out the class Triangle
00051 
00052 Revision 1.13  2007/04/04 23:23:55  pietroni
00053 - corrected and renamed distance to point ( function TrianglePointDistance)
00054 
00055 Revision 1.12  2007/01/13 00:25:23  cignoni
00056 Added (Normalized) Normal version templated on three points (instead forcing the creation of a new triangle)
00057 
00058 Revision 1.11  2006/10/17 06:51:33  fiorin
00059 In function Barycenter, replaced calls to (the inexistent) cP(i) with P(i)
00060 
00061 Revision 1.10  2006/10/10 09:33:47  cignoni
00062 added quality for triangle wrap
00063 
00064 Revision 1.9  2006/09/14 08:44:07  ganovelli
00065 changed t.P(*) in t.cP() nella funzione Barycenter
00066 
00067 Revision 1.8  2006/06/01 08:38:58  pietroni
00068 added PointDistance function
00069 
00070 Revision 1.7  2006/03/01 15:35:09  pietroni
00071 compiled InterspolationParameters function
00072 
00073 Revision 1.6  2006/01/22 10:00:56  cignoni
00074 Very Important Change: Area->DoubleArea (and no more Area function)
00075 
00076 Revision 1.5  2005/09/23 14:18:27  ganovelli
00077 added constructor
00078 
00079 Revision 1.4  2005/04/14 11:35:09  ponchio
00080 *** empty log message ***
00081 
00082 Revision 1.3  2004/07/15 13:22:37  cignoni
00083 Added the standard P() access function instead of the shortcut P0()
00084 
00085 Revision 1.2  2004/07/15 10:17:42  pietroni
00086 correct access to point funtions call in usage of triangle3 (ex. t.P(0)  in t.P0(0))
00087 
00088 Revision 1.1  2004/03/08 01:13:31  cignoni
00089 Initial commit
00090 
00091 
00092 ****************************************************************************/
00093 #ifndef __VCG_TRIANGLE3
00094 #define __VCG_TRIANGLE3
00095 
00096 #include <vcg/space/box3.h>
00097 #include <vcg/space/point2.h>
00098 #include <vcg/space/point3.h>
00099 #include <vcg/space/plane3.h>
00100 #include <vcg/space/segment3.h>
00101 
00102 namespace vcg {
00103 
00110 template <class ScalarTriangleType> class Triangle3
00111 {
00112 public:
00113   typedef ScalarTriangleType ScalarType;
00114         typedef Point3< ScalarType > CoordType;
00116         typedef Box3<ScalarType> BoxType;
00117 
00118 /*********************************************
00119     blah
00120     blah
00121 **/
00122         Triangle3(){}
00123         Triangle3(const CoordType & c0,const CoordType & c1,const CoordType & c2){_v[0]=c0;_v[1]=c1;_v[2]=c2;}
00124 protected:
00126         Point3<ScalarType> _v[3];
00127 public:
00128 
00130         inline CoordType & P( const int j ) { return _v[j];}
00131         inline CoordType & P0( const int j ) { return _v[j];}
00132         inline CoordType & P1( const int j ) { return _v[(j+1)%3];}
00133         inline CoordType & P2( const int j ) { return _v[(j+2)%3];}
00134         inline const CoordType &  P( const int j ) const { return _v[j];}
00135         inline const CoordType &  P0( const int j ) const { return _v[j];}
00136         inline const CoordType &  P1( const int j ) const { return _v[(j+1)%3];}
00137         inline const CoordType &  P2( const int j ) const { return _v[(j+2)%3];}
00138         inline const CoordType & cP0( const int j ) const { return _v[j];}
00139         inline const CoordType & cP1( const int j ) const { return _v[(j+1)%3];}
00140         inline const CoordType & cP2( const int j ) const { return _v[(j+2)%3];}
00141 
00142 
00143 
00144         bool InterpolationParameters(const CoordType & bq, ScalarType &a, ScalarType &b, ScalarType &_c ) const{
00145                 return InterpolationParameters(*this, bq, a, b,_c ); 
00146         }
00147 
00148 
00149 
00150 
00151 
00153 ScalarType QualityFace( ) const
00154 {
00155         return Quality(P(0), P(1), P(2));
00156 }
00157 
00158 }; //end Class
00159 
00161 template<class TriangleType>
00162 typename TriangleType::ScalarType QualityFace(const TriangleType &t)
00163 {
00164   return Quality(t.cP(0), t.cP(1), t.cP(2));
00165 }
00166 
00167 // More robust function to computing barycentric coords of a point inside a triangle.
00168 // it requires the knowledge of what is the direction that is more orthogonal to the face plane. 
00169 // Usually this info can be stored in a bit of the face flags (see updateFlags::FaceProjection(MeshType &m) )
00170 // and accessing the field with
00171 //        if(fp->Flags() & FaceType::NORMX )   axis = 0;
00172 //        else if(fp->Flags() & FaceType::NORMY )   axis = 1;
00173 //        else axis =2;
00174 //        InterpolationParameters(*fp,axis,Point,Bary);
00175 // This direction is used to project the triangle in 2D and solve the problem in 2D where it is well defined.
00176 
00177 template<class TriangleType, class ScalarType>
00178 bool InterpolationParameters(const TriangleType t, const int Axis, const Point3<ScalarType> & P,  Point3<ScalarType> & L)
00179 {
00180   Point2<ScalarType> test;
00181         typedef Point2<ScalarType> P2;
00182         if(Axis==0) return InterpolationParameters2( P2(t.P(0)[1],t.P(0)[2]), P2(t.P(1)[1],t.P(1)[2]), P2(t.P(2)[1],t.P(2)[2]), P2(P[1],P[2]), L);
00183         if(Axis==1) return InterpolationParameters2( P2(t.P(0)[0],t.P(0)[2]), P2(t.P(1)[0],t.P(1)[2]), P2(t.P(2)[0],t.P(2)[2]), P2(P[0],P[2]), L);
00184         if(Axis==2) return InterpolationParameters2( P2(t.P(0)[0],t.P(0)[1]), P2(t.P(1)[0],t.P(1)[1]), P2(t.P(2)[0],t.P(2)[1]), P2(P[0],P[1]), L);
00185         return false; 
00186 }
00188 template<class TriangleType, class ScalarType>
00189 bool InterpolationParameters(const TriangleType t, const Point3<ScalarType> & N, const Point3<ScalarType> & P,  Point3<ScalarType> & L)
00190 {
00191         if(N[0]>N[1])
00192         {
00193                 if(N[0]>N[2])
00194                         return InterpolationParameters(t,0,P,L); /* 0 > 1 ? 2 */
00195                 else 
00196                         return InterpolationParameters(t,2,P,L); /* 2 > 1 ? 2 */
00197                 }
00198         else 
00199         {
00200                 if(N[1]>N[2])
00201                         return InterpolationParameters(t,1,P,L); /* 1 > 0 ? 2 */
00202                 else 
00203                         return InterpolationParameters(t,2,P,L); /* 2 > 1 ? 2 */
00204         }
00205 }
00206                 
00207 
00208 
00209 // Function that computes the barycentric coords of a 2D triangle. Used by the above function. 
00210 // Algorithm: simply find a base for the frame of the triangle, assuming v3 as origin (matrix T) invert it and apply to P-v3.
00211 
00212 template<class ScalarType>
00213 bool InterpolationParameters2(const Point2<ScalarType> &V1,
00214                                                                                                            const Point2<ScalarType> &V2,
00215                                                                                                                  const Point2<ScalarType> &V3,
00216                                                                                                                  const Point2<ScalarType> &P, Point3<ScalarType> &L)
00217 {
00218         ScalarType T00 = V1[0]-V3[0];  ScalarType T01 = V2[0]-V3[0];  
00219         ScalarType T10 = V1[1]-V3[1];  ScalarType T11 = V2[1]-V3[1]; 
00220         ScalarType Det = T00 * T11 - T01*T10;
00221         if(fabs(Det) < 0.0000001) 
00222                                 return false;
00223         
00224         ScalarType IT00 =  T11/Det;     ScalarType IT01 = -T01/Det;     
00225         ScalarType IT10 = -T10/Det;     ScalarType IT11 =  T00/Det;     
00226 
00227         Point2<ScalarType> Delta = P-V3;
00228         
00229         L[0] = IT00*Delta[0] + IT01*Delta[1];
00230         L[1] = IT10*Delta[0] + IT11*Delta[1]; 
00231  
00232         if(L[0]<0) L[0]=0;
00233         if(L[1]<0) L[1]=0;
00234         if(L[0]>1.) L[0]=1;
00235         if(L[1]>1.) L[1]=1;
00236         
00237         L[2] = 1. - L[1] - L[0];
00238         if(L[2]<0) L[2]=0;
00239         
00240         assert(L[2] >= -0.00001);
00241 
00242         return true;
00243 }
00244 
00245 
00253 template<class TriangleType, class ScalarType>
00254 bool InterpolationParameters(const TriangleType t,const Point3<ScalarType> & N,const Point3<ScalarType> & bq, ScalarType &a, ScalarType &b, ScalarType &c ) 
00255 {       
00256         Point3<ScalarType> bary;
00257         bool done= InterpolationParameters(t,N,bq,bary);
00258         a=bary[0];
00259         b=bary[1];
00260         c=bary[2];
00261         return done;
00262 //        else if(fp->Flags() & FaceType::NORMY )   axis = 1;
00263 //        else axis =2;
00264 //        InterpolationParameters(*fp,axis,Point,Bary);
00265 //const ScalarType _EPSILON = ScalarType(0.000001);
00266 //#define x1 (t.P(0).X())
00267 //#define y1 (t.P(0).Y())
00268 //#define z1 (t.P(0).Z())
00269 //#define x2 (t.P(1).X())
00270 //#define y2 (t.P(1).Y())
00271 //#define z2 (t.P(1).Z())
00272 //#define x3 (t.P(2).X())
00273 //#define y3 (t.P(2).Y())
00274 //#define z3 (t.P(2).Z())
00275 //#define px (bq[0])
00276 //#define py (bq[1])
00277 //#define pz (bq[2])
00278 //
00279 //     ScalarType t1  = px*y2;
00280 //     ScalarType t2  = px*y3;
00281 //     ScalarType t3  = py*x2;
00282 //     ScalarType t4  = py*x3;
00283 //     ScalarType t5  = x2*y3;
00284 //     ScalarType t6  = x3*y2;
00285 //     ScalarType t8  = x1*y2;
00286 //     ScalarType t9  = x1*y3;
00287 //     ScalarType t10 = y1*x2;
00288 //     ScalarType t11 = y1*x3;
00289 //     ScalarType t13 = t8-t9-t10+t11+t5-t6;
00290 //     if(fabs(t13)>=_EPSILON)
00291 //       {
00292 //         ScalarType t15 = px*y1;
00293 //         ScalarType t16 = py*x1;
00294 //         a =  (t1 -t2-t3 +t4+t5-t6 )/t13;
00295 //         b = -(t15-t2-t16+t4+t9-t11)/t13;
00296 //         _c =  (t15-t1-t16+t3+t8-t10)/t13;
00297 //              return true;
00298 //     }
00299 //
00300 //     t1  = px*z2;
00301 //     t2  = px*z3;
00302 //     t3  = pz*x2;
00303 //     t4  = pz*x3;
00304 //     t5  = x2*z3;
00305 //     t6  = x3*z2;
00306 //     t8  = x1*z2;
00307 //     t9  = x1*z3;
00308 //     t10 = z1*x2;
00309 //     t11 = z1*x3;
00310 //     t13 = t8-t9-t10+t11+t5-t6;
00311 //     if(fabs(t13)>=_EPSILON)
00312 //       {
00313 //              ScalarType t15 = px*z1;
00314 //              ScalarType t16 = pz*x1;
00315 //              a =  (t1 -t2-t3 +t4+t5-t6 )/t13;
00316 //              b = -(t15-t2-t16+t4+t9-t11)/t13;
00317 //              _c =  (t15-t1-t16+t3+t8-t10)/t13;
00318 //              return true;
00319 //     }
00320 //
00321 //     t1  = pz*y2; t2  = pz*y3;
00322 //     t3  = py*z2; t4  = py*z3;
00323 //     t5  = z2*y3; t6  = z3*y2;
00324 //     t8  = z1*y2; t9  = z1*y3;
00325 //     t10 = y1*z2; t11 = y1*z3;
00326 //     t13 = t8-t9-t10+t11+t5-t6;
00327 //     if(fabs(t13)>=_EPSILON)
00328 //       {
00329 //         ScalarType t15 = pz*y1;
00330 //         ScalarType t16 = py*z1;
00331 //         a =  (t1 -t2-t3 +t4+t5-t6 )/t13;
00332 //         b = -(t15-t2-t16+t4+t9-t11)/t13;
00333 //         _c =  (t15-t1-t16+t3+t8-t10)/t13;
00334 //              return true;
00335 //     }
00336 //       
00337 //#undef x1
00338 //#undef y1
00339 //#undef z1
00340 //#undef x2
00341 //#undef y2
00342 //#undef z2
00343 //#undef x3
00344 //#undef y3
00345 //#undef z3
00346 //#undef px
00347 //#undef py
00348 //#undef pz
00349 //
00350 //     return false;
00351 }
00352 
00353 
00358 template<class P3ScalarType>
00359 P3ScalarType Quality( Point3<P3ScalarType> const &p0, Point3<P3ScalarType> const & p1,  Point3<P3ScalarType> const & p2)
00360 {
00361         Point3<P3ScalarType> d10=p1-p0;
00362         Point3<P3ScalarType> d20=p2-p0;
00363         Point3<P3ScalarType> d12=p1-p2;
00364         Point3<P3ScalarType> x = d10^d20;
00365 
00366         P3ScalarType a = Norm( x );
00367         if(a==0) return 0; // Area zero triangles have surely quality==0;
00368         P3ScalarType b = SquaredNorm( d10 );
00369         P3ScalarType t = b;
00370         t = SquaredNorm( d20 ); if ( b<t ) b = t;
00371         t = SquaredNorm( d12 ); if ( b<t ) b = t;
00372         assert(b!=0.0);
00373         return a/b;
00374 }
00375 
00376 
00381 template<class P3ScalarType>
00382 P3ScalarType QualityRadii(Point3<P3ScalarType> const &p0,
00383                                                                                                         Point3<P3ScalarType> const &p1,
00384                                                                                                         Point3<P3ScalarType> const &p2) {
00385 
00386         P3ScalarType a=(p1-p0).Norm();
00387         P3ScalarType b=(p2-p0).Norm();
00388         P3ScalarType c=(p1-p2).Norm();
00389 
00390         P3ScalarType sum = (a + b + c)*0.5;
00391         P3ScalarType area2 =  sum*(a+b-sum)*(a+c-sum)*(b+c-sum);
00392         if(area2 <= 0) return 0;
00393         //circumradius: (a*b*c)/(4*sqrt(area2))
00394         //inradius: (a*b*c)/(4*circumradius*sum) => sqrt(area2)/sum;
00395         return (8*area2)/(a*b*c*sum);
00396 }
00397 
00402 template<class P3ScalarType>
00403 P3ScalarType QualityMeanRatio(Point3<P3ScalarType> const &p0,
00404                                                                                                         Point3<P3ScalarType> const &p1,
00405                                                                                                         Point3<P3ScalarType> const &p2) {
00406 
00407         P3ScalarType a=(p1-p0).Norm();
00408         P3ScalarType b=(p2-p0).Norm();
00409         P3ScalarType c=(p1-p2).Norm();
00410         P3ScalarType sum = (a + b + c)*0.5; //semiperimeter
00411         P3ScalarType area2 =  sum*(a+b-sum)*(a+c-sum)*(b+c-sum);
00412         if(area2 <= 0) return 0;
00413         return (4.0*sqrt(3.0)*sqrt(area2))/(a*a + b*b + c*c);
00414 }
00415 
00417 template<class TriangleType>
00418 Point3<typename TriangleType::ScalarType> Normal(const TriangleType &t)
00419 {
00420         return (( t.P(1) - t.P(0)) ^ (t.P(2) - t.P(0)));
00421 }
00422 template<class Point3Type>
00423 Point3Type Normal( Point3Type const &p0, Point3Type const & p1,  Point3Type const & p2)
00424 {
00425         return (( p1 - p0) ^ (p2 - p0));
00426 }
00427 
00428 
00430 template<class TriangleType>
00431 Point3<typename TriangleType::ScalarType> NormalizedNormal(const TriangleType &t)
00432 {
00433         return (( t.P(1) - t.P(0)) ^ (t.P(2) - t.P(0))).Normalize();
00434 }
00435 template<class Point3Type>
00436 Point3Type NormalizedNormal( Point3Type const &p0, Point3Type const & p1,  Point3Type const & p2)
00437 {
00438         return (( p1 - p0) ^ (p2 - p0)).Normalize();
00439 }
00440 
00441 
00442 
00444 // NOTE the old Area function has been removed to intentionally 
00445 // cause compiling error that will help people to check their code...
00446 // A some  people used Area assuming that it returns the double and some not. 
00447 // So please check your codes!!!
00448 // And please DO NOT Insert any Area named function here!
00449 
00450 template<class TriangleType>
00451 typename TriangleType::ScalarType DoubleArea(const TriangleType &t) 
00452 {
00453         return Norm( (t.P(1) - t.P(0)) ^ (t.P(2) - t.P(0)) );
00454 }
00455 
00456 template<class TriangleType>
00457 typename TriangleType::ScalarType CosWedge(const TriangleType &t, int k)
00458 {
00459   typename TriangleType::CoordType 
00460     e0 = t.P((k+1)%3) - t.P(k),
00461     e1 = t.P((k+2)%3) - t.P(k);
00462   return (e0*e1)/(e0.Norm()*e1.Norm());
00463 }
00464 
00465 template<class TriangleType>
00466 Point3<typename TriangleType::ScalarType> Barycenter(const TriangleType &t) 
00467 {
00468         return ((t.P(0)+t.P(1)+t.P(2))/(typename TriangleType::ScalarType) 3.0);
00469 }
00470 
00471 template<class TriangleType>
00472 typename TriangleType::ScalarType Perimeter(const TriangleType &t) 
00473 {
00474         return Distance(t.P(0),t.P(1))+
00475                      Distance(t.P(1),t.P(2))+
00476                                  Distance(t.P(2),t.P(0));
00477 }
00478 
00479 template<class TriangleType>
00480 Point3<typename TriangleType::ScalarType> Circumcenter(const TriangleType &t)
00481 {
00482    typename TriangleType::ScalarType a2 = (t.P(1) - t.P(2)).SquaredNorm();
00483    typename TriangleType::ScalarType b2 = (t.P(2) - t.P(0)).SquaredNorm();
00484    typename TriangleType::ScalarType c2 = (t.P(0) - t.P(1)).SquaredNorm();      
00485    Point3<typename TriangleType::ScalarType>c = t.P(0)*a2*(-a2 + b2 + c2) + 
00486                                                 t.P(1)*b2*( a2 - b2 + c2) + 
00487                                                 t.P(2)*c2*( a2 + b2 - c2);
00488    c /= 2*(a2*b2 + a2*c2 + b2*c2) - a2*a2 - b2*b2 - c2*c2;
00489    return c;
00490 }
00491 
00501 template<class TriangleType>
00502 void TrianglePointDistance(const  TriangleType &t,
00503                                                         const typename TriangleType::CoordType & q,
00504                                                         typename TriangleType::ScalarType & dist, 
00505                                                         typename TriangleType::CoordType & closest )
00506 {
00507         typedef typename TriangleType::CoordType CoordType;
00508         typedef typename TriangleType::ScalarType ScalarType;
00509 
00510         CoordType clos[3];
00511         ScalarType distv[3];
00512         CoordType clos_proj;
00513         ScalarType distproj;
00514 
00516         vcg::Plane3<ScalarType> plane;
00517         plane.Init(t.P(0),t.P(1),t.P(2));
00518         clos_proj=plane.Projection(q);
00519 
00521         CoordType n=(t.P(1)-t.P(0))^(t.P(2)-t.P(0));
00522         CoordType n0=(t.P(0)-clos_proj)^(t.P(1)-clos_proj);
00523         CoordType n1=(t.P(1)-clos_proj)^(t.P(2)-clos_proj);
00524         CoordType n2=(t.P(2)-clos_proj)^(t.P(0)-clos_proj);
00525         distproj=(clos_proj-q).Norm();
00526         if (((n*n0)>=0)&&((n*n1)>=0)&&((n*n2)>=0))
00527         {
00528                 closest=clos_proj;
00529                 dist=distproj;
00530                 return;
00531         }
00532         
00533 
00534         //distance from the edges
00535         vcg::Segment3<ScalarType> e0=vcg::Segment3<ScalarType>(t.P(0),t.P(1));
00536         vcg::Segment3<ScalarType> e1=vcg::Segment3<ScalarType>(t.P(1),t.P(2));
00537         vcg::Segment3<ScalarType> e2=vcg::Segment3<ScalarType>(t.P(2),t.P(0));
00538         clos[0]=ClosestPoint<ScalarType>( e0, q);
00539         clos[1]=ClosestPoint<ScalarType>( e1, q);
00540         clos[2]=ClosestPoint<ScalarType>( e2, q);
00541         
00542         distv[0]=(clos[0]-q).Norm();
00543         distv[1]=(clos[1]-q).Norm();
00544         distv[2]=(clos[2]-q).Norm();
00545         int min=0;
00546 
00548         for (int i=1;i<3;i++)
00549         {
00550                 if (distv[i]<distv[min])
00551                         min=i;
00552         }
00553 
00554         closest=clos[min];
00555         dist=distv[min];
00556 }
00557 
00558 
00559 }        // end namespace
00560 
00561 
00562 #endif
00563