ITK
4.0.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 __itkQuadEdgeMeshDecimationQuadricElementHelper_h 00019 #define __itkQuadEdgeMeshDecimationQuadricElementHelper_h 00020 00021 #include "itkPoint.h" 00022 #include "vnl/vnl_vector_fixed.h" 00023 #include "vnl/vnl_matrix.h" 00024 #include "vnl/algo/vnl_matrix_inverse.h" 00025 00026 #include "itkTriangleHelper.h" 00027 00028 namespace itk 00029 { 00031 template< class TPoint > 00032 class ITK_EXPORT QuadEdgeMeshDecimationQuadricElementHelper 00033 { 00034 public: 00035 00036 typedef QuadEdgeMeshDecimationQuadricElementHelper Self; 00037 typedef TPoint PointType; 00038 typedef typename PointType::CoordRepType CoordType; 00039 00040 itkStaticConstMacro(PointDimension, unsigned int, PointType::PointDimension); 00041 itkStaticConstMacro(NumberOfCoefficients, unsigned int, 00042 PointDimension * ( PointDimension + 1 ) / 2 + PointDimension + 1); 00043 00044 typedef typename PointType::VectorType VectorType; 00045 typedef vnl_matrix< CoordType > VNLMatrixType; 00046 typedef vnl_vector_fixed< CoordType, 00047 itkGetStaticConstMacro(PointDimension) > VNLVectorType; 00048 typedef vnl_vector_fixed< CoordType, 00049 itkGetStaticConstMacro(NumberOfCoefficients) > CoefficientVectorType; 00050 typedef TriangleHelper< PointType > TriangleType; 00051 00052 // ***************************************************************** 00053 QuadEdgeMeshDecimationQuadricElementHelper(): 00054 m_Coefficients(itk::NumericTraits< CoordType >::Zero), 00055 m_A(PointDimension, PointDimension, itk::NumericTraits< CoordType >::Zero), 00056 m_B(itk::NumericTraits< CoordType >::Zero), 00057 m_SVDAbsoluteThreshold( static_cast< CoordType >( 1e-6 ) ), 00058 m_SVDRelativeThreshold( static_cast< CoordType >( 1e-3 ) ) 00059 { 00060 this->m_Rank = PointDimension; 00061 } 00062 00063 QuadEdgeMeshDecimationQuadricElementHelper(const CoefficientVectorType & iCoefficients): 00064 m_Coefficients(iCoefficients), 00065 m_A(PointDimension, PointDimension, itk::NumericTraits< CoordType >::Zero), 00066 m_B(itk::NumericTraits< CoordType >::Zero), 00067 m_SVDAbsoluteThreshold( static_cast< CoordType >( 1e-3 ) ), 00068 m_SVDRelativeThreshold( static_cast< CoordType >( 1e-3 ) ) 00069 { 00070 this->m_Rank = PointDimension; 00071 this->ComputeAMatrixAndBVector(); 00072 } 00073 00074 ~QuadEdgeMeshDecimationQuadricElementHelper() 00075 {} 00076 00077 CoefficientVectorType GetCoefficients() const 00078 { 00079 return this->m_Coefficients; 00080 } 00081 00082 VNLMatrixType GetAMatrix() 00083 { 00084 this->ComputeAMatrixAndBVector(); 00085 return m_A; 00086 } 00087 00088 VNLVectorType GetBVector() 00089 { 00090 ComputeAMatrixAndBVector(); 00091 return m_B; 00092 } 00093 00094 unsigned int GetRank() const 00095 { 00096 return m_Rank; 00097 } 00098 00100 inline CoordType ComputeError(const PointType & iP) const 00101 { 00102 // ComputeAMatrixAndBVector(); 00103 vnl_svd< CoordType > svd(m_A, m_SVDAbsoluteThreshold); 00104 svd.zero_out_relative(m_SVDRelativeThreshold); 00105 CoordType oError = inner_product( iP.GetVnlVector(), svd.recompose() * iP.GetVnlVector() ); 00106 00107 return this->m_Coefficients[this->m_Coefficients.size() - 1] - oError; 00108 /* 00109 CoordType oError( 0. ); 00110 00111 std::vector< CoordType > pt( PointDimension + 1, 1. ); 00112 00113 unsigned int dim1( 0 ), dim2, k( 0 ); 00114 00115 while( dim1 < PointDimension ) 00116 { 00117 pt[dim1] = iP[dim1]; 00118 ++dim1; 00119 } 00120 00121 for( dim1 = 0; dim1 < PointDimension + 1; ++dim1 ) 00122 { 00123 oError += this->m_Coefficients[k++] * pt[dim1] * pt[dim1]; 00124 00125 for( dim2 = dim1 + 1; dim2 < PointDimension + 1; ++dim2 ) 00126 { 00127 oError += 2. * this->m_Coefficients[k++] * pt[dim1] * pt[dim2]; 00128 } 00129 } 00130 oError += this->m_Coefficients[k++]; 00131 00132 return oError;*/ 00133 } 00134 00136 inline CoordType ComputeErrorAtOptimalLocation(const PointType & iP) 00137 { 00138 PointType optimal_location = ComputeOptimalLocation(iP); 00139 00140 return ComputeError(optimal_location); 00141 } 00142 00143 PointType ComputeOptimalLocation(const PointType & iP) 00144 { 00145 ComputeAMatrixAndBVector(); 00146 00147 vnl_svd< CoordType > svd(m_A, m_SVDAbsoluteThreshold); 00148 svd.zero_out_relative(m_SVDRelativeThreshold); 00149 00150 m_Rank = svd.rank(); 00151 00152 VNLVectorType y = m_B.as_vector() - m_A *iP.GetVnlVector(); 00153 00154 VNLVectorType displacement = svd.solve(y); 00155 PointType oP; 00156 00157 for ( unsigned int dim = 0; dim < PointDimension; dim++ ) 00158 { 00159 oP[dim] = iP[dim] + displacement[dim]; 00160 } 00161 00162 return oP; 00163 } 00164 00166 PointType ComputeOptimalLocation( 00167 const unsigned int & iNumberOfEigenValues) 00168 {} 00169 00170 // PointType ComputeOptimalLocation( 00171 // const CoordType& iValue ) 00172 // { 00173 // ComputeAMatrixAndBVector(); 00174 // vnl_svd< CoordType > svd( m_A ); 00175 // svd.zero.zero_out_relative( iValue ); 00176 // m_Rank = svd.rank(); 00177 // 00178 // VNLVectorType location = svd.solve( m_B ); 00179 // PointType oP; 00180 // 00181 // for( unsigned int dim = 0; dim < PointDimension; dim++ ) 00182 // oP[dim] = location[dim]; 00183 // 00184 // return oP; 00185 // } 00186 00187 void AddTriangle( const PointType & iP1, 00188 const PointType & iP2, 00189 const PointType & iP3, 00190 const CoordType & iWeight = static_cast< CoordType >( 1. ) ) 00191 { 00192 VectorType N = TriangleType::ComputeNormal(iP1, iP2, iP3); 00193 00194 AddPoint(iP1, N, iWeight); 00195 } 00196 00197 void AddPoint( const PointType & iP, 00198 const VectorType & iN, 00199 const CoordType & iWeight = static_cast< CoordType >( 1. ) ) 00200 { 00201 unsigned int k(0), dim1, dim2; 00202 00203 CoordType d = -iN *iP.GetVectorFromOrigin(); 00204 00205 for ( dim1 = 0; dim1 < PointDimension; ++dim1 ) 00206 { 00207 for ( dim2 = dim1; dim2 < PointDimension; ++dim2 ) 00208 { 00209 this->m_Coefficients[k++] += iWeight * iN[dim1] * iN[dim2]; 00210 } 00211 this->m_Coefficients[k++] += iWeight * iN[dim1] * d; 00212 } 00213 00214 this->m_Coefficients[k++] += iWeight * d * d; 00215 } 00216 00217 // *********************************************************************** 00218 // operators 00219 Self & operator=(const Self & iRight) 00220 { 00221 this->m_Coefficients = iRight.m_Coefficients; 00222 return *this; 00223 } 00224 00225 Self operator+(const Self & iRight) const 00226 { 00227 return Self(this->m_Coefficients + iRight.m_Coefficients); 00228 } 00229 00230 Self & operator+=(const Self & iRight) 00231 { 00232 this->m_Coefficients += iRight.m_Coefficients; 00233 return *this; 00234 } 00235 00236 Self operator-(const Self & iRight) const 00237 { 00238 return Self(this->m_Coefficients - iRight.m_Coefficients); 00239 } 00240 00241 Self & operator-=(const Self & iRight) 00242 { 00243 this->m_Coefficients -= iRight.m_Coefficients; 00244 return *this; 00245 } 00246 00247 Self operator*(const CoordType & iV) const 00248 { 00249 Self oElement = Self(this->m_Coefficients * iV); 00250 00251 return oElement; 00252 } 00253 00254 Self & operator*=(const CoordType & iV) 00255 { 00256 this->m_Coefficients *= iV; 00257 return *this; 00258 } 00259 00260 protected: 00261 00262 CoefficientVectorType m_Coefficients; 00263 VNLMatrixType m_A; 00264 VNLVectorType m_B; 00265 unsigned int m_Rank; 00266 CoordType m_SVDAbsoluteThreshold; 00267 CoordType m_SVDRelativeThreshold; 00268 //bool m_MatrixFilled; 00269 00270 void ComputeAMatrixAndBVector() 00271 { 00272 unsigned int k(0), dim1, dim2; 00273 00274 for ( dim1 = 0; dim1 < PointDimension; ++dim1 ) 00275 { 00276 for ( dim2 = dim1; dim2 < PointDimension; ++dim2 ) 00277 { 00278 m_A[dim1][dim2] = m_A[dim2][dim1] = m_Coefficients[k++]; 00279 } 00280 m_B[dim1] = -m_Coefficients[k++]; 00281 } 00282 //m_MatrixFilled = true; 00283 } 00284 }; 00285 } 00286 #endif 00287