ANIMA  4.0
animaBaseTractographyImageFilter.cxx
Go to the documentation of this file.
2 
3 #include <itkImageRegionIteratorWithIndex.h>
4 #include <itkMultiThreaderBase.h>
5 #include <itkProgressReporter.h>
6 
8 
9 #include <vtkPointData.h>
10 #include <vtkCellData.h>
11 
12 namespace anima
13 {
14 
16 {
17  m_PointsToProcess.clear();
18 
19  m_NumberOfFibersPerPixel = 1;
20 
21  m_StepProgression = 1;
22  m_MaxFiberAngle = M_PI / 2.0;
23 
24  m_ComputeLocalColors = true;
25  m_HighestProcessedSeed = 0;
26  m_ProgressReport = ITK_NULLPTR;
27 }
28 
30 {
31  if (m_ProgressReport)
32  delete m_ProgressReport;
33 }
34 
36 {
37  this->PrepareTractography();
38  m_Output = vtkPolyData::New();
39 
40  if (m_ProgressReport)
41  delete m_ProgressReport;
42 
43  unsigned int stepData = std::min((int)m_PointsToProcess.size(),100);
44  if (stepData == 0)
45  stepData = 1;
46 
47  unsigned int numSteps = std::floor(m_PointsToProcess.size() / (double)stepData);
48  if (m_PointsToProcess.size() % stepData != 0)
49  numSteps++;
50 
51  m_ProgressReport = new itk::ProgressReporter(this,0,numSteps);
52 
53  std::vector < FiberType > resultFibers;
54 
55  trackerArguments tmpStr;
56  tmpStr.trackerPtr = this;
57  for (unsigned int i = 0;i < this->GetNumberOfWorkUnits();++i)
58  tmpStr.resultFibersFromThreads.push_back(resultFibers);
59 
60  this->GetMultiThreader()->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
61  this->GetMultiThreader()->SetSingleMethod(this->ThreadTracker,&tmpStr);
62  this->GetMultiThreader()->SingleMethodExecute();
63 
64  for (unsigned int j = 0;j < this->GetNumberOfWorkUnits();++j)
65  {
66  resultFibers.insert(resultFibers.end(),tmpStr.resultFibersFromThreads[j].begin(),
67  tmpStr.resultFibersFromThreads[j].end());
68  }
69 
70  std::cout << "\nTracked a total of " << resultFibers.size() << " fibers" << std::endl;
71  resultFibers = this->FilterOutputFibers(resultFibers);
72  std::cout << "Kept " << resultFibers.size() << " fibers after filtering" << std::endl;
73  this->createVTKOutput(resultFibers);
74 }
75 
77 {
78  typedef itk::ImageRegionIteratorWithIndex <MaskImageType> MaskImageIteratorType;
79 
80  MaskImageIteratorType seedItr(m_SeedingImage, m_InputImage->GetLargestPossibleRegion());
81  m_PointsToProcess.clear();
82 
83  IndexType tmpIndex;
84  PointType tmpPoint;
85  ContinuousIndexType realIndex;
86 
87  m_FilteringValues.clear();
88  double startN = -0.5 + 1.0 / (2.0 * m_NumberOfFibersPerPixel);
89  double stepN = 1.0 / m_NumberOfFibersPerPixel;
90 
91  bool is2d = (m_InputImage->GetLargestPossibleRegion().GetSize()[2] <= 1);
92  FiberType tmpFiber(1);
93 
94  while (!seedItr.IsAtEnd())
95  {
96  if (seedItr.Get() == 0)
97  {
98  ++seedItr;
99  continue;
100  }
101 
102  tmpIndex = seedItr.GetIndex();
103 
104  if (is2d)
105  {
106  realIndex[2] = tmpIndex[2];
107  for (unsigned int j = 0;j < m_NumberOfFibersPerPixel;++j)
108  {
109  realIndex[1] = tmpIndex[1] + startN + j * stepN;
110  for (unsigned int i = 0;i < m_NumberOfFibersPerPixel;++i)
111  {
112  realIndex[0] = tmpIndex[0] + startN + i * stepN;
113  m_SeedingImage->TransformContinuousIndexToPhysicalPoint(realIndex,tmpPoint);
114  tmpFiber[0] = tmpPoint;
115  m_PointsToProcess.push_back(std::pair<FiberProgressType,FiberType>(Both,tmpFiber));
116  }
117  }
118  }
119  else
120  {
121  for (unsigned int k = 0;k < m_NumberOfFibersPerPixel;++k)
122  {
123  realIndex[2] = tmpIndex[2] + startN + k * stepN;
124  for (unsigned int j = 0;j < m_NumberOfFibersPerPixel;++j)
125  {
126  realIndex[1] = tmpIndex[1] + startN + j * stepN;
127  for (unsigned int i = 0;i < m_NumberOfFibersPerPixel;++i)
128  {
129  realIndex[0] = tmpIndex[0] + startN + i * stepN;
130  m_SeedingImage->TransformContinuousIndexToPhysicalPoint(realIndex,tmpPoint);
131  tmpFiber[0] = tmpPoint;
132  m_PointsToProcess.push_back(std::pair<FiberProgressType,FiberType>(Both,tmpFiber));
133  }
134  }
135  }
136  }
137 
138  ++seedItr;
139  }
140 
141  if (m_FilteringImage.IsNotNull())
142  {
143  MaskImageIteratorType filterItr(m_FilteringImage, m_InputImage->GetLargestPossibleRegion());
144 
145  while (!filterItr.IsAtEnd())
146  {
147  if (filterItr.Get() == 0)
148  {
149  ++filterItr;
150  continue;
151  }
152 
153  bool isAlreadyIn = false;
154  for (unsigned int i = 0;i < m_FilteringValues.size();++i)
155  {
156  if (m_FilteringValues[i] == filterItr.Get())
157  {
158  isAlreadyIn = true;
159  break;
160  }
161  }
162 
163  if (!isAlreadyIn)
164  m_FilteringValues.push_back(filterItr.Get());
165 
166  ++filterItr;
167  }
168  }
169 
170  std::cout << "Generated " << m_PointsToProcess.size() << " seed points from ROI mask" << std::endl;
171 }
172 
173 ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION BaseTractographyImageFilter::ThreadTracker(void *arg)
174 {
175  itk::MultiThreaderBase::WorkUnitInfo *threadArgs = (itk::MultiThreaderBase::WorkUnitInfo *)arg;
176  unsigned int nbThread = threadArgs->WorkUnitID;
177 
178  trackerArguments *tmpArg = (trackerArguments *)threadArgs->UserData;
179  tmpArg->trackerPtr->ThreadTrack(nbThread,tmpArg->resultFibersFromThreads[nbThread]);
180 
181  return ITK_THREAD_RETURN_DEFAULT_VALUE;
182 }
183 
184 void BaseTractographyImageFilter::ThreadTrack(unsigned int numThread, std::vector <FiberType> &resultFibers)
185 {
186  bool continueLoop = true;
187  unsigned int highestToleratedSeedIndex = m_PointsToProcess.size();
188 
189  unsigned int stepData = std::min((int)m_PointsToProcess.size(),100);
190  if (stepData == 0)
191  stepData = 1;
192 
193  while (continueLoop)
194  {
195  m_LockHighestProcessedSeed.lock();
196 
197  if (m_HighestProcessedSeed >= highestToleratedSeedIndex)
198  {
199  m_LockHighestProcessedSeed.unlock();
200  continueLoop = false;
201  continue;
202  }
203 
204  unsigned int startPoint = m_HighestProcessedSeed;
205  unsigned int endPoint = m_HighestProcessedSeed + stepData;
206  if (endPoint > highestToleratedSeedIndex)
207  endPoint = highestToleratedSeedIndex;
208 
209  m_HighestProcessedSeed = endPoint;
210 
211  m_LockHighestProcessedSeed.unlock();
212 
213  this->ThreadedTrackComputer(numThread,resultFibers,startPoint,endPoint);
214 
215  m_LockHighestProcessedSeed.lock();
216  m_ProgressReport->CompletedPixel();
217  m_LockHighestProcessedSeed.unlock();
218  }
219 }
220 
221 void BaseTractographyImageFilter::ThreadedTrackComputer(unsigned int numThread, std::vector <FiberType> &resultFibers,
222  unsigned int startSeedIndex, unsigned int endSeedIndex)
223 {
224  FiberType tmpFiber;
225  for (unsigned int i = startSeedIndex;i < endSeedIndex;++i)
226  {
227  tmpFiber = m_PointsToProcess[i].second;
228  this->ComputeFiber(tmpFiber,m_PointsToProcess[i].first,numThread);
229 
230  if (tmpFiber.size() > m_MinLengthFiber / m_StepProgression)
231  resultFibers.push_back(tmpFiber);
232  }
233 }
234 
235 std::vector < std::vector <BaseTractographyImageFilter::PointType> >
236 BaseTractographyImageFilter::FilterOutputFibers(std::vector < FiberType > &fibers)
237 {
238  std::vector < FiberType > resVal;
239 
240  if (m_FilteringValues.size() > 1)
241  {
242  std::vector <bool> touchingLabels(m_FilteringValues.size());
243  IndexType tmpIndex;
244  PointType tmpPoint;
245  for (unsigned int i = 0;i < fibers.size();++i)
246  {
247  std::fill(touchingLabels.begin(),touchingLabels.end(),false);
248  for (unsigned int j = 0;j < fibers[i].size();++j)
249  {
250  tmpPoint = fibers[i][j];
251  m_FilteringImage->TransformPhysicalPointToIndex(tmpPoint,tmpIndex);
252 
253  unsigned int maskValue = m_FilteringImage->GetPixel(tmpIndex);
254  if (maskValue != 0)
255  {
256  for (unsigned int k = 0;k < m_FilteringValues.size();++k)
257  {
258  if (maskValue == m_FilteringValues[k])
259  {
260  touchingLabels[k] = true;
261  break;
262  }
263  }
264  }
265  }
266 
267  bool keepFiber = true;
268  for (unsigned int k = 0;k < touchingLabels.size();++k)
269  {
270  if (!touchingLabels[k])
271  {
272  keepFiber = false;
273  break;
274  }
275  }
276 
277  if (keepFiber)
278  resVal.push_back(fibers[i]);
279  }
280  }
281  else
282  resVal = fibers;
283 
284  return resVal;
285 }
286 
287 void BaseTractographyImageFilter::createVTKOutput(std::vector < FiberType > &filteredFibers)
288 {
289  m_Output = vtkPolyData::New();
290  m_Output->Initialize();
291  m_Output->Allocate();
292 
293  vtkSmartPointer <vtkPoints> myPoints = vtkPoints::New();
294  for (unsigned int i = 0;i < filteredFibers.size();++i)
295  {
296  unsigned int npts = filteredFibers[i].size();
297  vtkIdType* ids = new vtkIdType[npts];
298 
299  for (unsigned int j = 0;j < npts;++j)
300  ids[j] = myPoints->InsertNextPoint(filteredFibers[i][j][0],filteredFibers[i][j][1],filteredFibers[i][j][2]);
301 
302  m_Output->InsertNextCell (VTK_POLY_LINE, npts, ids);
303  delete[] ids;
304  }
305 
306  m_Output->SetPoints (myPoints);
307  if (m_ComputeLocalColors)
308  {
309  std::cout << "Computing local colors and microstructure maps" << std::endl;
311  }
312 }
313 
317  itk::ThreadIdType threadId)
318 {
319  FiberProcessVectorType resVal;
320  bool is2d = m_InputImage->GetLargestPossibleRegion().GetSize()[2] == 1;
321 
322  if ((ways == Forward)||(ways == Both))
323  {
324  PointType curPoint = fiber[fiber.size() - 1];
325  PointType oldDir, newDir;
326 
327  ContinuousIndexType curIndex;
328  IndexType curNearestIndex;
329  bool continueLoop = true;
330  m_SeedingImage->TransformPhysicalPointToContinuousIndex(curPoint,curIndex);
331  m_SeedingImage->TransformPhysicalPointToIndex(curPoint,curNearestIndex);
332 
333  VectorType modelValue;
334  modelValue.SetSize(1);
335 
336  if (!m_CutMaskImage.IsNull())
337  {
338  if (!this->CheckIndexInImageBounds(curNearestIndex,m_CutMaskImage))
339  continueLoop = false;
340  else if (m_CutMaskImage->GetPixel(curNearestIndex) != 0)
341  continueLoop = false;
342  }
343 
344  if (!m_ForbiddenMaskImage.IsNull())
345  {
346  if (!this->CheckIndexInImageBounds(curNearestIndex,m_ForbiddenMaskImage))
347  {
348  fiber.clear();
349  return resVal;
350  }
351  if (m_ForbiddenMaskImage->GetPixel(curNearestIndex) != 0)
352  {
353  fiber.clear();
354  return resVal;
355  }
356  }
357 
358  this->GetModelValue(curIndex,modelValue);
359  if (isZero(modelValue))
360  continueLoop = false;
361 
362  if (!this->CheckModelCompatibility(modelValue,threadId))
363  continueLoop = false;
364 
365  while (continueLoop)
366  {
367  if (fiber.size() <= 1)
368  oldDir = this->GetModelPrincipalDirection(modelValue,is2d,threadId);
369  else
370  {
371  for (unsigned int i = 0;i < 3;++i)
372  oldDir[i] = fiber[fiber.size() - 1][i] - fiber[fiber.size() - 2][i];
373  }
374 
375  // Do some stuff to progress one step
376  newDir = this->GetNextDirection(oldDir,modelValue,is2d,threadId);
377 
378  if (anima::ComputeOrientationAngle(oldDir, newDir) > m_MaxFiberAngle)
379  {
380  continueLoop = false;
381  break;
382  }
383 
384  if (anima::ComputeScalarProduct(oldDir, newDir) < 0)
385  anima::Revert(newDir,newDir);
386 
387  for (unsigned int i = 0;i < 3;++i)
388  curPoint[i] = fiber[fiber.size() - 1][i] + m_StepProgression * newDir[i];
389 
390  m_SeedingImage->TransformPhysicalPointToContinuousIndex(curPoint,curIndex);
391  m_SeedingImage->TransformPhysicalPointToIndex(curPoint,curNearestIndex);
392 
393  if (!m_CutMaskImage.IsNull())
394  {
395  if (!this->CheckIndexInImageBounds(curNearestIndex,m_CutMaskImage))
396  continueLoop = false;
397  else if (m_CutMaskImage->GetPixel(curNearestIndex) != 0)
398  continueLoop = false;
399  }
400 
401  if (!m_ForbiddenMaskImage.IsNull())
402  {
403  if (!this->CheckIndexInImageBounds(curNearestIndex,m_ForbiddenMaskImage))
404  {
405  fiber.clear();
406  return resVal;
407  }
408  if (m_ForbiddenMaskImage->GetPixel(curNearestIndex) != 0)
409  {
410  fiber.clear();
411  return resVal;
412  }
413  }
414 
415  if (!this->CheckIndexInImageBounds(curIndex))
416  continueLoop = false;
417  else
418  {
419  this->GetModelValue(curIndex,modelValue);
420 
421  if (isZero(modelValue))
422  continueLoop = false;
423  }
424 
425  if (!this->CheckModelCompatibility(modelValue,threadId))
426  continueLoop = false;
427 
428  // Add new point to fiber
429  if (continueLoop)
430  fiber.push_back(curPoint);
431 
432  if (fiber.size() > m_MaxLengthFiber / m_StepProgression)
433  {
434  fiber.clear();
435  return resVal;
436  }
437  }
438  }
439 
440  if ((ways == Backward)||(ways == Both))
441  {
442  PointType curPoint = fiber[0];
443  PointType oldDir, newDir;
444 
445  ContinuousIndexType curIndex;
446  IndexType curNearestIndex;
447  bool continueLoop = true;
448  m_SeedingImage->TransformPhysicalPointToContinuousIndex(curPoint,curIndex);
449  m_SeedingImage->TransformPhysicalPointToIndex(curPoint,curNearestIndex);
450 
451  VectorType modelValue;
452  modelValue.SetSize(1);
453 
454  if (!m_CutMaskImage.IsNull())
455  {
456  if (!this->CheckIndexInImageBounds(curNearestIndex,m_CutMaskImage))
457  continueLoop = false;
458  else if (m_CutMaskImage->GetPixel(curNearestIndex) != 0)
459  continueLoop = false;
460  }
461 
462  if (!m_ForbiddenMaskImage.IsNull())
463  {
464  if (!this->CheckIndexInImageBounds(curNearestIndex,m_ForbiddenMaskImage))
465  {
466  fiber.clear();
467  return resVal;
468  }
469  if (m_ForbiddenMaskImage->GetPixel(curNearestIndex) != 0)
470  {
471  fiber.clear();
472  return resVal;
473  }
474  }
475 
476  this->GetModelValue(curIndex,modelValue);
477  if (isZero(modelValue))
478  continueLoop = false;
479 
480  if (!this->CheckModelCompatibility(modelValue,threadId))
481  continueLoop = false;
482 
483  while (continueLoop)
484  {
485  if (fiber.size() <= 1)
486  {
487  oldDir = this->GetModelPrincipalDirection(modelValue,is2d,threadId);
488  for (unsigned int i = 0;i < 3;++i)
489  oldDir[i] *= -1.0;
490  }
491  else
492  {
493  for (unsigned int i = 0;i < 3;++i)
494  oldDir[i] = fiber[0][i] - fiber[1][i];
495  }
496 
497  // Do some stuff to progress one step
498  newDir = this->GetNextDirection(oldDir,modelValue,is2d,threadId);
499 
500  if (anima::ComputeOrientationAngle(oldDir, newDir) > m_MaxFiberAngle)
501  {
502  continueLoop = false;
503  break;
504  }
505 
506  if (anima::ComputeScalarProduct(oldDir, newDir) < 0)
507  anima::Revert(newDir,newDir);
508 
509  for (unsigned int i = 0;i < 3;++i)
510  curPoint[i] = fiber[0][i] + m_StepProgression * newDir[i];
511 
512  m_SeedingImage->TransformPhysicalPointToContinuousIndex(curPoint,curIndex);
513  m_SeedingImage->TransformPhysicalPointToIndex(curPoint,curNearestIndex);
514 
515  if (!m_CutMaskImage.IsNull())
516  {
517  if (!this->CheckIndexInImageBounds(curNearestIndex,m_CutMaskImage))
518  continueLoop = false;
519  else if (m_CutMaskImage->GetPixel(curNearestIndex) != 0)
520  continueLoop = false;
521  }
522 
523  if (!m_ForbiddenMaskImage.IsNull())
524  {
525  if (!this->CheckIndexInImageBounds(curNearestIndex,m_ForbiddenMaskImage))
526  {
527  fiber.clear();
528  return resVal;
529  }
530  if (m_ForbiddenMaskImage->GetPixel(curNearestIndex) != 0)
531  {
532  fiber.clear();
533  return resVal;
534  }
535  }
536 
537  if (!this->CheckIndexInImageBounds(curIndex))
538  continueLoop = false;
539  else
540  {
541  this->GetModelValue(curIndex,modelValue);
542 
543  if (isZero(modelValue))
544  continueLoop = false;
545  }
546 
547  if (!this->CheckModelCompatibility(modelValue,threadId))
548  continueLoop = false;
549 
550  // Now add point to fiber
551  if (continueLoop)
552  fiber.insert(fiber.begin(),1,curPoint);
553 
554  if (fiber.size() > m_MaxLengthFiber / m_StepProgression)
555  {
556  fiber.clear();
557  return resVal;
558  }
559  }
560  }
561 
562  return resVal;
563 }
564 
566 {
567  RegionType imageRegion = testImage->GetLargestPossibleRegion();
568  for (unsigned int i = 0;i < ImageBaseType::ImageDimension;++i)
569  {
570  if ((index[i] < imageRegion.GetIndex()[i]) ||
571  (index[i] >= imageRegion.GetIndex()[i] + imageRegion.GetSize()[i]))
572  return false;
573  }
574 
575  return true;
576 }
577 
579 {
580  for (unsigned int i = 0;i < value.GetSize();++i)
581  {
582  if (value[i] != 0)
583  return false;
584  }
585 
586  return true;
587 }
588 
589 } // end of namespace anima
void Revert(const VectorType &v, const unsigned int vSize, VectorType &resVec)
itk::ImageBase< ModelImageType::ImageDimension > ImageBaseType
std::vector< std::pair< FiberProgressType, FiberType > > FiberProcessVectorType
virtual void GetModelValue(ContinuousIndexType &index, VectorType &modelValue)=0
Computes value of model from data. May use SNR and previous model value to perform smart interpolatio...
virtual bool CheckModelCompatibility(VectorType &modelValue, itk::ThreadIdType threadId)=0
virtual PointType GetNextDirection(PointType &previousDirection, VectorType &modelValue, bool is2d, itk::ThreadIdType threadId)=0
void createVTKOutput(std::vector< std::vector< PointType > > &filteredFibers)
FiberProcessVectorType ComputeFiber(FiberType &fiber, FiberProgressType ways, itk::ThreadIdType threadId)
itk::ContinuousIndex< double, 3 > ContinuousIndexType
double ComputeOrientationAngle(const VectorType &v1, const VectorType &v2)
static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION ThreadTracker(void *arg)
bool CheckIndexInImageBounds(IndexType &index, ImageBaseType *testImage)
double ComputeScalarProduct(const VectorType &v1, const VectorType &v2, const unsigned int NDimension)
void ThreadTrack(unsigned int numThread, std::vector< FiberType > &resultFibers)
void ThreadedTrackComputer(unsigned int numThread, std::vector< FiberType > &resultFibers, unsigned int startSeedIndex, unsigned int endSeedIndex)
std::vector< FiberType > FilterOutputFibers(std::vector< FiberType > &fibers)
virtual PointType GetModelPrincipalDirection(VectorType &modelValue, bool is2d, itk::ThreadIdType threadId)=0