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