ITK  4.2.0
Insight Segmentation and Registration Toolkit
itkQuadEdgeMeshDecimationQuadricElementHelper.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef __itkQuadEdgeMeshDecimationQuadricElementHelper_h
19 #define __itkQuadEdgeMeshDecimationQuadricElementHelper_h
20 
21 #include "itkPoint.h"
22 #include "vnl/vnl_vector_fixed.h"
23 #include "vnl/vnl_matrix.h"
24 #include "vnl/algo/vnl_matrix_inverse.h"
25 
26 #include "itkTriangleHelper.h"
27 
28 namespace itk
29 {
31 template< class TPoint >
33 {
34 public:
35 
37  typedef TPoint PointType;
38  typedef typename PointType::CoordRepType CoordType;
39 
40  itkStaticConstMacro(PointDimension, unsigned int, PointType::PointDimension);
41  itkStaticConstMacro(NumberOfCoefficients, unsigned int,
42  PointDimension * ( PointDimension + 1 ) / 2 + PointDimension + 1);
43 
44  typedef typename PointType::VectorType VectorType;
46  typedef vnl_vector_fixed< CoordType,
47  itkGetStaticConstMacro(PointDimension) > VNLVectorType;
48  typedef vnl_vector_fixed< CoordType,
49  itkGetStaticConstMacro(NumberOfCoefficients) > CoefficientVectorType;
51 
52  // *****************************************************************
54  m_Coefficients(itk::NumericTraits< CoordType >::Zero),
55  m_A(PointDimension, PointDimension, itk::NumericTraits< CoordType >::Zero),
56  m_B(itk::NumericTraits< CoordType >::Zero),
57  m_SVDAbsoluteThreshold( static_cast< CoordType >( 1e-6 ) ),
58  m_SVDRelativeThreshold( static_cast< CoordType >( 1e-3 ) )
59  {
60  this->m_Rank = PointDimension;
61  }
62 
64  m_Coefficients(iCoefficients),
65  m_A(PointDimension, PointDimension, itk::NumericTraits< CoordType >::Zero),
66  m_B(itk::NumericTraits< CoordType >::Zero),
67  m_SVDAbsoluteThreshold( static_cast< CoordType >( 1e-3 ) ),
68  m_SVDRelativeThreshold( static_cast< CoordType >( 1e-3 ) )
69  {
70  this->m_Rank = PointDimension;
71  this->ComputeAMatrixAndBVector();
72  }
73 
75  {}
76 
77  CoefficientVectorType GetCoefficients() const
78  {
79  return this->m_Coefficients;
80  }
81 
82  VNLMatrixType GetAMatrix()
83  {
84  this->ComputeAMatrixAndBVector();
85  return m_A;
86  }
87 
88  VNLVectorType GetBVector()
89  {
90  ComputeAMatrixAndBVector();
91  return m_B;
92  }
93 
94  unsigned int GetRank() const
95  {
96  return m_Rank;
97  }
98 
100  inline CoordType ComputeError(const PointType & iP) const
101  {
102  // ComputeAMatrixAndBVector();
103  vnl_svd< CoordType > svd(m_A, m_SVDAbsoluteThreshold);
104  svd.zero_out_relative(m_SVDRelativeThreshold);
105  CoordType oError = inner_product( iP.GetVnlVector(), svd.recompose() * iP.GetVnlVector() );
106 
107  return this->m_Coefficients[this->m_Coefficients.size() - 1] - oError;
108  /*
109  CoordType oError( 0. );
110 
111  std::vector< CoordType > pt( PointDimension + 1, 1. );
112 
113  unsigned int dim1( 0 ), dim2, k( 0 );
114 
115  while( dim1 < PointDimension )
116  {
117  pt[dim1] = iP[dim1];
118  ++dim1;
119  }
120 
121  for( dim1 = 0; dim1 < PointDimension + 1; ++dim1 )
122  {
123  oError += this->m_Coefficients[k++] * pt[dim1] * pt[dim1];
124 
125  for( dim2 = dim1 + 1; dim2 < PointDimension + 1; ++dim2 )
126  {
127  oError += 2. * this->m_Coefficients[k++] * pt[dim1] * pt[dim2];
128  }
129  }
130  oError += this->m_Coefficients[k++];
131 
132  return oError;*/
133  }
134 
136  inline CoordType ComputeErrorAtOptimalLocation(const PointType & iP)
137  {
138  PointType optimal_location = ComputeOptimalLocation(iP);
139 
140  return ComputeError(optimal_location);
141  }
142 
143  PointType ComputeOptimalLocation(const PointType & iP)
144  {
145  ComputeAMatrixAndBVector();
146 
147  vnl_svd< CoordType > svd(m_A, m_SVDAbsoluteThreshold);
148  svd.zero_out_relative(m_SVDRelativeThreshold);
149 
150  m_Rank = svd.rank();
151 
152  VNLVectorType y = m_B.as_vector() - m_A *iP.GetVnlVector();
153 
154  VNLVectorType displacement = svd.solve(y);
155  PointType oP;
156 
157  for ( unsigned int dim = 0; dim < PointDimension; dim++ )
158  {
159  oP[dim] = iP[dim] + displacement[dim];
160  }
161 
162  return oP;
163  }
164 
166  PointType ComputeOptimalLocation(
167  const unsigned int & )
168  {}
169 
170  void AddTriangle( const PointType & iP1,
171  const PointType & iP2,
172  const PointType & iP3,
173  const CoordType & iWeight = static_cast< CoordType >( 1. ) )
174  {
175  VectorType N = TriangleType::ComputeNormal(iP1, iP2, iP3);
176 
177  AddPoint(iP1, N, iWeight);
178  }
179 
180  void AddPoint( const PointType & iP,
181  const VectorType & iN,
182  const CoordType & iWeight = static_cast< CoordType >( 1. ) )
183  {
184  unsigned int k(0), dim1, dim2;
185 
186  CoordType d = -iN *iP.GetVectorFromOrigin();
187 
188  for ( dim1 = 0; dim1 < PointDimension; ++dim1 )
189  {
190  for ( dim2 = dim1; dim2 < PointDimension; ++dim2 )
191  {
192  this->m_Coefficients[k++] += iWeight * iN[dim1] * iN[dim2];
193  }
194  this->m_Coefficients[k++] += iWeight * iN[dim1] * d;
195  }
196 
197  this->m_Coefficients[k++] += iWeight * d * d;
198  }
199 
200  // ***********************************************************************
201  // operators
202  Self & operator=(const Self & iRight)
203  {
204  this->m_Coefficients = iRight.m_Coefficients;
205  return *this;
206  }
207 
208  Self operator+(const Self & iRight) const
209  {
210  return Self(this->m_Coefficients + iRight.m_Coefficients);
211  }
212 
213  Self & operator+=(const Self & iRight)
214  {
215  this->m_Coefficients += iRight.m_Coefficients;
216  return *this;
217  }
218 
219  Self operator-(const Self & iRight) const
220  {
221  return Self(this->m_Coefficients - iRight.m_Coefficients);
222  }
223 
224  Self & operator-=(const Self & iRight)
225  {
226  this->m_Coefficients -= iRight.m_Coefficients;
227  return *this;
228  }
229 
230  Self operator*(const CoordType & iV) const
231  {
232  Self oElement = Self(this->m_Coefficients * iV);
233 
234  return oElement;
235  }
236 
237  Self & operator*=(const CoordType & iV)
238  {
239  this->m_Coefficients *= iV;
240  return *this;
241  }
242 
243 protected:
244 
248  unsigned int m_Rank;
251  //bool m_MatrixFilled;
252 
253  void ComputeAMatrixAndBVector()
254  {
255  unsigned int k(0), dim1, dim2;
256 
257  for ( dim1 = 0; dim1 < PointDimension; ++dim1 )
258  {
259  for ( dim2 = dim1; dim2 < PointDimension; ++dim2 )
260  {
261  m_A[dim1][dim2] = m_A[dim2][dim1] = m_Coefficients[k++];
262  }
263  m_B[dim1] = -m_Coefficients[k++];
264  }
265  //m_MatrixFilled = true;
266  }
267 };
268 }
269 #endif
270