Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
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.