Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

itkTriangleHelper.h

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Program:   Insight Segmentation & Registration Toolkit
00004   Module:    $RCSfile: itkTriangleHelper.h,v $
00005   Language:  C++
00006   Date:      $Date: 2008-10-06 00:02:06 $
00007   Version:   $Revision: 1.6 $
00008 
00009   Copyright (c) Insight Software Consortium. All rights reserved.
00010   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
00011 
00012      This software is distributed WITHOUT ANY WARRANTY; without even
00013      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00014      PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 #ifndef __itkTriangleHelper_h
00018 #define __itkTriangleHelper_h
00019 
00020 #include "itkCrossHelper.h"
00021 
00022 namespace itk
00023 {
00029 template< typename TPoint >
00030 class TriangleHelper
00031 {
00032 public:
00033   typedef TriangleHelper                          Self;
00034   typedef TPoint                                  PointType;
00035   typedef typename PointType::CoordRepType        CoordRepType;
00036   typedef typename PointType::VectorType          VectorType;
00037   typedef CrossHelper< VectorType >               CrossVectorType;
00038 
00039   itkStaticConstMacro( PointDimension, unsigned int, PointType::PointDimension );
00040 
00041   static bool IsObtuse( const PointType& iA, const PointType& iB, const PointType& iC )
00042     {
00043     VectorType v01 = iB - iA;
00044     VectorType v02 = iC - iA;
00045     VectorType v12 = iC - iB;
00046 
00047     if( v01 * v02 < 0.0 )
00048       {
00049       return true;
00050       }
00051     else
00052       {
00053       if( v02 * v12 < 0.0 )
00054         {
00055         return true;
00056           }
00057       else
00058         {
00059         if( v01 * -v12 < 0.0 )
00060           {
00061           return true;
00062           }
00063         else
00064           {
00065           return false;
00066           }
00067         }
00068       }
00069     }
00070   
00071   static VectorType ComputeNormal ( const PointType& iA,
00072     const PointType& iB,
00073     const PointType& iC )
00074   {
00075     CrossVectorType cross;
00076     VectorType w = cross ( iB - iA, iC - iA );
00077     CoordRepType l2 = w.GetSquaredNorm();
00078 
00079     if( l2 != 0.0 )
00080       {
00081       w /= vcl_sqrt( l2 );
00082       }
00083 
00084     return w;
00085   }
00086 
00087   static CoordRepType Cotangent ( const PointType& iA,
00088                                   const PointType& iB,
00089                                   const PointType& iC )
00090     {
00091     VectorType v21 = iA - iB;
00092     v21.Normalize();
00093     
00094     VectorType v23 = iC - iB;
00095     v23.Normalize();
00096     
00097     CoordRepType bound( 0.999999 );
00098 
00099     CoordRepType cos_theta = vnl_math_max( -bound, 
00100       vnl_math_min( bound, v21 * v23 ) );
00101             
00102     return 1.0 / vcl_tan( vcl_acos( cos_theta ) );
00103     }
00104 
00105   static PointType ComputeBarycenter (
00106       const CoordRepType& iA1, const PointType& iP1,
00107       const CoordRepType& iA2, const PointType& iP2,
00108       const CoordRepType& iA3, const PointType& iP3 )
00109     {
00110     PointType oPt;
00111 
00112     for ( unsigned int dim = 0; dim < PointDimension; dim++ )
00113       {
00114       oPt[dim] = iA1 * iP1[dim] + iA2 * iP2[dim] + iA3 * iP3[dim];
00115       }
00116 
00117     return oPt;
00118     }
00119 
00120   static CoordRepType ComputeAngle( const PointType& iP1, const PointType& iP2,
00121       const PointType& iP3 )
00122     {
00123     VectorType v21 = iP1 - iP2;
00124     VectorType v23 = iP3 - iP2;
00125 
00126     CoordRepType v21_l2 = v21.GetSquaredNorm();
00127     CoordRepType v23_l2 = v23.GetSquaredNorm();
00128 
00129     if( v21_l2 != 0.0 )
00130       v21 /= v21_l2;
00131     if( v23_l2 != 0.0 )
00132       v23 /= v23_l2;
00133 
00134     CoordRepType bound( 0.999999 );
00135 
00136     CoordRepType cos_theta = vnl_math_max( -bound, 
00137       vnl_math_min( bound, v21 * v23 ) );
00138  
00139     return vcl_acos( cos_theta );
00140     }
00141 
00142   static PointType ComputeGravityCenter (
00143       const PointType& iP1,
00144       const PointType& iP2,
00145       const PointType& iP3 )
00146     {
00147     PointType oPt;
00148     CoordRepType inv_3 = 1.0 / 3.;
00149 
00150     for ( unsigned int dim = 0; dim < PointDimension; dim++ )
00151       {
00152       oPt[dim] = ( iP1[dim] + iP2[dim] + iP3[dim] ) * inv_3;
00153       }
00154 
00155     return oPt;
00156     }
00157 
00158   static PointType ComputeCircumCenter (
00159       const PointType& iP1,
00160       const PointType& iP2,
00161       const PointType& iP3 )
00162   {
00163     PointType oPt;
00164     oPt.Fill ( 0.0 );
00165 
00166     CoordRepType a = iP2.SquaredEuclideanDistanceTo ( iP3 );
00167     CoordRepType b = iP1.SquaredEuclideanDistanceTo ( iP3 );
00168     CoordRepType c = iP2.SquaredEuclideanDistanceTo ( iP1 );
00169 
00170     CoordRepType Weight[3];
00171     Weight[0] = a * ( b + c - a );
00172     Weight[1] = b * ( c + a - b );
00173     Weight[2] = c * ( a + b - c );
00174 
00175     CoordRepType SumWeight = Weight[0] + Weight[1] + Weight[2];
00176 
00177     if ( SumWeight != 0.0 )
00178       {
00179       SumWeight = 1.0 / SumWeight;
00180 
00181       for ( unsigned int dim = 0; dim < PointDimension; dim++ )
00182         {
00183         oPt[dim] = ( Weight[0] * iP1[dim] + Weight[1] * iP2[dim] + Weight[2] * iP3[dim] ) * SumWeight;
00184         }
00185       }
00186 
00187 
00188     return oPt;
00189   }
00190 
00191   static PointType ComputeConstrainedCircumCenter ( const PointType& iP1,
00192       const PointType& iP2, const PointType& iP3 )
00193     {
00194     PointType oPt;
00195     CoordRepType a = iP2.SquaredEuclideanDistanceTo ( iP3 );
00196     CoordRepType b = iP1.SquaredEuclideanDistanceTo ( iP3 );
00197     CoordRepType c = iP2.SquaredEuclideanDistanceTo ( iP1 );
00198 
00199     CoordRepType Weight[3];
00200     Weight[0] = a * ( b + c - a );
00201     Weight[1] = b * ( c + a - b );
00202     Weight[2] = c * ( a + b - c );
00203 
00204     for ( unsigned int i = 0; i < 3; i++ )
00205       {
00206       if ( Weight[i] < 0.0 )
00207         {
00208         Weight[i] = 0.;
00209         }
00210       }
00211 
00212     CoordRepType SumWeight = Weight[0] + Weight[1] + Weight[2];
00213 
00214     if ( SumWeight != 0.0 )
00215       {
00216       SumWeight = 1.0 / SumWeight;
00217 
00218       for ( unsigned int dim = 0; dim < PointDimension; dim++ )
00219         {
00220         oPt[dim] = ( Weight[0] * iP1[dim] + Weight[1] * iP2[dim] + Weight[2] * iP3[dim] ) * SumWeight;
00221         }
00222       }
00223 
00224     return oPt;
00225     }
00226 
00227   static CoordRepType ComputeArea ( const PointType& iP1, const PointType& iP2, const PointType& iP3 )
00228     {
00229     CoordRepType a = iP2.EuclideanDistanceTo ( iP3 );
00230     CoordRepType b = iP1.EuclideanDistanceTo ( iP3 );
00231     CoordRepType c = iP2.EuclideanDistanceTo ( iP1 );
00232 
00233     CoordRepType s = 0.5 * ( a + b + c );
00234     return static_cast< CoordRepType > ( vcl_sqrt ( s * ( s - a ) * ( s - b ) * ( s - c ) ) );
00235     }
00236 
00237 protected:
00238   TriangleHelper( );
00239   virtual ~TriangleHelper( );
00240 
00241   void PrintSelf ( std::ostream& os, Indent indent ) const;
00242 
00243 
00244 private:
00245   TriangleHelper( const Self& ); // purposely not implemented
00246   void operator = ( const Self& ); // purposely not implemented
00247 };
00248 }
00249 
00250 #endif
00251 

Generated at Sat Feb 28 13:41:09 2009 for ITK by doxygen 1.5.6 written by Dimitri van Heesch, © 1997-2000