Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
itkStreamlineTrackingFilter.cpp
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 #ifndef __itkStreamlineTrackingFilter_txx
18 #define __itkStreamlineTrackingFilter_txx
19 
20 #include <time.h>
21 #include <stdio.h>
22 #include <stdlib.h>
23 
24 #include "itkStreamlineTrackingFilter.h"
25 #include <itkImageRegionConstIterator.h>
26 #include <itkImageRegionConstIteratorWithIndex.h>
27 #include <itkImageRegionIterator.h>
28 
29 #define _USE_MATH_DEFINES
30 #include <math.h>
31 
32 namespace itk {
33 
34 template< class TTensorPixelType, class TPDPixelType>
35 StreamlineTrackingFilter< TTensorPixelType,
36 TPDPixelType>
37 ::StreamlineTrackingFilter()
38  : m_FiberPolyData(NULL)
39  , m_Points(NULL)
40  , m_Cells(NULL)
41  , m_FaImage(NULL)
42  , m_NumberOfInputs(1)
43  , m_FaThreshold(0.2)
44  , m_MinCurvatureRadius(0)
45  , m_StepSize(0)
46  , m_MaxLength(10000)
47  , m_MinTractLength(0.0)
48  , m_SeedsPerVoxel(1)
49  , m_F(1.0)
50  , m_G(0.0)
51  , m_Interpolate(true)
52  , m_PointPistance(0.0)
53  , m_ResampleFibers(false)
54  , m_SeedImage(NULL)
55  , m_MaskImage(NULL)
56 {
57  // At least 1 inputs is necessary for a vector image.
58  // For images added one at a time we need at least six
59  this->SetNumberOfRequiredInputs( 1 );
60  this->SetNumberOfIndexedInputs(3);
61 }
62 
63 template< class TTensorPixelType,
64  class TPDPixelType>
65 double StreamlineTrackingFilter< TTensorPixelType,
66 TPDPixelType>
67 ::RoundToNearest(double num) {
68  return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5);
69 }
70 
71 template< class TTensorPixelType,
72  class TPDPixelType>
73 void StreamlineTrackingFilter< TTensorPixelType,
74 TPDPixelType>
75 ::BeforeThreadedGenerateData()
76 {
77  m_FiberPolyData = FiberPolyDataType::New();
78  m_Points = vtkSmartPointer< vtkPoints >::New();
79  m_Cells = vtkSmartPointer< vtkCellArray >::New();
80 
81  InputImageType* inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
82  m_ImageSize.resize(3);
83  m_ImageSize[0] = inputImage->GetLargestPossibleRegion().GetSize()[0];
84  m_ImageSize[1] = inputImage->GetLargestPossibleRegion().GetSize()[1];
85  m_ImageSize[2] = inputImage->GetLargestPossibleRegion().GetSize()[2];
86 
87  if (m_ImageSize[0]<3 || m_ImageSize[1]<3 || m_ImageSize[2]<3)
88  m_Interpolate = false;
89 
90  m_ImageSpacing.resize(3);
91  m_ImageSpacing[0] = inputImage->GetSpacing()[0];
92  m_ImageSpacing[1] = inputImage->GetSpacing()[1];
93  m_ImageSpacing[2] = inputImage->GetSpacing()[2];
94 
95  double minSpacing;
96  if(m_ImageSpacing[0]<m_ImageSpacing[1] && m_ImageSpacing[0]<m_ImageSpacing[2])
97  minSpacing = m_ImageSpacing[0];
98  else if (m_ImageSpacing[1] < m_ImageSpacing[2])
99  minSpacing = m_ImageSpacing[1];
100  else
101  minSpacing = m_ImageSpacing[2];
102  if (m_StepSize<0.1*minSpacing)
103  m_StepSize = 0.1*minSpacing;
104 
105  if (m_ResampleFibers)
106  m_PointPistance = 0.5*minSpacing;
107 
108  m_PolyDataContainer = itk::VectorContainer< int, FiberPolyDataType >::New();
109  for (unsigned int i=0; i<this->GetNumberOfThreads(); i++)
110  {
111  FiberPolyDataType poly = FiberPolyDataType::New();
112  m_PolyDataContainer->InsertElement(i, poly);
113  }
114 
115  if (m_SeedImage.IsNull())
116  {
117  // initialize mask image
118  m_SeedImage = ItkUcharImgType::New();
119  m_SeedImage->SetSpacing( inputImage->GetSpacing() );
120  m_SeedImage->SetOrigin( inputImage->GetOrigin() );
121  m_SeedImage->SetDirection( inputImage->GetDirection() );
122  m_SeedImage->SetRegions( inputImage->GetLargestPossibleRegion() );
123  m_SeedImage->Allocate();
124  m_SeedImage->FillBuffer(1);
125  }
126 
127  if (m_MaskImage.IsNull())
128  {
129  // initialize mask image
130  m_MaskImage = ItkUcharImgType::New();
131  m_MaskImage->SetSpacing( inputImage->GetSpacing() );
132  m_MaskImage->SetOrigin( inputImage->GetOrigin() );
133  m_MaskImage->SetDirection( inputImage->GetDirection() );
134  m_MaskImage->SetRegions( inputImage->GetLargestPossibleRegion() );
135  m_MaskImage->Allocate();
136  m_MaskImage->FillBuffer(1);
137  }
138 
139  bool useUserFaImage = true;
140  if (m_FaImage.IsNull())
141  {
142  m_FaImage = ItkFloatImgType::New();
143  m_FaImage->SetSpacing( inputImage->GetSpacing() );
144  m_FaImage->SetOrigin( inputImage->GetOrigin() );
145  m_FaImage->SetDirection( inputImage->GetDirection() );
146  m_FaImage->SetRegions( inputImage->GetLargestPossibleRegion() );
147  m_FaImage->Allocate();
148  m_FaImage->FillBuffer(0.0);
149  useUserFaImage = false;
150  }
151 
152  m_NumberOfInputs = 0;
153  for (unsigned int i=0; i<this->GetNumberOfIndexedInputs(); i++)
154  {
155  if (this->ProcessObject::GetInput(i)==NULL)
156  break;
157 
158  ItkPDImgType::Pointer pdImage = ItkPDImgType::New();
159  pdImage->SetSpacing( inputImage->GetSpacing() );
160  pdImage->SetOrigin( inputImage->GetOrigin() );
161  pdImage->SetDirection( inputImage->GetDirection() );
162  pdImage->SetRegions( inputImage->GetLargestPossibleRegion() );
163  pdImage->Allocate();
164  m_PdImage.push_back(pdImage);
165 
166  ItkDoubleImgType::Pointer emaxImage = ItkDoubleImgType::New();
167  emaxImage->SetSpacing( inputImage->GetSpacing() );
168  emaxImage->SetOrigin( inputImage->GetOrigin() );
169  emaxImage->SetDirection( inputImage->GetDirection() );
170  emaxImage->SetRegions( inputImage->GetLargestPossibleRegion() );
171  emaxImage->Allocate();
172  emaxImage->FillBuffer(1.0);
173  m_EmaxImage.push_back(emaxImage);
174 
175  typename InputImageType::Pointer inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(i) );
176  m_InputImage.push_back(inputImage);
177 
178  m_NumberOfInputs++;
179  }
180  MITK_INFO << "Processing " << m_NumberOfInputs << " tensor files";
181 
182  typedef itk::DiffusionTensor3D<TTensorPixelType> TensorType;
183  typename TensorType::EigenValuesArrayType eigenvalues;
184  typename TensorType::EigenVectorsMatrixType eigenvectors;
185 
186 
187  for (int x=0; x<m_ImageSize[0]; x++)
188  for (int y=0; y<m_ImageSize[1]; y++)
189  for (int z=0; z<m_ImageSize[2]; z++)
190  {
191  typename InputImageType::IndexType index;
192  index[0] = x; index[1] = y; index[2] = z;
193  for (int i=0; i<m_NumberOfInputs; i++)
194  {
195  typename InputImageType::PixelType tensor = m_InputImage.at(i)->GetPixel(index);
196  vnl_vector_fixed<double,3> dir;
197  tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors);
198  dir[0] = eigenvectors(2, 0);
199  dir[1] = eigenvectors(2, 1);
200  dir[2] = eigenvectors(2, 2);
201  dir.normalize();
202  m_PdImage.at(i)->SetPixel(index, dir);
203  if (!useUserFaImage)
204  m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)+tensor.GetFractionalAnisotropy());
205  m_EmaxImage.at(i)->SetPixel(index, 2/eigenvalues[2]);
206  }
207  if (!useUserFaImage)
208  m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)/m_NumberOfInputs);
209  }
210 
211  if (m_Interpolate)
212  std::cout << "StreamlineTrackingFilter: using trilinear interpolation" << std::endl;
213  else
214  {
215  if (m_MinCurvatureRadius<0.0)
216  m_MinCurvatureRadius = 0.1*minSpacing;
217 
218  std::cout << "StreamlineTrackingFilter: using nearest neighbor interpolation" << std::endl;
219  }
220 
221  if (m_MinCurvatureRadius<0.0)
222  m_MinCurvatureRadius = 0.5*minSpacing;
223 
224  std::cout << "StreamlineTrackingFilter: Min. curvature radius: " << m_MinCurvatureRadius << std::endl;
225  std::cout << "StreamlineTrackingFilter: FA threshold: " << m_FaThreshold << std::endl;
226  std::cout << "StreamlineTrackingFilter: stepsize: " << m_StepSize << " mm" << std::endl;
227  std::cout << "StreamlineTrackingFilter: f: " << m_F << std::endl;
228  std::cout << "StreamlineTrackingFilter: g: " << m_G << std::endl;
229  std::cout << "StreamlineTrackingFilter: starting streamline tracking using " << this->GetNumberOfThreads() << " threads." << std::endl;
230 }
231 
232 template< class TTensorPixelType, class TPDPixelType>
233 void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType>
234 ::CalculateNewPosition(itk::ContinuousIndex<double, 3>& pos, vnl_vector_fixed<double,3>& dir, typename InputImageType::IndexType& index)
235 {
236  vnl_matrix_fixed< double, 3, 3 > rot = m_InputImage.at(0)->GetDirection().GetTranspose();
237  dir = rot*dir;
238  if (true)
239  {
240  dir *= m_StepSize;
241  pos[0] += dir[0]/m_ImageSpacing[0];
242  pos[1] += dir[1]/m_ImageSpacing[1];
243  pos[2] += dir[2]/m_ImageSpacing[2];
244  index[0] = RoundToNearest(pos[0]);
245  index[1] = RoundToNearest(pos[1]);
246  index[2] = RoundToNearest(pos[2]);
247  }
248  else
249  {
250  dir[0] /= m_ImageSpacing[0];
251  dir[1] /= m_ImageSpacing[1];
252  dir[2] /= m_ImageSpacing[2];
253 
254  int smallest = 0;
255  double x = 100000;
256  if (dir[0]>0)
257  {
258  if (fabs(fabs(RoundToNearest(pos[0])-pos[0])-0.5)>mitk::eps)
259  x = fabs(pos[0]-RoundToNearest(pos[0])-0.5)/dir[0];
260  else
261  x = fabs(pos[0]-std::ceil(pos[0])-0.5)/dir[0];
262  }
263  else if (dir[0]<0)
264  {
265  if (fabs(fabs(RoundToNearest(pos[0])-pos[0])-0.5)>mitk::eps)
266  x = -fabs(pos[0]-RoundToNearest(pos[0])+0.5)/dir[0];
267  else
268  x = -fabs(pos[0]-std::floor(pos[0])+0.5)/dir[0];
269  }
270  double s = x;
271 
272  double y = 100000;
273  if (dir[1]>0)
274  {
275  if (fabs(fabs(RoundToNearest(pos[1])-pos[1])-0.5)>mitk::eps)
276  y = fabs(pos[1]-RoundToNearest(pos[1])-0.5)/dir[1];
277  else
278  y = fabs(pos[1]-std::ceil(pos[1])-0.5)/dir[1];
279  }
280  else if (dir[1]<0)
281  {
282  if (fabs(fabs(RoundToNearest(pos[1])-pos[1])-0.5)>mitk::eps)
283  y = -fabs(pos[1]-RoundToNearest(pos[1])+0.5)/dir[1];
284  else
285  y = -fabs(pos[1]-std::floor(pos[1])+0.5)/dir[1];
286  }
287  if (s>y)
288  {
289  s=y;
290  smallest = 1;
291  }
292 
293  double z = 100000;
294  if (dir[2]>0)
295  {
296  if (fabs(fabs(RoundToNearest(pos[2])-pos[2])-0.5)>mitk::eps)
297  z = fabs(pos[2]-RoundToNearest(pos[2])-0.5)/dir[2];
298  else
299  z = fabs(pos[2]-std::ceil(pos[2])-0.5)/dir[2];
300  }
301  else if (dir[2]<0)
302  {
303  if (fabs(fabs(RoundToNearest(pos[2])-pos[2])-0.5)>mitk::eps)
304  z = -fabs(pos[2]-RoundToNearest(pos[2])+0.5)/dir[2];
305  else
306  z = -fabs(pos[2]-std::floor(pos[2])+0.5)/dir[2];
307  }
308  if (s>z)
309  {
310  s=z;
311  smallest = 2;
312  }
313 
314  // MITK_INFO << "---------------------------------------------";
315  // MITK_INFO << "s: " << s;
316  // MITK_INFO << "dir: " << dir;
317  // MITK_INFO << "old: " << pos[0] << ", " << pos[1] << ", " << pos[2];
318 
319  pos[0] += dir[0]*s;
320  pos[1] += dir[1]*s;
321  pos[2] += dir[2]*s;
322 
323  switch (smallest)
324  {
325  case 0:
326  if (dir[0]<0)
327  index[0] = std::floor(pos[0]);
328  else
329  index[0] = std::ceil(pos[0]);
330  index[1] = RoundToNearest(pos[1]);
331  index[2] = RoundToNearest(pos[2]);
332  break;
333 
334  case 1:
335  if (dir[1]<0)
336  index[1] = std::floor(pos[1]);
337  else
338  index[1] = std::ceil(pos[1]);
339  index[0] = RoundToNearest(pos[0]);
340  index[2] = RoundToNearest(pos[2]);
341  break;
342 
343  case 2:
344  if (dir[2]<0)
345  index[2] = std::floor(pos[2]);
346  else
347  index[2] = std::ceil(pos[2]);
348  index[1] = RoundToNearest(pos[1]);
349  index[0] = RoundToNearest(pos[0]);
350  }
351 
352  // double x = 100000;
353  // if (dir[0]>0)
354  // x = fabs(pos[0]-RoundToNearest(pos[0])-0.5)/dir[0];
355  // else if (dir[0]<0)
356  // x = -fabs(pos[0]-RoundToNearest(pos[0])+0.5)/dir[0];
357  // double s = x;
358 
359  // double y = 100000;
360  // if (dir[1]>0)
361  // y = fabs(pos[1]-RoundToNearest(pos[1])-0.5)/dir[1];
362  // else if (dir[1]<0)
363  // y = -fabs(pos[1]-RoundToNearest(pos[1])+0.5)/dir[1];
364  // if (s>y)
365  // s=y;
366 
367  // double z = 100000;
368  // if (dir[2]>0)
369  // z = fabs(pos[2]-RoundToNearest(pos[2])-0.5)/dir[2];
370  // else if (dir[2]<0)
371  // z = -fabs(pos[2]-RoundToNearest(pos[2])+0.5)/dir[2];
372 
373  // if (s>z)
374  // s=z;
375  // s *= 1.001;
376 
377  // pos[0] += dir[0]*s;
378  // pos[1] += dir[1]*s;
379  // pos[2] += dir[2]*s;
380 
381  // index[0] = RoundToNearest(pos[0]);
382  // index[1] = RoundToNearest(pos[1]);
383  // index[2] = RoundToNearest(pos[2]);
384 
385  // MITK_INFO << "new: " << pos[0] << ", " << pos[1] << ", " << pos[2];
386  }
387 }
388 
389 template< class TTensorPixelType, class TPDPixelType>
390 bool StreamlineTrackingFilter< TTensorPixelType, TPDPixelType>
391 ::IsValidPosition(itk::ContinuousIndex<double, 3>& pos, typename InputImageType::IndexType &index, vnl_vector_fixed< double, 8 >& interpWeights, int imageIdx)
392 {
393  if (!m_InputImage.at(imageIdx)->GetLargestPossibleRegion().IsInside(index) || m_MaskImage->GetPixel(index)==0)
394  return false;
395 
396  if (m_Interpolate)
397  {
398  double frac_x = pos[0] - index[0];
399  double frac_y = pos[1] - index[1];
400  double frac_z = pos[2] - index[2];
401 
402  if (frac_x<0)
403  {
404  index[0] -= 1;
405  frac_x += 1;
406  }
407  if (frac_y<0)
408  {
409  index[1] -= 1;
410  frac_y += 1;
411  }
412  if (frac_z<0)
413  {
414  index[2] -= 1;
415  frac_z += 1;
416  }
417 
418  frac_x = 1-frac_x;
419  frac_y = 1-frac_y;
420  frac_z = 1-frac_z;
421 
422  // int coordinates inside image?
423  if (index[0] < 0 || index[0] >= m_ImageSize[0]-1)
424  return false;
425  if (index[1] < 0 || index[1] >= m_ImageSize[1]-1)
426  return false;
427  if (index[2] < 0 || index[2] >= m_ImageSize[2]-1)
428  return false;
429 
430  interpWeights[0] = ( frac_x)*( frac_y)*( frac_z);
431  interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z);
432  interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z);
433  interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z);
434  interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z);
435  interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z);
436  interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z);
437  interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z);
438 
439  typename InputImageType::IndexType tmpIdx;
440  double FA = m_FaImage->GetPixel(index) * interpWeights[0];
441  tmpIdx = index; tmpIdx[0]++;
442  FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[1];
443  tmpIdx = index; tmpIdx[1]++;
444  FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[2];
445  tmpIdx = index; tmpIdx[2]++;
446  FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[3];
447  tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++;
448  FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[4];
449  tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++;
450  FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[5];
451  tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++;
452  FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[6];
453  tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++;
454  FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[7];
455 
456  if (FA<m_FaThreshold)
457  return false;
458  }
459  else if (m_FaImage->GetPixel(index)<m_FaThreshold)
460  return false;
461 
462  return true;
463 }
464 
465 template< class TTensorPixelType, class TPDPixelType>
466 double StreamlineTrackingFilter< TTensorPixelType, TPDPixelType>
467 ::FollowStreamline(itk::ContinuousIndex<double, 3> pos, int dirSign, vtkPoints* points, std::vector< vtkIdType >& ids, int imageIdx)
468 {
469  double tractLength = 0;
470  typedef itk::DiffusionTensor3D<TTensorPixelType> TensorType;
471  typename TensorType::EigenValuesArrayType eigenvalues;
472  typename TensorType::EigenVectorsMatrixType eigenvectors;
473  vnl_vector_fixed< double, 8 > interpWeights;
474 
475  typename InputImageType::IndexType index, indexOld;
476  indexOld[0] = -1; indexOld[1] = -1; indexOld[2] = -1;
477  itk::Point<double> worldPos;
478  double distance = 0;
479  double distanceInVoxel = 0;
480 
481  // starting index and direction
482  index[0] = RoundToNearest(pos[0]);
483  index[1] = RoundToNearest(pos[1]);
484  index[2] = RoundToNearest(pos[2]);
485 
486  vnl_vector_fixed<double,3> dir = m_PdImage.at(imageIdx)->GetPixel(index);
487  dir *= dirSign; // reverse direction
488  vnl_vector_fixed<double,3> dirOld = dir;
489  if (dir.magnitude()<mitk::eps)
490  return tractLength;
491 
492  for (int step=0; step< m_MaxLength/2; step++)
493  {
494  // get new position
495  CalculateNewPosition(pos, dir, index);
496  distance += m_StepSize;
497 
498  // is new position valid (inside image, above FA threshold etc.)
499  if (!IsValidPosition(pos, index, interpWeights, imageIdx)) // if not end streamline
500  {
501 // m_SeedImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos );
502 // ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer()));
503  return tractLength;
504  }
505  else if (distance>=m_PointPistance)
506  {
507  tractLength += m_StepSize;
508  distanceInVoxel += m_StepSize;
509  m_SeedImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos );
510  ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer()));
511  distance = 0;
512  }
513 
514  if (!m_Interpolate) // use nearest neighbour interpolation
515  {
516  if (indexOld!=index) // did we enter a new voxel? if yes, calculate new direction
517  {
518  double minAngle = 0;
519  for (int img=0; img<m_NumberOfInputs; img++)
520  {
521  vnl_vector_fixed<double,3> newDir = m_PdImage.at(img)->GetPixel(index); // get principal direction
522  if (newDir.magnitude()<mitk::eps)
523  continue;
524 
525  typename InputImageType::PixelType tensor = m_InputImage.at(img)->GetPixel(index);
526  double scale = m_EmaxImage.at(img)->GetPixel(index);
527  newDir[0] = m_F*newDir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + scale*m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2]));
528  newDir[1] = m_F*newDir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + scale*m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2]));
529  newDir[2] = m_F*newDir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + scale*m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2]));
530  newDir.normalize();
531 
532  double angle = dot_product(dirOld, newDir);
533  if (angle<0)
534  {
535  newDir *= -1;
536  angle *= -1;
537  }
538 
539  if (angle>minAngle)
540  {
541  minAngle = angle;
542  dir = newDir;
543  }
544  }
545 
546  //double r = m_StepSize/(2*std::asin(std::acos(minAngle)/2));
547  vnl_vector_fixed<double,3> v3 = dir+dirOld; v3 *= m_StepSize;
548  double a = m_StepSize;
549  double b = m_StepSize;
550  double c = v3.magnitude();
551  double r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle)
552 
553  if (r<m_MinCurvatureRadius)
554  {
555  return tractLength;
556  }
557  dirOld = dir;
558  indexOld = index;
559  distanceInVoxel = 0;
560  }
561  else
562  dir = dirOld;
563  }
564  else // use trilinear interpolation (weights calculated in IsValidPosition())
565  {
566  typename InputImageType::PixelType tensor;
567 
568  typename InputImageType::IndexType tmpIdx = index;
569  typename InputImageType::PixelType tmpTensor;
570 
571  if (m_NumberOfInputs>1)
572  {
573  double minAngle = 0;
574  for (int img=0; img<m_NumberOfInputs; img++)
575  {
576  double angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
577  if (fabs(angle)>minAngle)
578  {
579  minAngle = angle;
580  tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
581  }
582  }
583  tensor = tmpTensor * interpWeights[0];
584 
585  minAngle = 0;
586  tmpIdx = index; tmpIdx[0]++;
587  for (int img=0; img<m_NumberOfInputs; img++)
588  {
589  double angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
590  if (fabs(angle)>minAngle)
591  {
592  minAngle = angle;
593  tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
594  }
595  }
596  tensor += tmpTensor * interpWeights[1];
597 
598  minAngle = 0;
599  tmpIdx = index; tmpIdx[1]++;
600  for (int img=0; img<m_NumberOfInputs; img++)
601  {
602  double angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
603  if (fabs(angle)>minAngle)
604  {
605  minAngle = angle;
606  tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
607  }
608  }
609  tensor += tmpTensor * interpWeights[2];
610 
611  minAngle = 0;
612  tmpIdx = index; tmpIdx[2]++;
613  for (int img=0; img<m_NumberOfInputs; img++)
614  {
615  double angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
616  if (fabs(angle)>minAngle)
617  {
618  minAngle = angle;
619  tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
620  }
621  }
622  tensor += tmpTensor * interpWeights[3];
623 
624  minAngle = 0;
625  tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++;
626  for (int img=0; img<m_NumberOfInputs; img++)
627  {
628  double angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
629  if (fabs(angle)>minAngle)
630  {
631  minAngle = angle;
632  tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
633  }
634  }
635  tensor += tmpTensor * interpWeights[4];
636 
637  minAngle = 0;
638  tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++;
639  for (int img=0; img<m_NumberOfInputs; img++)
640  {
641  double angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
642  if (fabs(angle)>minAngle)
643  {
644  minAngle = angle;
645  tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
646  }
647  }
648  tensor += tmpTensor * interpWeights[5];
649 
650  minAngle = 0;
651  tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++;
652  for (int img=0; img<m_NumberOfInputs; img++)
653  {
654  double angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
655  if (fabs(angle)>minAngle)
656  {
657  minAngle = angle;
658  tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
659  }
660  }
661  tensor += tmpTensor * interpWeights[6];
662 
663  minAngle = 0;
664  tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++;
665  for (int img=0; img<m_NumberOfInputs; img++)
666  {
667  double angle = dot_product(dirOld, m_PdImage.at(img)->GetPixel(tmpIdx));
668  if (fabs(angle)>minAngle)
669  {
670  minAngle = angle;
671  tmpTensor = m_InputImage.at(img)->GetPixel(tmpIdx);
672  }
673  }
674  tensor += tmpTensor * interpWeights[7];
675  }
676  else
677  {
678  tensor = m_InputImage.at(0)->GetPixel(index) * interpWeights[0];
679  typename InputImageType::IndexType tmpIdx = index; tmpIdx[0]++;
680  tensor += m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[1];
681  tmpIdx = index; tmpIdx[1]++;
682  tensor += m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[2];
683  tmpIdx = index; tmpIdx[2]++;
684  tensor += m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[3];
685  tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++;
686  tensor += m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[4];
687  tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++;
688  tensor += m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[5];
689  tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++;
690  tensor += m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[6];
691  tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++;
692  tensor += m_InputImage.at(0)->GetPixel(tmpIdx) * interpWeights[7];
693  }
694 
695  tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors);
696  dir[0] = eigenvectors(2, 0);
697  dir[1] = eigenvectors(2, 1);
698  dir[2] = eigenvectors(2, 2);
699  if (dir.magnitude()<mitk::eps)
700  continue;
701  dir.normalize();
702 
703  double scale = 2/eigenvalues[2];
704  dir[0] = m_F*dir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + scale*m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2]));
705  dir[1] = m_F*dir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + scale*m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2]));
706  dir[2] = m_F*dir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + scale*m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2]));
707  dir.normalize();
708 
709  double angle = dot_product(dirOld, dir);
710  if (angle<0)
711  {
712  dir *= -1;
713  angle *= -1;
714  }
715 
716  vnl_vector_fixed<double,3> v3 = dir+dirOld; v3 *= m_StepSize;
717  double a = m_StepSize;
718  double b = m_StepSize;
719  double c = v3.magnitude();
720  double r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle)
721 
722  if (r<m_MinCurvatureRadius)
723  return tractLength;
724 
725  dirOld = dir;
726  indexOld = index;
727  }
728  }
729  return tractLength;
730 }
731 
732 template< class TTensorPixelType,
733  class TPDPixelType>
734 void StreamlineTrackingFilter< TTensorPixelType,
735 TPDPixelType>
736 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
737  ThreadIdType threadId)
738 {
739  FiberPolyDataType poly = m_PolyDataContainer->GetElement(threadId);
740  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
741  vtkSmartPointer<vtkCellArray> Cells = vtkSmartPointer<vtkCellArray>::New();
742 
743  typedef itk::DiffusionTensor3D<TTensorPixelType> TensorType;
744  typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
745  typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType;
746  typedef ImageRegionConstIterator< ItkFloatImgType > doubleIteratorType;
747  typedef typename InputImageType::PixelType InputTensorType;
748 
749  MaskIteratorType sit(m_SeedImage, outputRegionForThread );
750  doubleIteratorType fit(m_FaImage, outputRegionForThread );
751  MaskIteratorType mit(m_MaskImage, outputRegionForThread );
752 
753  for (int img=0; img<m_NumberOfInputs; img++)
754  {
755  sit.GoToBegin();
756  mit.GoToBegin();
757  fit.GoToBegin();
758  itk::Point<double> worldPos;
759  while( !sit.IsAtEnd() )
760  {
761  if (sit.Value()==0 || fit.Value()<m_FaThreshold || mit.Value()==0)
762  {
763  ++sit;
764  ++mit;
765  ++fit;
766  continue;
767  }
768 
769  for (int s=0; s<m_SeedsPerVoxel; s++)
770  {
771  vtkSmartPointer<vtkPolyLine> line = vtkSmartPointer<vtkPolyLine>::New();
772  std::vector< vtkIdType > pointIDs;
773  typename InputImageType::IndexType index = sit.GetIndex();
774  itk::ContinuousIndex<double, 3> start;
775  unsigned int counter = 0;
776 
777  if (m_SeedsPerVoxel>1)
778  {
779  start[0] = index[0]+(double)(rand()%99-49)/100;
780  start[1] = index[1]+(double)(rand()%99-49)/100;
781  start[2] = index[2]+(double)(rand()%99-49)/100;
782  }
783  else
784  {
785  start[0] = index[0];
786  start[1] = index[1];
787  start[2] = index[2];
788  }
789 
790  // forward tracking
791  double tractLength = FollowStreamline(start, 1, points, pointIDs, img);
792 
793  // add ids to line
794  counter += pointIDs.size();
795  while (!pointIDs.empty())
796  {
797  line->GetPointIds()->InsertNextId(pointIDs.back());
798  pointIDs.pop_back();
799  }
800 
801  // insert start point
802  m_SeedImage->TransformContinuousIndexToPhysicalPoint( start, worldPos );
803  line->GetPointIds()->InsertNextId(points->InsertNextPoint(worldPos.GetDataPointer()));
804 
805  // backward tracking
806  tractLength += FollowStreamline(start, -1, points, pointIDs, img);
807 
808  counter += pointIDs.size();
809 
810  //MITK_INFO << "Tract length " << tractLength;
811 
812  if (tractLength<m_MinTractLength || counter<2)
813  continue;
814 
815  // add ids to line
816  for (unsigned int i=0; i<pointIDs.size(); i++)
817  line->GetPointIds()->InsertNextId(pointIDs.at(i));
818 
819  Cells->InsertNextCell(line);
820  }
821  ++sit;
822  ++mit;
823  ++fit;
824  }
825  }
826  poly->SetPoints(points);
827  poly->SetLines(Cells);
828 
829  std::cout << "Thread " << threadId << " finished tracking" << std::endl;
830 }
831 
832 template< class TTensorPixelType,
833  class TPDPixelType>
834 vtkSmartPointer< vtkPolyData > StreamlineTrackingFilter< TTensorPixelType,
835 TPDPixelType>
836 ::AddPolyData(FiberPolyDataType poly1, FiberPolyDataType poly2)
837 {
838  vtkSmartPointer<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
839  vtkSmartPointer<vtkCellArray> vNewLines = poly1->GetLines();
840  vtkSmartPointer<vtkPoints> vNewPoints = poly1->GetPoints();
841 
842  for( int i=0; i<poly2->GetNumberOfLines(); i++ )
843  {
844  vtkCell* cell = poly2->GetCell(i);
845  int numPoints = cell->GetNumberOfPoints();
846  vtkPoints* points = cell->GetPoints();
847 
848  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
849  for (int j=0; j<numPoints; j++)
850  {
851  double p[3];
852  points->GetPoint(j, p);
853 
854  vtkIdType id = vNewPoints->InsertNextPoint(p);
855  container->GetPointIds()->InsertNextId(id);
856  }
857  vNewLines->InsertNextCell(container);
858  }
859 
860  // initialize polydata
861  vNewPolyData->SetPoints(vNewPoints);
862  vNewPolyData->SetLines(vNewLines);
863 
864  return vNewPolyData;
865 }
866 template< class TTensorPixelType,
867  class TPDPixelType>
868 void StreamlineTrackingFilter< TTensorPixelType,
869 TPDPixelType>
870 ::AfterThreadedGenerateData()
871 {
872  MITK_INFO << "Generating polydata ";
873  m_FiberPolyData = m_PolyDataContainer->GetElement(0);
874  for (unsigned int i=1; i<this->GetNumberOfThreads(); i++)
875  {
876  m_FiberPolyData = AddPolyData(m_FiberPolyData, m_PolyDataContainer->GetElement(i));
877  }
878  MITK_INFO << "done";
879 }
880 
881 template< class TTensorPixelType,
882  class TPDPixelType>
883 void StreamlineTrackingFilter< TTensorPixelType,
884 TPDPixelType>
885 ::PrintSelf(std::ostream&, Indent) const
886 {
887 }
888 
889 }
890 #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx
line
static char * line
Definition: svm.cpp:2884
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
itk::StreamlineTrackingFilter
Performes deterministic streamline tracking on the input tensor image.
Definition: itkStreamlineTrackingFilter.h:42
itk::StreamlineTrackingFilter::CalculateNewPosition
void CalculateNewPosition(itk::ContinuousIndex< double, 3 > &pos, vnl_vector_fixed< double, 3 > &dir, typename InputImageType::IndexType &index)
Calculate next integration step.
Definition: itkStreamlineTrackingFilter.cpp:234
MITK_INFO
#define MITK_INFO
Definition: mitkLogMacros.h:22
itk::StreamlineTrackingFilter::FollowStreamline
double FollowStreamline(itk::ContinuousIndex< double, 3 > pos, int dirSign, vtkPoints *points, std::vector< vtkIdType > &ids, int imageIdx)
Start streamline in one direction.
Definition: itkStreamlineTrackingFilter.cpp:467
itk::StreamlineTrackingFilter::InputImageType
Superclass::InputImageType InputImageType
Definition: itkStreamlineTrackingFilter.h:63
TTensorPixelType
double TTensorPixelType
Definition: TensorDerivedMapsExtraction.cpp:40
itk::StreamlineTrackingFilter::FiberPolyDataType
vtkSmartPointer< vtkPolyData > FiberPolyDataType
Definition: itkStreamlineTrackingFilter.h:70
itk::StreamlineTrackingFilter::IsValidPosition
bool IsValidPosition(itk::ContinuousIndex< double, 3 > &pos, typename InputImageType::IndexType &index, vnl_vector_fixed< double, 8 > &interpWeights, int imageIdx)
Are we outside of the mask image? Is the FA too low?
Definition: itkStreamlineTrackingFilter.cpp:391
mitk::eps
MITKCORE_EXPORT const ScalarType eps
Definition: mitkNumericConstants.cpp:21
PixelType
unsigned short PixelType
Definition: TumorProgressionMapping.cpp:55
itk::StreamlineTrackingFilter::OutputImageRegionType
Superclass::OutputImageRegionType OutputImageRegionType
Definition: itkStreamlineTrackingFilter.h:65
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.