00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef __itkQuadEdgeMeshDiscreteCurvatureEstimator_h
00019 #define __itkQuadEdgeMeshDiscreteCurvatureEstimator_h
00020
00021 #include <itkQuadEdgeMeshToQuadEdgeMeshFilter.h>
00022 #include "itkTriangleHelper.h"
00023
00024 namespace itk
00025 {
00032 template< class TInputMesh, class TOutputMesh >
00033 class QuadEdgeMeshDiscreteCurvatureEstimator :
00034 public QuadEdgeMeshToQuadEdgeMeshFilter< TInputMesh, TOutputMesh >
00035 {
00036 public:
00037 typedef QuadEdgeMeshDiscreteCurvatureEstimator Self;
00038 typedef SmartPointer< Self > Pointer;
00039 typedef SmartPointer< const Self > ConstPointer;
00040 typedef QuadEdgeMeshToQuadEdgeMeshFilter< TInputMesh, TOutputMesh > Superclass;
00041
00042 typedef TInputMesh InputMeshType;
00043 typedef typename InputMeshType::Pointer InputMeshPointer;
00044
00045 typedef TOutputMesh OutputMeshType;
00046 typedef typename OutputMeshType::Pointer OutputMeshPointer;
00047 typedef typename OutputMeshType::PointsContainerPointer OutputPointsContainerPointer;
00048 typedef typename OutputMeshType::PointsContainerIterator OutputPointsContainerIterator;
00049 typedef typename OutputMeshType::PointType OutputPointType;
00050 typedef typename OutputPointType::CoordRepType OutputCoordType;
00051 typedef typename OutputMeshType::PointIdentifier OutputPointIdentifier;
00052 typedef typename OutputMeshType::CellIdentifier OutputCellIdentifier;
00053 typedef typename OutputMeshType::QEType OutputQEType;
00054 typedef typename OutputMeshType::MeshTraits OutputMeshTraits;
00055 typedef typename OutputMeshTraits::PixelType OutputCurvatureType;
00056
00057 typedef TriangleHelper< OutputPointType > TriangleType;
00058
00060 itkTypeMacro( QuadEdgeMeshDiscreteCurvatureEstimator, QuadEdgeMeshToQuadEdgeMeshFilter );
00061
00062 protected:
00063 QuadEdgeMeshDiscreteCurvatureEstimator() {}
00064 ~QuadEdgeMeshDiscreteCurvatureEstimator() {}
00065
00066 virtual OutputCurvatureType EstimateCurvature( const OutputPointType& iP ) = 0;
00067
00068 OutputCurvatureType ComputeMixedArea( OutputQEType* iQE1, OutputQEType* iQE2 )
00069 {
00070 OutputMeshPointer output = this->GetOutput();
00071
00072 OutputPointIdentifier id[3];
00073 id[0] = iQE1->GetOrigin();
00074 id[1] = iQE1->GetDestination();
00075 id[2] = iQE2->GetDestination();
00076
00077 OutputPointType p[3];
00078
00079 for( int i = 0; i < 3; i++ )
00080 p[i] = output->GetPoint( id[i] );
00081
00082 if( !TriangleType::IsObtuse( p[0], p[1], p[2] ) )
00083 {
00084 OutputCurvatureType sq_d01 =
00085 static_cast< OutputCurvatureType >(
00086 p[0].SquaredEuclideanDistanceTo( p[1] ) );
00087 OutputCurvatureType sq_d02 =
00088 static_cast< OutputCurvatureType >(
00089 p[0].SquaredEuclideanDistanceTo( p[2] ) );
00090
00091 OutputCurvatureType cot_theta_210 =
00092 TriangleType::Cotangent( p[2], p[1], p[0] );
00093 OutputCurvatureType cot_theta_021 =
00094 TriangleType::Cotangent( p[0], p[2], p[1] );
00095
00096 return 0.125 * ( sq_d02 * cot_theta_210 + sq_d01 * cot_theta_021 );
00097 }
00098 else
00099 {
00100 OutputCurvatureType area =
00101 static_cast< OutputCurvatureType >(
00102 TriangleType::ComputeArea( p[0], p[1], p[2] ) );
00103 if( ( p[1] - p[0] ) * ( p[2] - p[0] ) < 0. )
00104 return 0.5 * area;
00105 else
00106 return 0.25 * area;
00107 }
00108 }
00109
00110 virtual void GenerateData()
00111 {
00112 Superclass::GenerateData();
00113 OutputMeshPointer output = this->GetOutput();
00114
00115 OutputPointsContainerPointer points = output->GetPoints();
00116
00117 for( OutputPointsContainerIterator p_it = points->Begin();
00118 p_it != points->End();
00119 p_it++ )
00120 {
00121 output->SetPointData( p_it->Index(), EstimateCurvature( p_it->Value() ) );
00122 }
00123 }
00124
00125 private:
00126 QuadEdgeMeshDiscreteCurvatureEstimator( const Self& );
00127 void operator = ( const Self& );
00128 };
00129
00130 }
00131
00132 #endif
00133