ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkDiscretePrincipalCurvaturesQuadEdgeMeshFilter.h
Go to the documentation of this file.
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