ITK
4.1.0
Insight Segmentation and Registration Toolkit
|
00001 /*========================================================================= 00002 * 00003 * Copyright Insight Software Consortium 00004 * 00005 * Licensed under the Apache License, Version 2.0 (the "License"); 00006 * you may not use this file except in compliance with the License. 00007 * You may obtain a copy of the License at 00008 * 00009 * http://www.apache.org/licenses/LICENSE-2.0.txt 00010 * 00011 * Unless required by applicable law or agreed to in writing, software 00012 * distributed under the License is distributed on an "AS IS" BASIS, 00013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 00014 * See the License for the specific language governing permissions and 00015 * limitations under the License. 00016 * 00017 *=========================================================================*/ 00018 #ifndef __itkDiscretePrincipalCurvaturesQuadEdgeMeshFilter_h 00019 #define __itkDiscretePrincipalCurvaturesQuadEdgeMeshFilter_h 00020 00021 #include "itkDiscreteCurvatureQuadEdgeMeshFilter.h" 00022 #include "itkQuadEdgeMeshParamMatrixCoefficients.h" 00023 00024 namespace itk 00025 { 00033 template< class TInputMesh, class TOutputMesh=TInputMesh > 00034 class ITK_EXPORT DiscretePrincipalCurvaturesQuadEdgeMeshFilter: 00035 public DiscreteCurvatureQuadEdgeMeshFilter< TInputMesh, TOutputMesh > 00036 { 00037 public: 00038 typedef DiscretePrincipalCurvaturesQuadEdgeMeshFilter Self; 00039 typedef SmartPointer< Self > Pointer; 00040 typedef SmartPointer< const Self > ConstPointer; 00041 typedef DiscreteCurvatureQuadEdgeMeshFilter< 00042 TInputMesh, TOutputMesh > Superclass; 00043 00044 typedef typename Superclass::InputMeshType InputMeshType; 00045 typedef typename Superclass::InputMeshPointer InputMeshPointer; 00046 typedef typename Superclass::OutputMeshType OutputMeshType; 00047 typedef typename Superclass::OutputMeshPointer OutputMeshPointer; 00048 typedef typename Superclass::OutputPointsContainerPointer OutputPointsContainerPointer; 00049 typedef typename Superclass::OutputPointsContainerIterator OutputPointsContainerIterator; 00050 typedef typename Superclass::OutputPointType OutputPointType; 00051 typedef typename Superclass::OutputVectorType OutputVectorType; 00052 typedef typename Superclass::OutputCoordType OutputCoordType; 00053 typedef typename Superclass::OutputPointIdentifier OutputPointIdentifier; 00054 typedef typename Superclass::OutputCellIdentifier OutputCellIdentifier; 00055 typedef typename Superclass::OutputQEType OutputQEType; 00056 typedef typename Superclass::OutputMeshTraits OutputMeshTraits; 00057 typedef typename Superclass::OutputCurvatureType OutputCurvatureType; 00058 00059 typedef typename Superclass::TriangleType TriangleType; 00060 00062 itkTypeMacro(DiscretePrincipalCurvaturesQuadEdgeMeshFilter, DiscreteCurvatureQuadEdgeMeshFilter); 00063 00064 typedef ConformalMatrixCoefficients< OutputMeshType > CoefficientType; 00065 00066 #ifdef ITK_USE_CONCEPT_CHECKING 00067 00068 itkConceptMacro( OutputIsFloatingPointCheck, 00069 ( Concept::IsFloatingPoint< OutputCurvatureType > ) ); 00070 00072 #endif 00073 00074 protected: 00075 DiscretePrincipalCurvaturesQuadEdgeMeshFilter(): 00076 m_Gaussian(0.0), m_Mean(0.0){} 00077 ~DiscretePrincipalCurvaturesQuadEdgeMeshFilter() {} 00078 00079 OutputCurvatureType m_Gaussian; 00080 OutputCurvatureType m_Mean; 00081 00082 void ComputeMeanAndGaussianCurvatures(const OutputPointType & iP) 00083 { 00084 OutputMeshPointer output = this->GetOutput(); 00085 00086 OutputQEType *qe = iP.GetEdge(); 00087 00088 m_Mean = 0.; 00089 m_Gaussian = 0.; 00090 00091 if ( qe != 0 ) 00092 { 00093 OutputVectorType Laplace; 00094 Laplace.Fill(0.); 00095 00096 OutputQEType *qe_it = qe; 00097 00098 OutputCurvatureType area(0.), sum_theta(0.); 00099 00100 if ( qe_it != qe_it->GetOnext() ) 00101 { 00102 qe_it = qe; 00103 OutputQEType *qe_it2; 00104 00105 OutputPointType q0, q1; 00106 OutputVectorType face_normal; 00107 00108 OutputVectorType normal; 00109 normal.Fill(0.); 00110 00111 OutputCurvatureType temp_area; 00112 OutputCoordType temp_coeff; 00113 00114 CoefficientType coefficent; 00115 00116 do 00117 { 00118 qe_it2 = qe_it->GetOnext(); 00119 q0 = output->GetPoint( qe_it->GetDestination() ); 00120 q1 = output->GetPoint( qe_it2->GetDestination() ); 00121 00122 temp_coeff = coefficent(output, qe_it); 00123 Laplace += temp_coeff * ( iP - q0 ); 00124 00125 // Compute Angle; 00126 sum_theta += static_cast< OutputCurvatureType >( 00127 TriangleType::ComputeAngle(q0, iP, q1) ); 00128 00129 temp_area = this->ComputeMixedArea(qe_it, qe_it2); 00130 area += temp_area; 00131 00132 face_normal = TriangleType::ComputeNormal(q0, iP, q1); 00133 normal += face_normal; 00134 00135 qe_it = qe_it2; 00136 } 00137 while ( qe_it != qe ); 00138 00139 if ( area > 1e-10 ) 00140 { 00141 area = 1. / area; 00142 Laplace *= 0.25 * area; 00143 m_Mean = Laplace * normal; 00144 m_Gaussian = ( 2. * vnl_math::pi - sum_theta ) * area; 00145 } 00146 } 00147 } 00148 } 00149 00150 virtual OutputCurvatureType ComputeDelta() 00151 { 00152 return vnl_math_max( static_cast<OutputCurvatureType>( 0. ), 00153 m_Mean * m_Mean - m_Gaussian ); 00154 } 00155 00156 private: 00157 DiscretePrincipalCurvaturesQuadEdgeMeshFilter(const Self &); // purposely 00158 // not 00159 // implemented 00160 void operator=(const Self &); // purposely 00161 // not 00162 // implemented 00163 }; 00164 } 00165 00166 #endif 00167