Medical Imaging Interaction Toolkit  2018.4.99-389bf124
Medical Imaging Interaction Toolkit
mitkItkNonUniformBSpline.txx
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 /*===================================================================
14 
15 This file is based heavily on a corresponding ITK filter.
16 
17 ===================================================================*/
18 
19 /*********************************
20  This file was taken from ITK, CVS version 1.13 to circumvent a bug in ITK release 3.18 (see http://public.kitware.com/Bug/view.php?id=10633
21  *********************************/
22 
23 #ifndef __itkNonUniformBSpline_txx
24 #define __itkNonUniformBSpline_txx
25 
26 #if defined(_MSC_VER)
27 #pragma warning ( disable : 4786 )
28 #endif
29 
30 #include "mitkItkNonUniformBSpline.h"
31 
32 #include "vnl/vnl_vector.h"
33 #include "vnl/vnl_matrix.h"
34 #include "vnl/algo/vnl_lsqr.h"
35 #include "vnl/vnl_linear_system.h"
36 
37 
38 // #define DEBUG_SPLINE
39 
40 namespace itk
41 {
42 
43 /** Constructor */
44 template< unsigned int TDimension >
45 NonUniformBSpline< TDimension >
46 ::NonUniformBSpline()
47 {
48  // Cubic bspline => 4th order
49  m_SplineOrder = 3;
50  m_SpatialDimension = TDimension;
51 }
52 
53 /** Destructor */
54 template< unsigned int TDimension >
55 NonUniformBSpline< TDimension >
56 ::~NonUniformBSpline()
57 {
58 }
59 
60 /** Print the object */
61 template< unsigned int TDimension >
62 void
63 NonUniformBSpline< TDimension >
64 ::PrintSelf( std::ostream& os, Indent indent ) const
65 {
66  Superclass::PrintSelf( os, indent );
67  os << indent << "NonUniformBSpline(" << this << ")" << std::endl;
68 
69  os << indent << "Chord lengths : " << std::endl;
70  for (ChordLengthListType::const_iterator iter = m_CumulativeChordLength.begin();
71  iter != m_CumulativeChordLength.end();
72  iter++)
73  {
74  os << indent << indent << *iter << std::endl;
75  }
76  os << indent << "Knots : " << std::endl;
77  for (KnotListType::const_iterator kiter = m_Knots.begin();
78  kiter != m_Knots.end();
79  kiter++)
80  {
81  os << indent << indent << *kiter << std::endl;
82  }
83  os << indent << "Control Points : " << std::endl;
84  for (typename ControlPointListType::const_iterator cpiter = m_ControlPoints.begin();
85  cpiter != m_ControlPoints.end();
86  cpiter++)
87  {
88  os << indent << indent << *cpiter << std::endl;
89  }
90 }
91 
92 
93 /** Set the list of points composing the tube */
94 template< unsigned int TDimension >
95 void
96 NonUniformBSpline< TDimension >
97 ::SetPoints( PointListType & points )
98 {
99  m_Points.clear();
100 
101  typename PointListType::iterator it,end;
102  it = points.begin();
103  end = points.end();
104  while(it != end)
105  {
106  m_Points.push_back(*it);
107  it++;
108  }
109 
110  this->Modified();
111 }
112 
113 /** Set the list of points composing the tube */
114 template< unsigned int TDimension >
115 void
116 NonUniformBSpline< TDimension >
117 ::SetKnots( KnotListType & knots )
118 {
119  m_Knots.clear();
120 
121  int len = knots.size();
122  double max_knot = knots[len - 1];
123 
124  typename KnotListType::iterator it;
125  typename KnotListType::iterator end;
126 
127  it = knots.begin();
128  end = knots.end();
129 
130  while(it != end)
131  {
132  m_Knots.push_back(*it/max_knot);
133  it++;
134  }
135 
136  this->Modified();
137 }
138 
139 template< unsigned int TDimension >
140 double
141 NonUniformBSpline< TDimension >
142 ::NonUniformBSplineFunctionRecursive(unsigned int order, unsigned int i, double t) const
143 {
144  if (order == 1)
145  {
146  if (m_Knots[i] <= t && t < m_Knots[i+1])
147  {
148  return 1;
149  }
150  else
151  {
152  return 0;
153  }
154  }
155 
156  //
157  // Be careful, we must use the passed in parameter for the order since this
158  // function is recursive.
159  //
160  double numer1 = (t - m_Knots[i]) * NonUniformBSplineFunctionRecursive(order-1, i, t);
161  double denom1 = (m_Knots[i+order-1] - m_Knots[i]);
162  double val1 = numer1 / denom1;
163  if (denom1 == 0 && numer1 == 0)
164  val1 = 0;
165  else if (denom1 == 0)
166  std::cout << "Error : " << denom1 << ", " << numer1 << std::endl;
167 
168  double numer2 = (m_Knots[i+order] - t) * NonUniformBSplineFunctionRecursive(order-1, i+1, t);
169  double denom2 = (m_Knots[i + order] - m_Knots[i+1]);
170  double val2 = numer2 / denom2;
171  if (denom2 == 0 && numer2 == 0)
172  val2 = 0;
173  else if (denom2 == 0)
174  std::cout << "Error : " << denom2 << ", " << numer2 << std::endl;
175 
176  return val1 + val2;
177 }
178 
179 template< unsigned int TDimension >
180 void
181 NonUniformBSpline< TDimension >
182 ::ComputeChordLengths()
183 {
184  m_ChordLength.clear();
185  m_CumulativeChordLength.clear();
186 
187  m_ChordLength.push_back(0);
188  m_CumulativeChordLength.push_back(0);
189 
190  double total_chord_length = 0.0;
191  ChordLengthListType temp;
192 
193  for (::size_t i = 0; i < m_Points.size()-1; i++)
194  {
195  PointType pt = m_Points[i];
196  PointType pt2 = m_Points[i+1];
197 
198  double chord = pt.EuclideanDistanceTo(pt2);
199  m_ChordLength.push_back(chord);
200  total_chord_length = total_chord_length + chord;
201  temp.push_back(total_chord_length);
202  }
203 
204  for (ChordLengthListType::iterator aiter = temp.begin();
205  aiter != temp.end();
206  aiter++)
207  {
208  m_CumulativeChordLength.push_back(*aiter/total_chord_length);
209  }
210 
211  //
212  // Debug printouts
213  //
214 #ifdef DEBUG_SPLINE
215  std::cout << "Total chord length : " << total_chord_length << std::endl;
216 
217  std::cout << "Chord length : " << std::endl;
218  for (ChordLengthListType::iterator aiter2 = m_ChordLength.begin();
219  aiter2 != m_ChordLength.end();
220  aiter2++)
221  {
222  std::cout << *aiter2 << std::endl;
223  }
224 
225  std::cout << "Cumulative chord length : " << std::endl;
226  for (ChordLengthListType::iterator aiter3 = m_CumulativeChordLength.begin();
227  aiter3 != m_CumulativeChordLength.end();
228  aiter3++)
229  {
230  std::cout << *aiter3 << std::endl;
231  }
232  std::cout << std::endl;
233 #endif
234 }
235 
236 template< unsigned int TDimension >
237 void
238 NonUniformBSpline< TDimension >
239 ::SetControlPoints( ControlPointListType& ctrlpts )
240 {
241  m_ControlPoints.clear();
242  for (typename ControlPointListType::iterator iter = ctrlpts.begin();
243  iter != ctrlpts.end();
244  iter++)
245  {
246  m_ControlPoints.push_back(*iter);
247  }
248  this->Modified();
249 }
250 
251 
252 template< unsigned int TDimension >
253 const typename
254 NonUniformBSpline< TDimension >::ControlPointListType &
255 NonUniformBSpline< TDimension >::GetControlPoints() const
256 {
257  return this->m_ControlPoints;
258 }
259 
260 
261 template< unsigned int TDimension >
262 const typename
263 NonUniformBSpline< TDimension >::KnotListType &
264 NonUniformBSpline< TDimension >::GetKnots() const
265 {
266  return this->m_Knots;
267 }
268 
269 
270 template< unsigned int TDimension >
271 const typename
272 NonUniformBSpline< TDimension >::PointListType &
273 NonUniformBSpline< TDimension >::GetPoints() const
274 {
275  return this->m_Points;
276 }
277 
278 
279 template< unsigned int TDimension >
280 void
281 NonUniformBSpline< TDimension >::ComputeControlPoints()
282 {
283  unsigned int dim = m_Points[0].GetPointDimension();
284 
285 #ifdef DEBUG_SPLINE
286  std::cout << "Points have dimension : " << dim << std::endl;
287 #endif
288 
289  //
290  // +1 in cols for radius
291  //
292  vnl_matrix<double> data_matrix(m_Points.size(), dim);
293 
294  //
295  // Form data point matrix
296  //
297  int rr = 0;
298  for (typename PointListType::iterator iter = m_Points.begin();
299  iter != m_Points.end();
300  iter++)
301  {
302  PointType pt = (*iter);
303  for (unsigned int i = 0; i < dim; i++)
304  {
305  data_matrix(rr, i) = pt.GetVnlVector()[i];
306  }
307  rr++;
308  }
309 
310 #ifdef DEBUG_SPLINE
311  std::cout << std::endl << "Data matrix" << std::endl;
312  std::cout << data_matrix << std::endl;
313 #endif
314 
315  //
316  // Form basis function matrix
317  //
318  //int num_basis_functions = 2 * m_SplineOrder - 1;
319  //int num_basis_functions = m_Points.size();
320  int num_rows = m_Points.size();
321 
322  //
323  // Assumes multiplicity k (m_SplineOrder at the ends).
324  //
325  int num_cols = m_Knots.size() - m_SplineOrder;
326 
327  vnl_matrix<double> N_matrix(num_rows, num_cols);
328 
329  //N_matrix(0, 0) = 1.0;
330 
331  for (int r = 0; r < num_rows; r++)
332  {
333  for (int c = 0; c < num_cols; c++)
334  {
335  double t = m_CumulativeChordLength[r];
336  N_matrix(r, c) = NonUniformBSplineFunctionRecursive(m_SplineOrder, c, t);
337  }
338  }
339 
340  N_matrix(num_rows-1, num_cols-1) = 1.0;
341 
342 #ifdef DEBUG_SPLINE
343  std::cout << "Basis function matrix : " << std::endl;
344  std::cout << N_matrix << std::endl;
345 #endif
346 
347 //FIXME: Use the LSQR linear solver here:
348  vnl_matrix<double> B;
349 
350 // = vnl_matrix_inverse<double>(N_matrix.transpose() * N_matrix) * N_matrix.transpose() * data_matrix;
351 
352 // vnl_linear_system ls( N_matrix.rows(), N_matrix.cols() );
353 
354 // vnl_lsqr solver( ls );
355 
356 
357 //#ifdef DEBUG_SPLINE
358  std::cout << "Control point matrix : " << std::endl;
359  std::cout << B << std::endl;
360 //#endif
361 
362  m_ControlPoints.clear();
363 
364  for ( unsigned int j = 0; j < B.rows(); j++ )
365  {
366  vnl_vector<double> v = B.get_row(j);
367  itk::Vector<double> iv;
368  iv.SetVnlVector(v);
369  itk::Point<double, TDimension> pt;
370  for ( unsigned int d = 0; d < dim; d++ )
371  {
372  pt[d] = v(d);
373  }
374  m_ControlPoints.push_back(pt);
375  }
376  return;
377 }
378 
379 template< unsigned int TDimension >
380 typename NonUniformBSpline<TDimension>::PointType
381 NonUniformBSpline< TDimension >
382 ::EvaluateSpline(const itk::Array<double> & p) const
383 {
384  double t = p[0];
385 
386  return EvaluateSpline(t);
387 }
388 
389 template< unsigned int TDimension >
390 typename NonUniformBSpline<TDimension>::PointType
391 NonUniformBSpline< TDimension >
392 ::EvaluateSpline(double t) const
393 {
394  int i = 0;
395 
396  vnl_vector<double> result(TDimension);
397  result.fill(0);
398 
399  for (typename ControlPointListType::const_iterator cpiter = m_ControlPoints.begin();
400  cpiter != m_ControlPoints.end();
401  cpiter++)
402  {
403  ControlPointType pt = *cpiter;
404  vnl_vector<double> v = pt.GetVnlVector();
405 
406  const double N = this->NonUniformBSplineFunctionRecursive(m_SplineOrder, i, t);
407 
408  for( unsigned j = 0; j < TDimension; j++ )
409  {
410  result[j] += N * v[j];
411  }
412 
413  i++;
414  }
415 
416  double array[TDimension];
417  for ( unsigned int d = 0; d < TDimension; d++ )
418  {
419  array[d] = result[d];
420  }
421 
422  ControlPointType sum(array);
423 #ifdef DEBUG_SPLINE
424  std::cout << "Result : " << result << std::endl;
425  std::cout << "Sum : " << sum << std::endl;
426 #endif
427 
428  return sum;
429 }
430 
431 } // end namespace itk
432 
433 #endif