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 __itkDiscreteCurvatureQuadEdgeMeshFilter_h 00019 #define __itkDiscreteCurvatureQuadEdgeMeshFilter_h 00020 00021 #include "itkQuadEdgeMeshToQuadEdgeMeshFilter.h" 00022 #include "itkConceptChecking.h" 00023 #include "itkTriangleHelper.h" 00024 00025 namespace itk 00026 { 00034 template< class TInputMesh, class TOutputMesh=TInputMesh > 00035 class ITK_EXPORT DiscreteCurvatureQuadEdgeMeshFilter: 00036 public QuadEdgeMeshToQuadEdgeMeshFilter< TInputMesh, TOutputMesh > 00037 { 00038 public: 00039 typedef DiscreteCurvatureQuadEdgeMeshFilter Self; 00040 typedef SmartPointer< Self > Pointer; 00041 typedef SmartPointer< const Self > ConstPointer; 00042 typedef QuadEdgeMeshToQuadEdgeMeshFilter< TInputMesh, TOutputMesh > Superclass; 00043 00044 typedef TInputMesh InputMeshType; 00045 typedef typename InputMeshType::Pointer InputMeshPointer; 00046 00047 typedef TOutputMesh OutputMeshType; 00048 typedef typename OutputMeshType::Pointer OutputMeshPointer; 00049 typedef typename OutputMeshType::PointsContainerPointer OutputPointsContainerPointer; 00050 typedef typename OutputMeshType::PointsContainerIterator OutputPointsContainerIterator; 00051 typedef typename OutputMeshType::PointType OutputPointType; 00052 typedef typename OutputPointType::CoordRepType OutputCoordType; 00053 typedef typename OutputMeshType::PointIdentifier OutputPointIdentifier; 00054 typedef typename OutputMeshType::CellIdentifier OutputCellIdentifier; 00055 typedef typename OutputMeshType::QEType OutputQEType; 00056 typedef typename OutputMeshType::MeshTraits OutputMeshTraits; 00057 typedef typename OutputMeshTraits::PixelType OutputCurvatureType; 00058 00059 typedef TriangleHelper< OutputPointType > TriangleType; 00060 00062 itkTypeMacro(DiscreteCurvatureQuadEdgeMeshFilter, QuadEdgeMeshToQuadEdgeMeshFilter); 00063 00064 #ifdef ITK_USE_CONCEPT_CHECKING 00065 00066 itkConceptMacro( OutputIsFloatingPointCheck, 00067 ( Concept::IsFloatingPoint< OutputCurvatureType > ) ); 00068 00070 #endif 00071 00072 protected: 00073 DiscreteCurvatureQuadEdgeMeshFilter() {} 00074 ~DiscreteCurvatureQuadEdgeMeshFilter() {} 00075 00076 virtual OutputCurvatureType EstimateCurvature(const OutputPointType & iP) = 0; 00077 00078 OutputCurvatureType ComputeMixedArea(OutputQEType *iQE1, OutputQEType *iQE2) 00079 { 00080 OutputMeshPointer output = this->GetOutput(); 00081 00082 OutputPointIdentifier id[3]; 00083 00084 id[0] = iQE1->GetOrigin(); 00085 id[1] = iQE1->GetDestination(); 00086 id[2] = iQE2->GetDestination(); 00087 00088 OutputPointType p[3]; 00089 00090 for ( int i = 0; i < 3; i++ ) 00091 { 00092 p[i] = output->GetPoint(id[i]); 00093 } 00094 00095 if ( !TriangleType::IsObtuse(p[0], p[1], p[2]) ) 00096 { 00097 OutputCurvatureType sq_d01 = 00098 static_cast< OutputCurvatureType >( 00099 p[0].SquaredEuclideanDistanceTo(p[1]) ); 00100 OutputCurvatureType sq_d02 = 00101 static_cast< OutputCurvatureType >( 00102 p[0].SquaredEuclideanDistanceTo(p[2]) ); 00103 00104 OutputCurvatureType cot_theta_210 = 00105 TriangleType::Cotangent(p[2], p[1], p[0]); 00106 OutputCurvatureType cot_theta_021 = 00107 TriangleType::Cotangent(p[0], p[2], p[1]); 00108 00109 return 0.125 * ( sq_d02 * cot_theta_210 + sq_d01 * cot_theta_021 ); 00110 } 00111 else 00112 { 00113 OutputCurvatureType area = 00114 static_cast< OutputCurvatureType >( 00115 TriangleType::ComputeArea(p[0], p[1], p[2]) ); 00116 if ( ( p[1] - p[0] ) * ( p[2] - p[0] ) < 0. ) 00117 { 00118 return 0.5 * area; 00119 } 00120 else 00121 { 00122 return 0.25 * area; 00123 } 00124 } 00125 } 00126 00127 virtual void GenerateData() 00128 { 00129 this->CopyInputMeshToOutputMesh(); 00130 00131 OutputMeshPointer output = this->GetOutput(); 00132 00133 OutputPointsContainerPointer points = output->GetPoints(); 00134 OutputPointsContainerIterator p_it = points->Begin(); 00135 00136 OutputCurvatureType curvature; 00137 00138 while ( p_it != points->End() ) 00139 { 00140 curvature = EstimateCurvature( p_it->Value() ); 00141 output->SetPointData(p_it->Index(), curvature); 00142 ++p_it; 00143 } 00144 } 00145 00146 private: 00147 DiscreteCurvatureQuadEdgeMeshFilter(const Self &); // purposely not 00148 // implemented 00149 void operator=(const Self &); // purposely not 00150 // implemented 00151 }; 00152 } // end namespace itk 00153 00154 #endif 00155