Go to the documentation of this file.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 {
00081 p[i] = output->GetPoint( id[i] );
00082 }
00083
00084 if( !TriangleType::IsObtuse( p[0], p[1], p[2] ) )
00085 {
00086 OutputCurvatureType sq_d01 =
00087 static_cast< OutputCurvatureType >(
00088 p[0].SquaredEuclideanDistanceTo( p[1] ) );
00089 OutputCurvatureType sq_d02 =
00090 static_cast< OutputCurvatureType >(
00091 p[0].SquaredEuclideanDistanceTo( p[2] ) );
00092
00093 OutputCurvatureType cot_theta_210 =
00094 TriangleType::Cotangent( p[2], p[1], p[0] );
00095 OutputCurvatureType cot_theta_021 =
00096 TriangleType::Cotangent( p[0], p[2], p[1] );
00097
00098 return 0.125 * ( sq_d02 * cot_theta_210 + sq_d01 * cot_theta_021 );
00099 }
00100 else
00101 {
00102 OutputCurvatureType area =
00103 static_cast< OutputCurvatureType >(
00104 TriangleType::ComputeArea( p[0], p[1], p[2] ) );
00105 if( ( p[1] - p[0] ) * ( p[2] - p[0] ) < 0. )
00106 {
00107 return 0.5 * area;
00108 }
00109 else
00110 {
00111 return 0.25 * area;
00112 }
00113 }
00114 }
00115
00116 virtual void GenerateData()
00117 {
00118 this->CopyInputMeshToOutputMesh();
00119
00120 OutputMeshPointer output = this->GetOutput();
00121
00122 OutputPointsContainerPointer points = output->GetPoints();
00123 OutputPointsContainerIterator p_it = points->Begin();
00124
00125 OutputCurvatureType curvature;
00126
00127 while( p_it != points->End() )
00128 {
00129 curvature = EstimateCurvature( p_it->Value() );
00130 output->SetPointData( p_it->Index(), curvature );
00131 ++p_it;
00132 }
00133 }
00134
00135 private:
00136 QuadEdgeMeshDiscreteCurvatureEstimator( const Self& );
00137 void operator = ( const Self& );
00138
00139 };
00140
00141 }
00142
00143 #endif
00144