1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
17 /*===================================================================
19 This file is based heavily on a corresponding ITK filter.
21 ===================================================================*/
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 *********************************/
27 #ifndef __itkNonUniformBSpline_txx
28 #define __itkNonUniformBSpline_txx
31 #pragma warning ( disable : 4786 )
34 #include "mitkItkNonUniformBSpline.h"
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"
42 // #define DEBUG_SPLINE
48 template< unsigned int TDimension >
49 NonUniformBSpline< TDimension >
52 // Cubic bspline => 4th order
54 m_SpatialDimension = TDimension;
58 template< unsigned int TDimension >
59 NonUniformBSpline< TDimension >
60 ::~NonUniformBSpline()
64 /** Print the object */
65 template< unsigned int TDimension >
67 NonUniformBSpline< TDimension >
68 ::PrintSelf( std::ostream& os, Indent indent ) const
70 Superclass::PrintSelf( os, indent );
71 os << indent << "NonUniformBSpline(" << this << ")" << std::endl;
73 os << indent << "Chord lengths : " << std::endl;
74 for (ChordLengthListType::const_iterator iter = m_CumulativeChordLength.begin();
75 iter != m_CumulativeChordLength.end();
78 os << indent << indent << *iter << std::endl;
80 os << indent << "Knots : " << std::endl;
81 for (KnotListType::const_iterator kiter = m_Knots.begin();
82 kiter != m_Knots.end();
85 os << indent << indent << *kiter << std::endl;
87 os << indent << "Control Points : " << std::endl;
88 for (typename ControlPointListType::const_iterator cpiter = m_ControlPoints.begin();
89 cpiter != m_ControlPoints.end();
92 os << indent << indent << *cpiter << std::endl;
97 /** Set the list of points composing the tube */
98 template< unsigned int TDimension >
100 NonUniformBSpline< TDimension >
101 ::SetPoints( PointListType & points )
105 typename PointListType::iterator it,end;
110 m_Points.push_back(*it);
117 /** Set the list of points composing the tube */
118 template< unsigned int TDimension >
120 NonUniformBSpline< TDimension >
121 ::SetKnots( KnotListType & knots )
125 int len = knots.size();
126 double max_knot = knots[len - 1];
128 typename KnotListType::iterator it;
129 typename KnotListType::iterator end;
136 m_Knots.push_back(*it/max_knot);
143 template< unsigned int TDimension >
145 NonUniformBSpline< TDimension >
146 ::NonUniformBSplineFunctionRecursive(unsigned int order, unsigned int i, double t) const
150 if (m_Knots[i] <= t && t < m_Knots[i+1])
161 // Be careful, we must use the passed in parameter for the order since this
162 // function is recursive.
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)
169 else if (denom1 == 0)
170 std::cout << "Error : " << denom1 << ", " << numer1 << std::endl;
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)
177 else if (denom2 == 0)
178 std::cout << "Error : " << denom2 << ", " << numer2 << std::endl;
183 template< unsigned int TDimension >
185 NonUniformBSpline< TDimension >
186 ::ComputeChordLengths()
188 m_ChordLength.clear();
189 m_CumulativeChordLength.clear();
191 m_ChordLength.push_back(0);
192 m_CumulativeChordLength.push_back(0);
194 double total_chord_length = 0.0;
195 ChordLengthListType temp;
197 for (::size_t i = 0; i < m_Points.size()-1; i++)
199 PointType pt = m_Points[i];
200 PointType pt2 = m_Points[i+1];
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);
208 for (ChordLengthListType::iterator aiter = temp.begin();
212 m_CumulativeChordLength.push_back(*aiter/total_chord_length);
219 std::cout << "Total chord length : " << total_chord_length << std::endl;
221 std::cout << "Chord length : " << std::endl;
222 for (ChordLengthListType::iterator aiter2 = m_ChordLength.begin();
223 aiter2 != m_ChordLength.end();
226 std::cout << *aiter2 << std::endl;
229 std::cout << "Cumulative chord length : " << std::endl;
230 for (ChordLengthListType::iterator aiter3 = m_CumulativeChordLength.begin();
231 aiter3 != m_CumulativeChordLength.end();
234 std::cout << *aiter3 << std::endl;
236 std::cout << std::endl;
240 template< unsigned int TDimension >
242 NonUniformBSpline< TDimension >
243 ::SetControlPoints( ControlPointListType& ctrlpts )
245 m_ControlPoints.clear();
246 for (typename ControlPointListType::iterator iter = ctrlpts.begin();
247 iter != ctrlpts.end();
250 m_ControlPoints.push_back(*iter);
256 template< unsigned int TDimension >
258 NonUniformBSpline< TDimension >::ControlPointListType &
259 NonUniformBSpline< TDimension >::GetControlPoints() const
261 return this->m_ControlPoints;
265 template< unsigned int TDimension >
267 NonUniformBSpline< TDimension >::KnotListType &
268 NonUniformBSpline< TDimension >::GetKnots() const
270 return this->m_Knots;
274 template< unsigned int TDimension >
276 NonUniformBSpline< TDimension >::PointListType &
277 NonUniformBSpline< TDimension >::GetPoints() const
279 return this->m_Points;
283 template< unsigned int TDimension >
285 NonUniformBSpline< TDimension >::ComputeControlPoints()
287 unsigned int dim = m_Points[0].GetPointDimension();
290 std::cout << "Points have dimension : " << dim << std::endl;
294 // +1 in cols for radius
296 vnl_matrix<double> data_matrix(m_Points.size(), dim);
299 // Form data point matrix
302 for (typename PointListType::iterator iter = m_Points.begin();
303 iter != m_Points.end();
306 PointType pt = (*iter);
307 for (unsigned int i = 0; i < dim; i++)
309 data_matrix(rr, i) = pt.GetVnlVector()[i];
315 std::cout << std::endl << "Data matrix" << std::endl;
316 std::cout << data_matrix << std::endl;
320 // Form basis function matrix
322 //int num_basis_functions = 2 * m_SplineOrder - 1;
323 //int num_basis_functions = m_Points.size();
324 int num_rows = m_Points.size();
327 // Assumes multiplicity k (m_SplineOrder at the ends).
329 int num_cols = m_Knots.size() - m_SplineOrder;
331 vnl_matrix<double> N_matrix(num_rows, num_cols);
333 //N_matrix(0, 0) = 1.0;
335 for (int r = 0; r < num_rows; r++)
337 for (int c = 0; c < num_cols; c++)
339 double t = m_CumulativeChordLength[r];
340 N_matrix(r, c) = NonUniformBSplineFunctionRecursive(m_SplineOrder, c, t);
344 N_matrix(num_rows-1, num_cols-1) = 1.0;
347 std::cout << "Basis function matrix : " << std::endl;
348 std::cout << N_matrix << std::endl;
351 //FIXME: Use the LSQR linear solver here:
352 vnl_matrix<double> B;
354 // = vnl_matrix_inverse<double>(N_matrix.transpose() * N_matrix) * N_matrix.transpose() * data_matrix;
356 // vnl_linear_system ls( N_matrix.rows(), N_matrix.cols() );
358 // vnl_lsqr solver( ls );
361 //#ifdef DEBUG_SPLINE
362 std::cout << "Control point matrix : " << std::endl;
363 std::cout << B << std::endl;
366 m_ControlPoints.clear();
368 for ( unsigned int j = 0; j < B.rows(); j++ )
370 vnl_vector<double> v = B.get_row(j);
371 itk::Vector<double> iv;
373 itk::Point<double, TDimension> pt;
374 for ( unsigned int d = 0; d < dim; d++ )
378 m_ControlPoints.push_back(pt);
383 template< unsigned int TDimension >
384 typename NonUniformBSpline<TDimension>::PointType
385 NonUniformBSpline< TDimension >
386 ::EvaluateSpline(const itk::Array<double> & p) const
390 return EvaluateSpline(t);
393 template< unsigned int TDimension >
394 typename NonUniformBSpline<TDimension>::PointType
395 NonUniformBSpline< TDimension >
396 ::EvaluateSpline(double t) const
400 vnl_vector<double> result(TDimension);
403 for (typename ControlPointListType::const_iterator cpiter = m_ControlPoints.begin();
404 cpiter != m_ControlPoints.end();
407 ControlPointType pt = *cpiter;
408 vnl_vector<double> v = pt.GetVnlVector();
410 const double N = this->NonUniformBSplineFunctionRecursive(m_SplineOrder, i, t);
412 for( unsigned j = 0; j < TDimension; j++ )
414 result[j] += N * v[j];
420 double array[TDimension];
421 for ( unsigned int d = 0; d < TDimension; d++ )
423 array[d] = result[d];
426 ControlPointType sum(array);
428 std::cout << "Result : " << result << std::endl;
429 std::cout << "Sum : " << sum << std::endl;
435 } // end namespace itk