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
itkTractsToDWIImageFilter.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 #include "itkTractsToDWIImageFilter.h"
18 #include <boost/progress.hpp>
19 #include <vtkSmartPointer.h>
20 #include <vtkPolyData.h>
21 #include <vtkCellArray.h>
22 #include <vtkPoints.h>
23 #include <vtkPolyLine.h>
24 #include <itkImageRegionIteratorWithIndex.h>
25 #include <itkResampleImageFilter.h>
26 #include <itkNearestNeighborInterpolateImageFunction.h>
27 #include <itkBSplineInterpolateImageFunction.h>
28 #include <itkCastImageFilter.h>
29 #include <itkImageFileWriter.h>
30 #include <itkRescaleIntensityImageFilter.h>
31 #include <itkWindowedSincInterpolateImageFunction.h>
32 #include <itkResampleDwiImageFilter.h>
33 #include <itkKspaceImageFilter.h>
34 #include <itkDftImageFilter.h>
35 #include <itkAddImageFilter.h>
36 #include <itkConstantPadImageFilter.h>
37 #include <itkCropImageFilter.h>
38 #include <mitkAstroStickModel.h>
39 #include <vtkTransform.h>
40 #include <iostream>
41 #include <fstream>
42 #include <exception>
43 #include <itkImageDuplicator.h>
44 #include <itksys/SystemTools.hxx>
45 #include <mitkIOUtil.h>
46 #include <itkDiffusionTensor3DReconstructionImageFilter.h>
47 #include <itkDiffusionTensor3D.h>
48 #include <itkTractDensityImageFilter.h>
49 #include <boost/lexical_cast.hpp>
50 #include <itkTractsToVectorImageFilter.h>
51 #include <itkInvertIntensityImageFilter.h>
52 #include <itkShiftScaleImageFilter.h>
53 #include <itkExtractImageFilter.h>
54 #include <itkResampleDwiImageFilter.h>
55 #include <boost/algorithm/string/replace.hpp>
56 
57 namespace itk
58 {
59 
60  template< class PixelType >
61  TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter()
62  : m_FiberBundle(NULL)
63  , m_StatusText("")
64  , m_UseConstantRandSeed(false)
65  , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New())
66  {
67  m_RandGen->SetSeed();
68  }
69 
70  template< class PixelType >
71  TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter()
72  {
73 
74  }
75 
76  template< class PixelType >
77  TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::
78  SimulateKspaceAcquisition( std::vector< DoubleDwiType::Pointer >& images )
79  {
80  int numFiberCompartments = m_Parameters.m_FiberModelList.size();
81  // create slice object
82  ImageRegion<2> sliceRegion;
83  sliceRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]);
84  sliceRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]);
85  Vector< double, 2 > sliceSpacing;
86  sliceSpacing[0] = m_WorkingSpacing[0];
87  sliceSpacing[1] = m_WorkingSpacing[1];
88 
89  DoubleDwiType::PixelType nullPix; nullPix.SetSize(images.at(0)->GetVectorLength()); nullPix.Fill(0.0);
90  auto magnitudeDwiImage = DoubleDwiType::New();
91  magnitudeDwiImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
92  magnitudeDwiImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin );
93  magnitudeDwiImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
94  magnitudeDwiImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
95  magnitudeDwiImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
96  magnitudeDwiImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
97  magnitudeDwiImage->SetVectorLength( images.at(0)->GetVectorLength() );
98  magnitudeDwiImage->Allocate();
99  magnitudeDwiImage->FillBuffer(nullPix);
100 
101  m_PhaseImage = DoubleDwiType::New();
102  m_PhaseImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
103  m_PhaseImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin );
104  m_PhaseImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
105  m_PhaseImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
106  m_PhaseImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
107  m_PhaseImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
108  m_PhaseImage->SetVectorLength( images.at(0)->GetVectorLength() );
109  m_PhaseImage->Allocate();
110  m_PhaseImage->FillBuffer(nullPix);
111 
112  m_KspaceImage = DoubleDwiType::New();
113  m_KspaceImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
114  m_KspaceImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin );
115  m_KspaceImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
116  m_KspaceImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
117  m_KspaceImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
118  m_KspaceImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
119  m_KspaceImage->SetVectorLength( m_Parameters.m_SignalGen.m_NumberOfCoils );
120  m_KspaceImage->Allocate();
121  m_KspaceImage->FillBuffer(nullPix);
122 
123  std::vector< unsigned int > spikeVolume;
124  for (unsigned int i=0; i<m_Parameters.m_SignalGen.m_Spikes; i++)
125  spikeVolume.push_back(m_RandGen->GetIntegerVariate()%(images.at(0)->GetVectorLength()));
126  std::sort (spikeVolume.begin(), spikeVolume.end());
127  std::reverse (spikeVolume.begin(), spikeVolume.end());
128 
129  // calculate coil positions
130  double a = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)*m_Parameters.m_SignalGen.m_ImageSpacing[0];
131  double b = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)*m_Parameters.m_SignalGen.m_ImageSpacing[1];
132  double c = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(2)*m_Parameters.m_SignalGen.m_ImageSpacing[2];
133  double diagonal = sqrt(a*a+b*b)/1000; // image diagonal in m
134 
135  m_CoilPointset = mitk::PointSet::New();
136  std::vector< itk::Vector<double, 3> > coilPositions;
137  itk::Vector<double, 3> pos; pos.Fill(0.0); pos[1] = -diagonal/2;
138  itk::Vector<double, 3> center;
139  center[0] = a/2-m_Parameters.m_SignalGen.m_ImageSpacing[0]/2;
140  center[1] = b/2-m_Parameters.m_SignalGen.m_ImageSpacing[2]/2;
141  center[2] = c/2-m_Parameters.m_SignalGen.m_ImageSpacing[1]/2;
142  for (int c=0; c<m_Parameters.m_SignalGen.m_NumberOfCoils; c++)
143  {
144  coilPositions.push_back(pos);
145  m_CoilPointset->InsertPoint(c, pos*1000 + m_Parameters.m_SignalGen.m_ImageOrigin.GetVectorFromOrigin() + center );
146 
147  double rz = 360.0/m_Parameters.m_SignalGen.m_NumberOfCoils * M_PI/180;
148  vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity();
149  rotZ[0][0] = cos(rz);
150  rotZ[1][1] = rotZ[0][0];
151  rotZ[0][1] = -sin(rz);
152  rotZ[1][0] = -rotZ[0][1];
153 
154  pos.SetVnlVector(rotZ*pos.GetVnlVector());
155  }
156 
157  PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false);
158  PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false);
159  unsigned long lastTick = 0;
160 
161  boost::progress_display disp(images.at(0)->GetVectorLength()*images.at(0)->GetLargestPossibleRegion().GetSize(2));
162  for (unsigned int g=0; g<images.at(0)->GetVectorLength(); g++)
163  {
164  std::vector< unsigned int > spikeSlice;
165  while (!spikeVolume.empty() && spikeVolume.back()==g)
166  {
167  spikeSlice.push_back(m_RandGen->GetIntegerVariate()%images.at(0)->GetLargestPossibleRegion().GetSize(2));
168  spikeVolume.pop_back();
169  }
170  std::sort (spikeSlice.begin(), spikeSlice.end());
171  std::reverse (spikeSlice.begin(), spikeSlice.end());
172 
173  for (unsigned int z=0; z<images.at(0)->GetLargestPossibleRegion().GetSize(2); z++)
174  {
175  std::vector< SliceType::Pointer > compartmentSlices;
176  std::vector< double > t2Vector;
177  std::vector< double > t1Vector;
178 
179  for (unsigned int i=0; i<images.size(); i++)
180  {
181  DiffusionSignalModel<double>* signalModel;
182  if (i<numFiberCompartments)
183  signalModel = m_Parameters.m_FiberModelList.at(i);
184  else
185  signalModel = m_Parameters.m_NonFiberModelList.at(i-numFiberCompartments);
186 
187  auto slice = SliceType::New();
188  slice->SetLargestPossibleRegion( sliceRegion );
189  slice->SetBufferedRegion( sliceRegion );
190  slice->SetRequestedRegion( sliceRegion );
191  slice->SetSpacing(sliceSpacing);
192  slice->Allocate();
193  slice->FillBuffer(0.0);
194 
195  // extract slice from channel g
196  for (unsigned int y=0; y<images.at(0)->GetLargestPossibleRegion().GetSize(1); y++)
197  for (unsigned int x=0; x<images.at(0)->GetLargestPossibleRegion().GetSize(0); x++)
198  {
199  SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y;
200  DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
201 
202  slice->SetPixel(index2D, images.at(i)->GetPixel(index3D)[g]);
203  }
204 
205  compartmentSlices.push_back(slice);
206  t2Vector.push_back(signalModel->GetT2());
207  t1Vector.push_back(signalModel->GetT1());
208  }
209 
210  int numSpikes = 0;
211  while (!spikeSlice.empty() && spikeSlice.back()==z)
212  {
213  numSpikes++;
214  spikeSlice.pop_back();
215  }
216  int spikeCoil = m_RandGen->GetIntegerVariate()%m_Parameters.m_SignalGen.m_NumberOfCoils;
217 
218  if (this->GetAbortGenerateData())
219  return NULL;
220 
221 #pragma omp parallel for
222  for (int c=0; c<m_Parameters.m_SignalGen.m_NumberOfCoils; c++)
223  {
224  // create k-sapce (inverse fourier transform slices)
225  auto idft = itk::KspaceImageFilter< SliceType::PixelType >::New();
226  idft->SetCompartmentImages(compartmentSlices);
227  idft->SetT2(t2Vector);
228  idft->SetT1(t1Vector);
229  idft->SetUseConstantRandSeed(m_UseConstantRandSeed);
230  idft->SetParameters(&m_Parameters);
231  idft->SetZ((double)z-(double)( images.at(0)->GetLargestPossibleRegion().GetSize(2)
232  -images.at(0)->GetLargestPossibleRegion().GetSize(2)%2 ) / 2.0);
233  idft->SetZidx(z);
234  idft->SetCoilPosition(coilPositions.at(c));
235  idft->SetFiberBundle(m_FiberBundleWorkingCopy);
236  idft->SetTranslation(m_Translations.at(g));
237  idft->SetRotation(m_Rotations.at(g));
238  idft->SetDiffusionGradientDirection(m_Parameters.m_SignalGen.GetGradientDirection(g));
239  if (c==spikeCoil)
240  idft->SetSpikesPerSlice(numSpikes);
241  idft->Update();
242 
243 #pragma omp critical
244  if (c==spikeCoil && numSpikes>0)
245  {
246  m_SpikeLog += "Volume " + boost::lexical_cast<std::string>(g) + " Coil " + boost::lexical_cast<std::string>(c) + "\n";
247  m_SpikeLog += idft->GetSpikeLog();
248  }
249 
250  ComplexSliceType::Pointer fSlice;
251  fSlice = idft->GetOutput();
252 
253  // fourier transform slice
254  ComplexSliceType::Pointer newSlice;
255  auto dft = itk::DftImageFilter< SliceType::PixelType >::New();
256  dft->SetInput(fSlice);
257  dft->SetParameters(m_Parameters);
258  dft->Update();
259  newSlice = dft->GetOutput();
260 
261  // put slice back into channel g
262  for (unsigned int y=0; y<fSlice->GetLargestPossibleRegion().GetSize(1); y++)
263  for (unsigned int x=0; x<fSlice->GetLargestPossibleRegion().GetSize(0); x++)
264  {
265  DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
266  ComplexSliceType::IndexType index2D; index2D[0]=x; index2D[1]=y;
267 
268  ComplexSliceType::PixelType cPix = newSlice->GetPixel(index2D);
269  double magn = sqrt(cPix.real()*cPix.real()+cPix.imag()*cPix.imag());
270  double phase = 0;
271  if (cPix.real()!=0)
272  phase = atan( cPix.imag()/cPix.real() );
273 
274 
275  DoubleDwiType::PixelType dwiPix = magnitudeDwiImage->GetPixel(index3D);
276  DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D);
277 
278  if (m_Parameters.m_SignalGen.m_NumberOfCoils>1)
279  {
280  dwiPix[g] += magn*magn;
281  phasePix[g] += phase*phase;
282  }
283  else
284  {
285  dwiPix[g] = magn;
286  phasePix[g] = phase;
287  }
288 
289 #pragma omp critical
290  {
291  magnitudeDwiImage->SetPixel(index3D, dwiPix);
292  m_PhaseImage->SetPixel(index3D, phasePix);
293 
294  // k-space image
295  if (g==0)
296  {
297  DoubleDwiType::PixelType kspacePix = m_KspaceImage->GetPixel(index3D);
298  kspacePix[c] = idft->GetKSpaceImage()->GetPixel(index2D);
299  m_KspaceImage->SetPixel(index3D, kspacePix);
300  }
301  }
302  }
303  }
304 
305  if (m_Parameters.m_SignalGen.m_NumberOfCoils>1)
306  {
307 #ifdef WIN32
308 #pragma omp parallel for
309 #else
310 #pragma omp parallel for collapse(2)
311 #endif
312  for (int y=0; y<magnitudeDwiImage->GetLargestPossibleRegion().GetSize(1); y++)
313  for (int x=0; x<magnitudeDwiImage->GetLargestPossibleRegion().GetSize(0); x++)
314  {
315  DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z;
316  DoubleDwiType::PixelType magPix = magnitudeDwiImage->GetPixel(index3D);
317  magPix[g] = sqrt(magPix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils);
318 
319  DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D);
320  phasePix[g] = sqrt(phasePix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils);
321 
322 #pragma omp critical
323  {
324  magnitudeDwiImage->SetPixel(index3D, magPix);
325  m_PhaseImage->SetPixel(index3D, phasePix);
326  }
327  }
328  }
329 
330  ++disp;
331  unsigned long newTick = 50*disp.count()/disp.expected_count();
332  for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
333  PrintToLog("*", false, false, false);
334  lastTick = newTick;
335  }
336  }
337  PrintToLog("\n", false);
338  return magnitudeDwiImage;
339  }
340 
341  template< class PixelType >
342  TractsToDWIImageFilter< PixelType >::ItkDoubleImgType::Pointer TractsToDWIImageFilter< PixelType >::
343  NormalizeInsideMask(ItkDoubleImgType::Pointer image)
344  {
345  double max = itk::NumericTraits< double >::min();
346  double min = itk::NumericTraits< double >::max();
347 
348  itk::ImageRegionIterator< ItkDoubleImgType > it(image, image->GetLargestPossibleRegion());
349  while(!it.IsAtEnd())
350  {
351  if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && m_Parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0)
352  {
353  it.Set(0.0);
354  ++it;
355  continue;
356  }
357  // if (it.Get()>900)
358  // it.Set(900);
359 
360  if (it.Get()>max)
361  max = it.Get();
362  if (it.Get()<min)
363  min = it.Get();
364  ++it;
365  }
366  auto scaler = itk::ShiftScaleImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New();
367  scaler->SetInput(image);
368  scaler->SetShift(-min);
369  scaler->SetScale(1.0/(max-min));
370  scaler->Update();
371  return scaler->GetOutput();
372  }
373 
374  template< class PixelType >
375  void TractsToDWIImageFilter< PixelType >::CheckVolumeFractionImages()
376  {
377  m_UseRelativeNonFiberVolumeFractions = false;
378 
379  // check for fiber volume fraction maps
380  int fibVolImages = 0;
381  for (int i=0; i<m_Parameters.m_FiberModelList.size(); i++)
382  {
383  if (m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage().IsNotNull())
384  {
385  PrintToLog("Using volume fraction map for fiber compartment " + boost::lexical_cast<std::string>(i+1));
386  fibVolImages++;
387  }
388  }
389 
390  // check for non-fiber volume fraction maps
391  int nonfibVolImages = 0;
392  for (int i=0; i<m_Parameters.m_NonFiberModelList.size(); i++)
393  {
394  if (m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage().IsNotNull())
395  {
396  PrintToLog("Using volume fraction map for non-fiber compartment " + boost::lexical_cast<std::string>(i+1));
397  nonfibVolImages++;
398  }
399  }
400 
401  // not all fiber compartments are using volume fraction maps
402  // --> non-fiber volume fractions are assumed to be relative to the
403  // non-fiber volume and not absolute voxel-volume fractions.
404  // this means if two non-fiber compartments are used but only one of them
405  // has an associated volume fraction map, the repesctive other volume fraction map
406  // can be determined as inverse (1-val) of the present volume fraction map-
407  if ( fibVolImages<m_Parameters.m_FiberModelList.size() && nonfibVolImages==1 && m_Parameters.m_NonFiberModelList.size()==2)
408  {
409  PrintToLog("Calculating missing non-fiber volume fraction image by inverting the other.\n"
410  "Assuming non-fiber volume fraction images to contain values relative to the"
411  " remaining non-fiber volume, not absolute values.");
412 
413  auto inverter = itk::InvertIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New();
414  inverter->SetMaximum(1.0);
415  if ( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNull()
416  && m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNotNull() )
417  {
418  // m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(
419  // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ) );
420  inverter->SetInput( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() );
421  inverter->Update();
422  m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(inverter->GetOutput());
423  }
424  else if ( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNull()
425  && m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNotNull() )
426  {
427  // m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(
428  // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ) );
429  inverter->SetInput( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() );
430  inverter->Update();
431  m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(inverter->GetOutput());
432  }
433  else
434  {
435  itkExceptionMacro("Something went wrong in automatically calculating the missing non-fiber volume fraction image!"
436  " Did you use two non fiber compartments but only one volume fraction image?"
437  " Then it should work and this error is really strange.");
438  }
439  m_UseRelativeNonFiberVolumeFractions = true;
440 
441  nonfibVolImages++;
442  }
443 
444  // Up to two fiber compartments are allowed without volume fraction maps since the volume fractions can then be determined automatically
445  if (m_Parameters.m_FiberModelList.size()>2
446  && fibVolImages!=m_Parameters.m_FiberModelList.size())
447  itkExceptionMacro("More than two fiber compartment selected but no corresponding volume fraction maps set!");
448 
449  // One non-fiber compartment is allowed without volume fraction map since the volume fraction can then be determined automatically
450  if (m_Parameters.m_NonFiberModelList.size()>1
451  && nonfibVolImages!=m_Parameters.m_NonFiberModelList.size())
452  itkExceptionMacro("More than one non-fiber compartment selected but no volume fraction maps set!");
453 
454  if (fibVolImages<m_Parameters.m_FiberModelList.size() && nonfibVolImages>0)
455  {
456  PrintToLog("Not all fiber compartments are using an associated volume fraction image.\n"
457  "Assuming non-fiber volume fraction images to contain values relative to the"
458  " remaining non-fiber volume, not absolute values.");
459  m_UseRelativeNonFiberVolumeFractions = true;
460 
461  // itk::ImageFileWriter<ItkDoubleImgType>::Pointer wr = itk::ImageFileWriter<ItkDoubleImgType>::New();
462  // wr->SetInput(m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage());
463  // wr->SetFileName("/local/volumefraction.nrrd");
464  // wr->Update();
465  }
466 
467  // initialize the images that store the output volume fraction of each compartment
468  m_VolumeFractions.clear();
469  for (int i=0; i<m_Parameters.m_FiberModelList.size()+m_Parameters.m_NonFiberModelList.size(); i++)
470  {
471  auto doubleImg = ItkDoubleImgType::New();
472  doubleImg->SetSpacing( m_WorkingSpacing );
473  doubleImg->SetOrigin( m_WorkingOrigin );
474  doubleImg->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
475  doubleImg->SetLargestPossibleRegion( m_WorkingImageRegion );
476  doubleImg->SetBufferedRegion( m_WorkingImageRegion );
477  doubleImg->SetRequestedRegion( m_WorkingImageRegion );
478  doubleImg->Allocate();
479  doubleImg->FillBuffer(0);
480  m_VolumeFractions.push_back(doubleImg);
481  }
482  }
483 
484  template< class PixelType >
485  void TractsToDWIImageFilter< PixelType >::InitializeData()
486  {
487  m_Rotations.clear();
488  m_Translations.clear();
489  m_MotionLog = "";
490  m_SpikeLog = "";
491 
492  // initialize output dwi image
493  m_Parameters.m_SignalGen.m_CroppedRegion = m_Parameters.m_SignalGen.m_ImageRegion;
494  m_Parameters.m_SignalGen.m_CroppedRegion.SetSize( 1, m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1)
495  *m_Parameters.m_SignalGen.m_CroppingFactor);
496 
497  itk::Point<double,3> shiftedOrigin = m_Parameters.m_SignalGen.m_ImageOrigin;
498  shiftedOrigin[1] += (m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)
499  -m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1))*m_Parameters.m_SignalGen.m_ImageSpacing[1]/2;
500 
501  m_OutputImage = OutputImageType::New();
502  m_OutputImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing );
503  m_OutputImage->SetOrigin( shiftedOrigin );
504  m_OutputImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
505  m_OutputImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
506  m_OutputImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
507  m_OutputImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion );
508  m_OutputImage->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() );
509  m_OutputImage->Allocate();
510  typename OutputImageType::PixelType temp;
511  temp.SetSize(m_Parameters.m_SignalGen.GetNumVolumes());
512  temp.Fill(0.0);
513  m_OutputImage->FillBuffer(temp);
514 
515  // Apply in-plane upsampling for Gibbs ringing artifact
516  double upsampling = 1;
517  if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging)
518  upsampling = 2;
519  m_WorkingSpacing = m_Parameters.m_SignalGen.m_ImageSpacing;
520  m_WorkingSpacing[0] /= upsampling;
521  m_WorkingSpacing[1] /= upsampling;
522  m_WorkingImageRegion = m_Parameters.m_SignalGen.m_ImageRegion;
523  m_WorkingImageRegion.SetSize(0, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[0]*upsampling);
524  m_WorkingImageRegion.SetSize(1, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[1]*upsampling);
525  m_WorkingOrigin = m_Parameters.m_SignalGen.m_ImageOrigin;
526  m_WorkingOrigin[0] -= m_Parameters.m_SignalGen.m_ImageSpacing[0]/2;
527  m_WorkingOrigin[0] += m_WorkingSpacing[0]/2;
528  m_WorkingOrigin[1] -= m_Parameters.m_SignalGen.m_ImageSpacing[1]/2;
529  m_WorkingOrigin[1] += m_WorkingSpacing[1]/2;
530  m_WorkingOrigin[2] -= m_Parameters.m_SignalGen.m_ImageSpacing[2]/2;
531  m_WorkingOrigin[2] += m_WorkingSpacing[2]/2;
532  m_VoxelVolume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2];
533 
534  // generate double images to store the individual compartment signals
535  m_CompartmentImages.clear();
536  int numFiberCompartments = m_Parameters.m_FiberModelList.size();
537  int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size();
538  for (int i=0; i<numFiberCompartments+numNonFiberCompartments; i++)
539  {
540  auto doubleDwi = DoubleDwiType::New();
541  doubleDwi->SetSpacing( m_WorkingSpacing );
542  doubleDwi->SetOrigin( m_WorkingOrigin );
543  doubleDwi->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
544  doubleDwi->SetLargestPossibleRegion( m_WorkingImageRegion );
545  doubleDwi->SetBufferedRegion( m_WorkingImageRegion );
546  doubleDwi->SetRequestedRegion( m_WorkingImageRegion );
547  doubleDwi->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() );
548  doubleDwi->Allocate();
549  DoubleDwiType::PixelType pix;
550  pix.SetSize(m_Parameters.m_SignalGen.GetNumVolumes());
551  pix.Fill(0.0);
552  doubleDwi->FillBuffer(pix);
553  m_CompartmentImages.push_back(doubleDwi);
554  }
555 
556  if (m_FiberBundle.IsNull() && m_InputImage.IsNotNull())
557  {
558  m_CompartmentImages.clear();
559 
560  m_Parameters.m_SignalGen.m_DoAddMotion = false;
561  m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false;
562 
563  PrintToLog("Simulating acquisition for input diffusion-weighted image.", false);
564  auto caster = itk::CastImageFilter< OutputImageType, DoubleDwiType >::New();
565  caster->SetInput(m_InputImage);
566  caster->Update();
567 
568  if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging)
569  {
570  PrintToLog("Upsampling input diffusion-weighted image for Gibbs ringing simulation.", false);
571 
572  auto resampler = itk::ResampleDwiImageFilter< double >::New();
573  resampler->SetInput(caster->GetOutput());
574  itk::Vector< double, 3 > samplingFactor;
575  samplingFactor[0] = upsampling;
576  samplingFactor[1] = upsampling;
577  samplingFactor[2] = 1;
578  resampler->SetSamplingFactor(samplingFactor);
579  resampler->SetInterpolation(itk::ResampleDwiImageFilter< double >::Interpolate_WindowedSinc);
580  resampler->Update();
581  m_CompartmentImages.push_back(resampler->GetOutput());
582  }
583  else
584  m_CompartmentImages.push_back(caster->GetOutput());
585 
586  for (unsigned int g=0; g<m_Parameters.m_SignalGen.GetNumVolumes(); g++)
587  {
588  m_Rotation.Fill(0.0);
589  m_Translation.Fill(0.0);
590  m_Rotations.push_back(m_Rotation);
591  m_Translations.push_back(m_Translation);
592  }
593  }
594 
595  // resample mask image and frequency map to fit upsampled geometry
596  if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging)
597  {
598  if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull())
599  {
600  // rescale mask image (otherwise there are problems with the resampling)
601  auto rescaler = itk::RescaleIntensityImageFilter<ItkUcharImgType,ItkUcharImgType>::New();
602  rescaler->SetInput(0,m_Parameters.m_SignalGen.m_MaskImage);
603  rescaler->SetOutputMaximum(100);
604  rescaler->SetOutputMinimum(0);
605  rescaler->Update();
606 
607  // resample mask image
608  auto resampler = itk::ResampleImageFilter<ItkUcharImgType, ItkUcharImgType>::New();
609  resampler->SetInput(rescaler->GetOutput());
610  resampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage);
611  resampler->SetSize(m_WorkingImageRegion.GetSize());
612  resampler->SetOutputSpacing(m_WorkingSpacing);
613  resampler->SetOutputOrigin(m_WorkingOrigin);
614  auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction<ItkUcharImgType, double>::New();
615  resampler->SetInterpolator(nn_interpolator);
616  resampler->Update();
617  m_Parameters.m_SignalGen.m_MaskImage = resampler->GetOutput();
618  }
619  // resample frequency map
620  if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull())
621  {
622  auto resampler = itk::ResampleImageFilter<ItkDoubleImgType, ItkDoubleImgType>::New();
623  resampler->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap);
624  resampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_FrequencyMap);
625  resampler->SetSize(m_WorkingImageRegion.GetSize());
626  resampler->SetOutputSpacing(m_WorkingSpacing);
627  resampler->SetOutputOrigin(m_WorkingOrigin);
628  auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction<ItkDoubleImgType, double>::New();
629  resampler->SetInterpolator(nn_interpolator);
630  resampler->Update();
631  m_Parameters.m_SignalGen.m_FrequencyMap = resampler->GetOutput();
632  }
633  }
634 
635  m_MaskImageSet = true;
636  if (m_Parameters.m_SignalGen.m_MaskImage.IsNull())
637  {
638  // no input tissue mask is set -> create default
639  PrintToLog("No tissue mask set", false);
640  m_Parameters.m_SignalGen.m_MaskImage = ItkUcharImgType::New();
641  m_Parameters.m_SignalGen.m_MaskImage->SetSpacing( m_WorkingSpacing );
642  m_Parameters.m_SignalGen.m_MaskImage->SetOrigin( m_WorkingOrigin );
643  m_Parameters.m_SignalGen.m_MaskImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
644  m_Parameters.m_SignalGen.m_MaskImage->SetLargestPossibleRegion( m_WorkingImageRegion );
645  m_Parameters.m_SignalGen.m_MaskImage->SetBufferedRegion( m_WorkingImageRegion );
646  m_Parameters.m_SignalGen.m_MaskImage->SetRequestedRegion( m_WorkingImageRegion );
647  m_Parameters.m_SignalGen.m_MaskImage->Allocate();
648  m_Parameters.m_SignalGen.m_MaskImage->FillBuffer(100);
649  m_MaskImageSet = false;
650  }
651  else
652  {
653  if (m_Parameters.m_SignalGen.m_MaskImage->GetLargestPossibleRegion()!=m_WorkingImageRegion)
654  {
655  itkExceptionMacro("Mask image and specified DWI geometry are not matching!");
656  }
657  PrintToLog("Using tissue mask", false);
658  }
659 
660  if (m_Parameters.m_SignalGen.m_DoAddMotion)
661  {
662  if (m_Parameters.m_SignalGen.m_DoRandomizeMotion)
663  {
664  PrintToLog("Random motion artifacts:", false);
665  PrintToLog("Maximum rotation: +/-" + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Rotation) + "°", false);
666  PrintToLog("Maximum translation: +/-" + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Translation) + "mm", false);
667  }
668  else
669  {
670  PrintToLog("Linear motion artifacts:", false);
671  PrintToLog("Maximum rotation: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Rotation) + "°", false);
672  PrintToLog("Maximum translation: " + boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_Translation) + "mm", false);
673  }
674  }
675  if ( m_Parameters.m_SignalGen.m_MotionVolumes.empty() )
676  {
677  // no motion in first volume
678  m_Parameters.m_SignalGen.m_MotionVolumes.push_back(false);
679 
680  // motion in all other volumes
681  while ( m_Parameters.m_SignalGen.m_MotionVolumes.size() < m_Parameters.m_SignalGen.GetNumVolumes() )
682  {
683  m_Parameters.m_SignalGen.m_MotionVolumes.push_back(true);
684  }
685  }
686  // we need to know for every volume if there is motion. if this information is missing, then set corresponding fal to false
687  while ( m_Parameters.m_SignalGen.m_MotionVolumes.size()<m_Parameters.m_SignalGen.GetNumVolumes() )
688  {
689  m_Parameters.m_SignalGen.m_MotionVolumes.push_back(false);
690  }
691 
692  m_NumMotionVolumes = 0;
693  for (int i=0; i<m_Parameters.m_SignalGen.GetNumVolumes(); i++)
694  {
695  if (m_Parameters.m_SignalGen.m_MotionVolumes[i])
696  {
697  m_NumMotionVolumes++;
698  }
699  }
700  m_MotionCounter = 0;
701 
702  // creat image to hold transformed mask (motion artifact)
703  m_TransformedMaskImage = ItkUcharImgType::New();
704  auto duplicator = itk::ImageDuplicator<ItkUcharImgType>::New();
705  duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage);
706  duplicator->Update();
707  m_TransformedMaskImage = duplicator->GetOutput();
708 
709  // second upsampling needed for motion artifacts
710  ImageRegion<3> upsampledImageRegion = m_WorkingImageRegion;
711  DoubleVectorType upsampledSpacing = m_WorkingSpacing;
712  upsampledSpacing[0] /= 4;
713  upsampledSpacing[1] /= 4;
714  upsampledSpacing[2] /= 4;
715  upsampledImageRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]*4);
716  upsampledImageRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]*4);
717  upsampledImageRegion.SetSize(2, m_WorkingImageRegion.GetSize()[2]*4);
718  itk::Point<double,3> upsampledOrigin = m_WorkingOrigin;
719  upsampledOrigin[0] -= m_WorkingSpacing[0]/2;
720  upsampledOrigin[0] += upsampledSpacing[0]/2;
721  upsampledOrigin[1] -= m_WorkingSpacing[1]/2;
722  upsampledOrigin[1] += upsampledSpacing[1]/2;
723  upsampledOrigin[2] -= m_WorkingSpacing[2]/2;
724  upsampledOrigin[2] += upsampledSpacing[2]/2;
725 
726  m_UpsampledMaskImage = ItkUcharImgType::New();
727  auto upsampler = itk::ResampleImageFilter<ItkUcharImgType, ItkUcharImgType>::New();
728  upsampler->SetInput(m_Parameters.m_SignalGen.m_MaskImage);
729  upsampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage);
730  upsampler->SetSize(upsampledImageRegion.GetSize());
731  upsampler->SetOutputSpacing(upsampledSpacing);
732  upsampler->SetOutputOrigin(upsampledOrigin);
733 
734  auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction<ItkUcharImgType, double>::New();
735  upsampler->SetInterpolator(nn_interpolator);
736  upsampler->Update();
737  m_UpsampledMaskImage = upsampler->GetOutput();
738  }
739 
740  template< class PixelType >
741  void TractsToDWIImageFilter< PixelType >::InitializeFiberData()
742  {
743  // resample fiber bundle for sufficient voxel coverage
744  PrintToLog("Resampling fibers ...");
745  m_SegmentVolume = 0.0001;
746  float minSpacing = 1;
747  if( m_WorkingSpacing[0]<m_WorkingSpacing[1]
748  && m_WorkingSpacing[0]<m_WorkingSpacing[2])
749  {
750  minSpacing = m_WorkingSpacing[0];
751  }
752  else if (m_WorkingSpacing[1] < m_WorkingSpacing[2])
753  {
754  minSpacing = m_WorkingSpacing[1];
755  }
756  else
757  {
758  minSpacing = m_WorkingSpacing[2];
759  }
760 
761  // working copy is needed because we need to resample the fibers but do not want to change the original bundle
762  m_FiberBundleWorkingCopy = m_FiberBundle->GetDeepCopy();
763  double volumeAccuracy = 10;
764  m_FiberBundleWorkingCopy->ResampleSpline(minSpacing/volumeAccuracy);
765  m_mmRadius = m_Parameters.m_SignalGen.m_AxonRadius/1000;
766 
767  if (m_mmRadius>0) { m_SegmentVolume = M_PI*m_mmRadius*m_mmRadius*minSpacing/volumeAccuracy; }
768  // a second fiber bundle is needed to store the transformed version of the m_FiberBundleWorkingCopy
769  m_FiberBundleTransformed = m_FiberBundleWorkingCopy;
770  }
771 
772 
773  template< class PixelType >
774  bool TractsToDWIImageFilter< PixelType >::PrepareLogFile()
775  {
776  assert( ! m_Logfile.is_open() );
777 
778  std::string filePath;
779  std::string fileName;
780 
781  // Get directory name:
782  if (m_Parameters.m_Misc.m_OutputPath.size() > 0)
783  {
784  filePath = m_Parameters.m_Misc.m_OutputPath;
785  if( *(--(filePath.cend())) != '/')
786  {
787  filePath.push_back('/');
788  }
789  }
790  else
791  {
792  filePath = mitk::IOUtil::GetTempPath() + '/';
793  }
794  // check if directory exists, else use /tmp/:
795  if( itksys::SystemTools::FileIsDirectory( filePath ) )
796  {
797  while( *(--(filePath.cend())) == '/')
798  {
799  filePath.pop_back();
800  }
801  filePath = filePath + '/';
802  }
803  else
804  {
805  filePath = mitk::IOUtil::GetTempPath() + '/';
806  }
807 
808  // Get file name:
809  if( ! m_Parameters.m_Misc.m_ResultNode->GetName().empty() )
810  {
811  fileName = m_Parameters.m_Misc.m_ResultNode->GetName();
812  }
813  else
814  {
815  fileName = "";
816  }
817 
818  if( ! m_Parameters.m_Misc.m_OutputPrefix.empty() )
819  {
820  fileName = m_Parameters.m_Misc.m_OutputPrefix + fileName;
821  }
822  else
823  {
824  fileName = "fiberfox";
825  }
826 
827  // check if file already exists and DO NOT overwrite existing files:
828  std::string NameTest = fileName;
829  int c = 0;
830  while( itksys::SystemTools::FileExists( filePath + '/' + fileName + ".log" )
831  && c <= std::numeric_limits<int>::max() )
832  {
833  fileName = NameTest + "_" + boost::lexical_cast<std::string>(c);
834  ++c;
835  }
836 
837  try
838  {
839  m_Logfile.open( ( filePath + '/' + fileName + ".log" ).c_str() );
840  }
841  catch (const std::ios_base::failure &fail)
842  {
843  MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Exception " << fail.what()
844  << " while trying to open file" << filePath << '/' << fileName << ".log";
845  return false;
846  }
847 
848  if ( m_Logfile.is_open() )
849  {
850  PrintToLog( "Logfile: " + filePath + '/' + fileName + ".log", false );
851  return true;
852  }
853  else
854  {
855  m_StatusText += "Logfile could not be opened!\n";
856  MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Logfile could not be opened!";
857  return false;
858  }
859  }
860 
861 
862  template< class PixelType >
863  void TractsToDWIImageFilter< PixelType >::GenerateData()
864  {
865  // prepare logfile
866  if ( ! PrepareLogFile() )
867  {
868  this->SetAbortGenerateData( true );
869  return;
870  }
871 
872  m_TimeProbe.Start();
873 
874  // check input data
875  if (m_FiberBundle.IsNull() && m_InputImage.IsNull())
876  itkExceptionMacro("Input fiber bundle and input diffusion-weighted image is NULL!");
877 
878  if (m_Parameters.m_FiberModelList.empty() && m_InputImage.IsNull())
879  itkExceptionMacro("No diffusion model for fiber compartments defined and input diffusion-weighted"
880  " image is NULL! At least one fiber compartment is necessary to simulate diffusion.");
881 
882  if (m_Parameters.m_NonFiberModelList.empty() && m_InputImage.IsNull())
883  itkExceptionMacro("No diffusion model for non-fiber compartments defined and input diffusion-weighted"
884  " image is NULL! At least one non-fiber compartment is necessary to simulate diffusion.");
885 
886  int baselineIndex = m_Parameters.m_SignalGen.GetFirstBaselineIndex();
887  if (baselineIndex<0) { itkExceptionMacro("No baseline index found!"); }
888 
889  if (!m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition) // No upsampling of input image needed if no k-space simulation is performed
890  {
891  m_Parameters.m_SignalGen.m_DoAddGibbsRinging = false;
892  }
893 
894  if (m_UseConstantRandSeed) // always generate the same random numbers?
895  { m_RandGen->SetSeed(0); }
896  else { m_RandGen->SetSeed(); }
897 
898  InitializeData();
899  if ( m_FiberBundle.IsNotNull() ) // if no fiber bundle is found, we directly proceed to the k-space acquisition simulation
900  {
901  CheckVolumeFractionImages();
902  InitializeFiberData();
903 
904  int numFiberCompartments = m_Parameters.m_FiberModelList.size();
905  int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size();
906 
907  double maxVolume = 0;
908  unsigned long lastTick = 0;
909  int signalModelSeed = m_RandGen->GetIntegerVariate();
910 
911  PrintToLog("\n", false, false);
912  PrintToLog("Generating " + boost::lexical_cast<std::string>(numFiberCompartments+numNonFiberCompartments)
913  + "-compartment diffusion-weighted signal.");
914 
915  int numFibers = m_FiberBundleWorkingCopy->GetNumFibers();
916  boost::progress_display disp(numFibers*m_Parameters.m_SignalGen.GetNumVolumes());
917 
918  PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false);
919  PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false);
920 
921  for (unsigned int g=0; g<m_Parameters.m_SignalGen.GetNumVolumes(); g++)
922  {
923  // move fibers
924  SimulateMotion(g);
925 
926  // Set signal model random generator seeds to get same configuration in each voxel
927  for (int i=0; i<m_Parameters.m_FiberModelList.size(); i++)
928  m_Parameters.m_FiberModelList.at(i)->SetSeed(signalModelSeed);
929  for (int i=0; i<m_Parameters.m_NonFiberModelList.size(); i++)
930  m_Parameters.m_NonFiberModelList.at(i)->SetSeed(signalModelSeed);
931 
932  // storing voxel-wise intra-axonal volume in mm³
933  auto intraAxonalVolumeImage = ItkDoubleImgType::New();
934  intraAxonalVolumeImage->SetSpacing( m_WorkingSpacing );
935  intraAxonalVolumeImage->SetOrigin( m_WorkingOrigin );
936  intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection );
937  intraAxonalVolumeImage->SetLargestPossibleRegion( m_WorkingImageRegion );
938  intraAxonalVolumeImage->SetBufferedRegion( m_WorkingImageRegion );
939  intraAxonalVolumeImage->SetRequestedRegion( m_WorkingImageRegion );
940  intraAxonalVolumeImage->Allocate();
941  intraAxonalVolumeImage->FillBuffer(0);
942  maxVolume = 0;
943 
944  vtkPolyData* fiberPolyData = m_FiberBundleTransformed->GetFiberPolyData();
945  // generate fiber signal (if there are any fiber models present)
946  if (!m_Parameters.m_FiberModelList.empty())
947  for( int i=0; i<numFibers; i++ )
948  {
949  float fiberWeight = m_FiberBundleTransformed->GetFiberWeight(i);
950  vtkCell* cell = fiberPolyData->GetCell(i);
951  int numPoints = cell->GetNumberOfPoints();
952  vtkPoints* points = cell->GetPoints();
953 
954  if (numPoints<2)
955  continue;
956 
957  for( int j=0; j<numPoints; j++)
958  {
959  if (this->GetAbortGenerateData())
960  {
961  PrintToLog("\n", false, false);
962  PrintToLog("Simulation aborted");
963  return;
964  }
965 
966  double* temp = points->GetPoint(j);
967  itk::Point<float, 3> vertex = GetItkPoint(temp);
968  itk::Vector<double> v = GetItkVector(temp);
969 
970  itk::Vector<double, 3> dir(3);
971  if (j<numPoints-1) { dir = GetItkVector(points->GetPoint(j+1))-v; }
972  else { dir = v-GetItkVector(points->GetPoint(j-1)); }
973 
974  if ( dir.GetSquaredNorm()<0.0001 || dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2] )
975  {
976  continue;
977  }
978 
979  itk::Index<3> idx;
980  itk::ContinuousIndex<float, 3> contIndex;
981  m_TransformedMaskImage->TransformPhysicalPointToIndex(vertex, idx);
982  m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex);
983 
984  if (!m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(idx) || m_TransformedMaskImage->GetPixel(idx)<=0)
985  {
986  continue;
987  }
988 
989  // generate signal for each fiber compartment
990  for (int k=0; k<numFiberCompartments; k++)
991  {
992  m_Parameters.m_FiberModelList[k]->SetFiberDirection(dir);
993  DoubleDwiType::PixelType pix = m_CompartmentImages.at(k)->GetPixel(idx);
994  pix[g] += fiberWeight*m_SegmentVolume*m_Parameters.m_FiberModelList[k]->SimulateMeasurement(g);
995 
996  m_CompartmentImages.at(k)->SetPixel(idx, pix);
997  }
998 
999  // update fiber volume image
1000  double vol = intraAxonalVolumeImage->GetPixel(idx) + m_SegmentVolume*fiberWeight;
1001  intraAxonalVolumeImage->SetPixel(idx, vol);
1002 
1003  // we assume that the first volume is always unweighted!
1004  if (vol>maxVolume) { maxVolume = vol; }
1005  }
1006 
1007  // progress report
1008  ++disp;
1009  unsigned long newTick = 50*disp.count()/disp.expected_count();
1010  for (unsigned int tick = 0; tick<(newTick-lastTick); tick++)
1011  {
1012  PrintToLog("*", false, false, false);
1013  }
1014  lastTick = newTick;
1015  }
1016 
1017  // generate non-fiber signal
1018  ImageRegionIterator<ItkUcharImgType> it3(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion());
1019  double fact = 1; // density correction factor in mm³
1020  if (m_Parameters.m_SignalGen.m_AxonRadius<0.0001 || maxVolume>m_VoxelVolume) // the fullest voxel is always completely full
1021  fact = m_VoxelVolume/maxVolume;
1022  while(!it3.IsAtEnd())
1023  {
1024  if (it3.Get()>0)
1025  {
1026  DoubleDwiType::IndexType index = it3.GetIndex();
1027  itk::Point<double, 3> point;
1028  m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, point);
1029  if ( m_Parameters.m_SignalGen.m_DoAddMotion && g>=0
1030  && m_Parameters.m_SignalGen.m_MotionVolumes[g] )
1031  {
1032  if (m_Parameters.m_SignalGen.m_DoRandomizeMotion)
1033  {
1034  point = m_FiberBundleWorkingCopy->TransformPoint( point.GetVnlVector(), -m_Rotation[0], -m_Rotation[1], -m_Rotation[2],
1035  -m_Translation[0], -m_Translation[1], -m_Translation[2] );
1036  }
1037  else
1038  {
1039  point = m_FiberBundleWorkingCopy->TransformPoint( point.GetVnlVector(),
1040  -m_Rotation[0]*m_MotionCounter, -m_Rotation[1]*m_MotionCounter, -m_Rotation[2]*m_MotionCounter,
1041  -m_Translation[0]*m_MotionCounter, -m_Translation[1]*m_MotionCounter, -m_Translation[2]*m_MotionCounter );
1042  }
1043  }
1044 
1045  double iAxVolume = intraAxonalVolumeImage->GetPixel(index);
1046 
1047  // if volume fraction image is set use it, otherwise use scaling factor to obtain one full fiber voxel
1048  double fact2 = fact;
1049  if ( m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr
1050  && iAxVolume>0.0001 )
1051  {
1052  double val = InterpolateValue(point, m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage());
1053  if (val>=0) { fact2 = m_VoxelVolume*val/iAxVolume; }
1054  }
1055 
1056  // adjust intra-axonal image value
1057  for (int i=0; i<numFiberCompartments; i++)
1058  {
1059  DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
1060  pix[g] *= fact2;
1061  m_CompartmentImages.at(i)->SetPixel(index, pix);
1062  }
1063 
1064  // simulate other compartments
1065  SimulateExtraAxonalSignal(index, iAxVolume*fact2, g);
1066  }
1067  ++it3;
1068  }
1069  }
1070 
1071  PrintToLog("\n", false);
1072  if (this->GetAbortGenerateData())
1073  {
1074  PrintToLog("\n", false, false);
1075  PrintToLog("Simulation aborted");
1076  return;
1077  }
1078  }
1079 
1080  DoubleDwiType::Pointer doubleOutImage;
1081  double signalScale = m_Parameters.m_SignalGen.m_SignalScale;
1082  if ( m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition ) // do k-space stuff
1083  {
1084  PrintToLog("\n", false, false);
1085  PrintToLog("Simulating k-space acquisition using "
1086  +boost::lexical_cast<std::string>(m_Parameters.m_SignalGen.m_NumberOfCoils)
1087  +" coil(s)");
1088 
1089  switch (m_Parameters.m_SignalGen.m_AcquisitionType)
1090  {
1091  case SignalGenerationParameters::SingleShotEpi:
1092  {
1093  PrintToLog("Acquisition type: single shot EPI", false);
1094  break;
1095  }
1096  case SignalGenerationParameters::SpinEcho:
1097  {
1098  PrintToLog("Acquisition type: classic spin echo with cartesian k-space trajectory", false);
1099  break;
1100  }
1101  default:
1102  {
1103  PrintToLog("Acquisition type: single shot EPI", false);
1104  break;
1105  }
1106  }
1107 
1108  if (m_Parameters.m_SignalGen.m_NoiseVariance>0)
1109  PrintToLog("Simulating complex Gaussian noise", false);
1110  if (m_Parameters.m_SignalGen.m_DoSimulateRelaxation)
1111  PrintToLog("Simulating signal relaxation", false);
1112  if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull())
1113  PrintToLog("Simulating distortions", false);
1114  if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging)
1115  PrintToLog("Simulating ringing artifacts", false);
1116  if (m_Parameters.m_SignalGen.m_EddyStrength>0)
1117  PrintToLog("Simulating eddy currents", false);
1118  if (m_Parameters.m_SignalGen.m_Spikes>0)
1119  PrintToLog("Simulating spikes", false);
1120  if (m_Parameters.m_SignalGen.m_CroppingFactor<1.0)
1121  PrintToLog("Simulating aliasing artifacts", false);
1122  if (m_Parameters.m_SignalGen.m_KspaceLineOffset>0)
1123  PrintToLog("Simulating ghosts", false);
1124 
1125  doubleOutImage = SimulateKspaceAcquisition(m_CompartmentImages);
1126  signalScale = 1; // already scaled in SimulateKspaceAcquisition()
1127  }
1128  else // don't do k-space stuff, just sum compartments
1129  {
1130  PrintToLog("Summing compartments");
1131  doubleOutImage = m_CompartmentImages.at(0);
1132 
1133  for (unsigned int i=1; i<m_CompartmentImages.size(); i++)
1134  {
1135  auto adder = itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::New();
1136  adder->SetInput1(doubleOutImage);
1137  adder->SetInput2(m_CompartmentImages.at(i));
1138  adder->Update();
1139  doubleOutImage = adder->GetOutput();
1140  }
1141  }
1142  if (this->GetAbortGenerateData())
1143  {
1144  PrintToLog("\n", false, false);
1145  PrintToLog("Simulation aborted");
1146  return;
1147  }
1148 
1149  PrintToLog("Finalizing image");
1150  if (signalScale>1)
1151  PrintToLog(" Scaling signal", false);
1152  if (m_Parameters.m_NoiseModel)
1153  PrintToLog(" Adding noise", false);
1154  unsigned int window = 0;
1155  unsigned int min = itk::NumericTraits<unsigned int>::max();
1156  ImageRegionIterator<OutputImageType> it4 (m_OutputImage, m_OutputImage->GetLargestPossibleRegion());
1157  DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.m_SignalGen.GetNumVolumes());
1158  boost::progress_display disp2(m_OutputImage->GetLargestPossibleRegion().GetNumberOfPixels());
1159 
1160  PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false);
1161  PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false);
1162  int lastTick = 0;
1163 
1164  while(!it4.IsAtEnd())
1165  {
1166  if (this->GetAbortGenerateData())
1167  {
1168  PrintToLog("\n", false, false);
1169  PrintToLog("Simulation aborted");
1170  return;
1171  }
1172 
1173  ++disp2;
1174  unsigned long newTick = 50*disp2.count()/disp2.expected_count();
1175  for (unsigned long tick = 0; tick<(newTick-lastTick); tick++)
1176  PrintToLog("*", false, false, false);
1177  lastTick = newTick;
1178 
1179  typename OutputImageType::IndexType index = it4.GetIndex();
1180  signal = doubleOutImage->GetPixel(index)*signalScale;
1181 
1182  if (m_Parameters.m_NoiseModel)
1183  m_Parameters.m_NoiseModel->AddNoise(signal);
1184 
1185  for (unsigned int i=0; i<signal.Size(); i++)
1186  {
1187  if (signal[i]>0)
1188  signal[i] = floor(signal[i]+0.5);
1189  else
1190  signal[i] = ceil(signal[i]-0.5);
1191 
1192  if ( (!m_Parameters.m_SignalGen.IsBaselineIndex(i) || signal.Size()==1) && signal[i]>window)
1193  window = signal[i];
1194  if ( (!m_Parameters.m_SignalGen.IsBaselineIndex(i) || signal.Size()==1) && signal[i]<min)
1195  min = signal[i];
1196  }
1197  it4.Set(signal);
1198  ++it4;
1199  }
1200  window -= min;
1201  unsigned int level = window/2 + min;
1202  m_LevelWindow.SetLevelWindow(level, window);
1203  this->SetNthOutput(0, m_OutputImage);
1204 
1205  PrintToLog("\n", false);
1206  PrintToLog("Finished simulation");
1207  m_TimeProbe.Stop();
1208 
1209  if (m_Parameters.m_SignalGen.m_DoAddMotion)
1210  {
1211  PrintToLog("\nHead motion log:", false);
1212  PrintToLog(m_MotionLog, false, false);
1213  }
1214 
1215  if (m_Parameters.m_SignalGen.m_Spikes>0)
1216  {
1217  PrintToLog("\nSpike log:", false);
1218  PrintToLog(m_SpikeLog, false, false);
1219  }
1220 
1221  if (m_Logfile.is_open())
1222  m_Logfile.close();
1223  } // Heilliger Spaghetti-Code, Batman! todo/mdh
1224 
1225 
1226  template< class PixelType >
1227  void TractsToDWIImageFilter< PixelType >::PrintToLog(string m, bool addTime, bool linebreak, bool stdOut)
1228  {
1229  // timestamp
1230  if (addTime)
1231  {
1232  m_Logfile << this->GetTime() << " > ";
1233  m_StatusText += this->GetTime() + " > ";
1234  if (stdOut)
1235  std::cout << this->GetTime() << " > ";
1236  }
1237 
1238  // message
1239  if (m_Logfile.is_open())
1240  m_Logfile << m;
1241  m_StatusText += m;
1242  if (stdOut)
1243  std::cout << m;
1244 
1245  // new line
1246  if (linebreak)
1247  {
1248  if (m_Logfile.is_open())
1249  m_Logfile << "\n";
1250  m_StatusText += "\n";
1251  if (stdOut)
1252  std::cout << "\n";
1253  }
1254  }
1255 
1256  template< class PixelType >
1257  void TractsToDWIImageFilter< PixelType >::SimulateMotion(int g)
1258  {
1259  // is motion artifact enabled?
1260  // is the current volume g affected by motion?
1261  if ( m_Parameters.m_SignalGen.m_DoAddMotion
1262  && m_Parameters.m_SignalGen.m_MotionVolumes[g]
1263  && g<m_Parameters.m_SignalGen.GetNumVolumes() )
1264  {
1265  if ( m_Parameters.m_SignalGen.m_DoRandomizeMotion )
1266  {
1267  // either undo last transform or work on fresh copy of untransformed fibers
1268  m_FiberBundleTransformed = m_FiberBundleWorkingCopy->GetDeepCopy();
1269 
1270  m_Rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[0]*2)
1271  -m_Parameters.m_SignalGen.m_Rotation[0];
1272  m_Rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[1]*2)
1273  -m_Parameters.m_SignalGen.m_Rotation[1];
1274  m_Rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[2]*2)
1275  -m_Parameters.m_SignalGen.m_Rotation[2];
1276 
1277  m_Translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[0]*2)
1278  -m_Parameters.m_SignalGen.m_Translation[0];
1279  m_Translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[1]*2)
1280  -m_Parameters.m_SignalGen.m_Translation[1];
1281  m_Translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[2]*2)
1282  -m_Parameters.m_SignalGen.m_Translation[2];
1283  }
1284  else
1285  {
1286  m_Rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes;
1287  m_Translation = m_Parameters.m_SignalGen.m_Translation / m_NumMotionVolumes;
1288  m_MotionCounter++;
1289  }
1290 
1291  // move mask image
1292  if (m_MaskImageSet)
1293  {
1294  ImageRegionIterator<ItkUcharImgType> maskIt(m_UpsampledMaskImage, m_UpsampledMaskImage->GetLargestPossibleRegion());
1295  m_TransformedMaskImage->FillBuffer(0);
1296 
1297  while(!maskIt.IsAtEnd())
1298  {
1299  if (maskIt.Get()<=0)
1300  {
1301  ++maskIt;
1302  continue;
1303  }
1304 
1305  DoubleDwiType::IndexType index = maskIt.GetIndex();
1306  itk::Point<double, 3> point;
1307  m_UpsampledMaskImage->TransformIndexToPhysicalPoint(index, point);
1308  if (m_Parameters.m_SignalGen.m_DoRandomizeMotion)
1309  {
1310  point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(),
1311  m_Rotation[0],m_Rotation[1],m_Rotation[2],
1312  m_Translation[0],m_Translation[1],m_Translation[2]);
1313  }
1314  else
1315  {
1316  point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(),
1317  m_Rotation[0]*m_MotionCounter,m_Rotation[1]*m_MotionCounter,m_Rotation[2]*m_MotionCounter,
1318  m_Translation[0]*m_MotionCounter,m_Translation[1]*m_MotionCounter,m_Translation[2]*m_MotionCounter);
1319  }
1320 
1321  m_TransformedMaskImage->TransformPhysicalPointToIndex(point, index);
1322 
1323  if (m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(index))
1324  { m_TransformedMaskImage->SetPixel(index,100); }
1325  ++maskIt;
1326  }
1327  }
1328 
1329  if (m_Parameters.m_SignalGen.m_DoRandomizeMotion)
1330  {
1331  m_Rotations.push_back(m_Rotation);
1332  m_Translations.push_back(m_Translation);
1333 
1334  m_MotionLog += boost::lexical_cast<std::string>(g) + " rotation: " + boost::lexical_cast<std::string>(m_Rotation[0])
1335  + "," + boost::lexical_cast<std::string>(m_Rotation[1])
1336  + "," + boost::lexical_cast<std::string>(m_Rotation[2]) + ";";
1337 
1338  m_MotionLog += " translation: " + boost::lexical_cast<std::string>(m_Translation[0])
1339  + "," + boost::lexical_cast<std::string>(m_Translation[1])
1340  + "," + boost::lexical_cast<std::string>(m_Translation[2]) + "\n";
1341  }
1342  else
1343  {
1344  m_Rotations.push_back(m_Rotation*m_MotionCounter);
1345  m_Translations.push_back(m_Translation*m_MotionCounter);
1346 
1347  m_MotionLog += boost::lexical_cast<std::string>(g) + " rotation: " + boost::lexical_cast<std::string>(m_Rotation[0]*m_MotionCounter)
1348  + "," + boost::lexical_cast<std::string>(m_Rotation[1]*m_MotionCounter)
1349  + "," + boost::lexical_cast<std::string>(m_Rotation[2]*m_MotionCounter) + ";";
1350 
1351  m_MotionLog += " translation: " + boost::lexical_cast<std::string>(m_Translation[0]*m_MotionCounter)
1352  + "," + boost::lexical_cast<std::string>(m_Translation[1]*m_MotionCounter)
1353  + "," + boost::lexical_cast<std::string>(m_Translation[2]*m_MotionCounter) + "\n";
1354  }
1355 
1356  m_FiberBundleTransformed->TransformFibers(m_Rotation[0],m_Rotation[1],m_Rotation[2],m_Translation[0],m_Translation[1],m_Translation[2]);
1357  }
1358  else
1359  {
1360  m_Rotation.Fill(0.0);
1361  m_Translation.Fill(0.0);
1362  m_Rotations.push_back(m_Rotation);
1363  m_Translations.push_back(m_Translation);
1364 
1365  m_MotionLog += boost::lexical_cast<std::string>(g) + " rotation: " + boost::lexical_cast<std::string>(m_Rotation[0])
1366  + "," + boost::lexical_cast<std::string>(m_Rotation[1])
1367  + "," + boost::lexical_cast<std::string>(m_Rotation[2]) + ";";
1368 
1369  m_MotionLog += " translation: " + boost::lexical_cast<std::string>(m_Translation[0])
1370  + "," + boost::lexical_cast<std::string>(m_Translation[1])
1371  + "," + boost::lexical_cast<std::string>(m_Translation[2]) + "\n";
1372  }
1373  }
1374 
1375  template< class PixelType >
1376  void TractsToDWIImageFilter< PixelType >::
1377  SimulateExtraAxonalSignal(ItkUcharImgType::IndexType index, double intraAxonalVolume, int g)
1378  {
1379  int numFiberCompartments = m_Parameters.m_FiberModelList.size();
1380  int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size();
1381 
1382  if (intraAxonalVolume>0.0001 && m_Parameters.m_SignalGen.m_DoDisablePartialVolume) // only fiber in voxel
1383  {
1384  DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index);
1385  if (g>=0)
1386  pix[g] *= m_VoxelVolume/intraAxonalVolume;
1387  else
1388  pix *= m_VoxelVolume/intraAxonalVolume;
1389  m_CompartmentImages.at(0)->SetPixel(index, pix);
1390  if (g==0)
1391  m_VolumeFractions.at(0)->SetPixel(index, 1);
1392  for (int i=1; i<numFiberCompartments; i++)
1393  {
1394  DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
1395  if (g>=0)
1396  pix[g] = 0.0;
1397  else
1398  pix.Fill(0.0);
1399  m_CompartmentImages.at(i)->SetPixel(index, pix);
1400  }
1401  }
1402  else
1403  {
1404  if (g==0)
1405  {
1406  m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/m_VoxelVolume);
1407  }
1408 
1409  // get non-transformed point (remove headmotion tranformation)
1410  // this point can then be transformed to each of the original images, regardless of their geometry
1411  itk::Point<double, 3> point;
1412  m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, point);
1413 
1414  if ( m_Parameters.m_SignalGen.m_DoAddMotion && g>=0
1415  && m_Parameters.m_SignalGen.m_MotionVolumes[g] )
1416  {
1417  if (m_Parameters.m_SignalGen.m_DoRandomizeMotion)
1418  {
1419  point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(),
1420  -m_Rotation[0],-m_Rotation[1],-m_Rotation[2],
1421  -m_Translation[0],-m_Translation[1],-m_Translation[2]);
1422  }
1423  else
1424  {
1425  point = m_FiberBundleWorkingCopy->TransformPoint(point.GetVnlVector(),
1426  -m_Rotation[0]*m_MotionCounter,-m_Rotation[1]*m_MotionCounter,-m_Rotation[2]*m_MotionCounter,
1427  -m_Translation[0]*m_MotionCounter,-m_Translation[1]*m_MotionCounter,-m_Translation[2]*m_MotionCounter);
1428  }
1429  }
1430 
1431  if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume)
1432  {
1433  int maxVolumeIndex = 0;
1434  double maxWeight = 0;
1435  for (int i=0; i<numNonFiberCompartments; i++)
1436  {
1437  double weight = 0;
1438  if (numNonFiberCompartments>1)
1439  {
1440  double val = InterpolateValue(point, m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage());
1441  if (val<0)
1442  continue;
1443  else
1444  weight = val;
1445  }
1446 
1447  if (weight>maxWeight)
1448  {
1449  maxWeight = weight;
1450  maxVolumeIndex = i;
1451  }
1452  }
1453 
1454  DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(maxVolumeIndex+numFiberCompartments);
1455  DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index);
1456 
1457  if (g>=0)
1458  pix[g] += m_Parameters.m_NonFiberModelList[maxVolumeIndex]->SimulateMeasurement(g)*m_VoxelVolume;
1459  else
1460  pix += m_Parameters.m_NonFiberModelList[maxVolumeIndex]->SimulateMeasurement()*m_VoxelVolume;
1461  doubleDwi->SetPixel(index, pix);
1462  if (g==0)
1463  m_VolumeFractions.at(maxVolumeIndex+numFiberCompartments)->SetPixel(index, 1);
1464  }
1465  else
1466  {
1467  double extraAxonalVolume = m_VoxelVolume-intraAxonalVolume; // non-fiber volume
1468  if (extraAxonalVolume<0)
1469  {
1470  MITK_ERROR << "Corrupted intra-axonal signal voxel detected. Fiber volume larger voxel volume!";
1471  extraAxonalVolume = 0;
1472  }
1473  double interAxonalVolume = 0;
1474  if (numFiberCompartments>1)
1475  interAxonalVolume = extraAxonalVolume * intraAxonalVolume/m_VoxelVolume; // inter-axonal fraction of non fiber compartment
1476  double other = extraAxonalVolume - interAxonalVolume; // rest of compartment
1477  if (other<0)
1478  {
1479  MITK_ERROR << "Corrupted signal voxel detected. Fiber volume larger voxel volume!";
1480  other = 0;
1481  }
1482 
1483  // adjust non-fiber and intra-axonal signal
1484  for (int i=1; i<numFiberCompartments; i++)
1485  {
1486  double weight = interAxonalVolume;
1487  DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index);
1488  if (intraAxonalVolume>0) // remove scaling by intra-axonal volume from inter-axonal compartment
1489  {
1490  if (g>=0)
1491  pix[g] /= intraAxonalVolume;
1492  else
1493  pix /= intraAxonalVolume;
1494  }
1495 
1496  if (m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()!=nullptr)
1497  {
1498  double val = InterpolateValue(point, m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage());
1499  if (val<0)
1500  continue;
1501  else
1502  weight = val*m_VoxelVolume;
1503  }
1504 
1505  if (g>=0)
1506  pix[g] *= weight;
1507  else
1508  pix *= weight;
1509  m_CompartmentImages.at(i)->SetPixel(index, pix);
1510  if (g==0)
1511  m_VolumeFractions.at(i)->SetPixel(index, weight/m_VoxelVolume);
1512  }
1513 
1514  for (int i=0; i<numNonFiberCompartments; i++)
1515  {
1516  double weight = other;
1517  DoubleDwiType::PixelType pix = m_CompartmentImages.at(i+numFiberCompartments)->GetPixel(index);
1518  if (m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()!=nullptr)
1519  {
1520  double val = InterpolateValue(point, m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage());
1521  if (val<0)
1522  continue;
1523  else
1524  weight = val*m_VoxelVolume;
1525 
1526  if (m_UseRelativeNonFiberVolumeFractions)
1527  weight *= other/m_VoxelVolume;
1528  }
1529 
1530  if (g>=0)
1531  pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g)*weight;
1532  else
1533  pix += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement()*weight;
1534  m_CompartmentImages.at(i+numFiberCompartments)->SetPixel(index, pix);
1535  if (g==0)
1536  m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, weight/m_VoxelVolume);
1537  }
1538  }
1539  }
1540  }
1541 
1542  template< class PixelType >
1543  double TractsToDWIImageFilter< PixelType >::
1544  InterpolateValue(itk::Point<float, 3> itkP, ItkDoubleImgType::Pointer img)
1545  {
1546  itk::Index<3> idx;
1547  itk::ContinuousIndex< double, 3> cIdx;
1548  img->TransformPhysicalPointToIndex(itkP, idx);
1549  img->TransformPhysicalPointToContinuousIndex(itkP, cIdx);
1550 
1551  double pix = -1;
1552  if ( img->GetLargestPossibleRegion().IsInside(idx) )
1553  pix = img->GetPixel(idx);
1554  else
1555  return pix;
1556 
1557  double frac_x = cIdx[0] - idx[0];
1558  double frac_y = cIdx[1] - idx[1];
1559  double frac_z = cIdx[2] - idx[2];
1560  if (frac_x<0)
1561  {
1562  idx[0] -= 1;
1563  frac_x += 1;
1564  }
1565  if (frac_y<0)
1566  {
1567  idx[1] -= 1;
1568  frac_y += 1;
1569  }
1570  if (frac_z<0)
1571  {
1572  idx[2] -= 1;
1573  frac_z += 1;
1574  }
1575  frac_x = 1-frac_x;
1576  frac_y = 1-frac_y;
1577  frac_z = 1-frac_z;
1578 
1579  // int coordinates inside image?
1580  if (idx[0] >= 0 && idx[0] < img->GetLargestPossibleRegion().GetSize(0)-1 &&
1581  idx[1] >= 0 && idx[1] < img->GetLargestPossibleRegion().GetSize(1)-1 &&
1582  idx[2] >= 0 && idx[2] < img->GetLargestPossibleRegion().GetSize(2)-1)
1583  {
1584  vnl_vector_fixed<double, 8> interpWeights;
1585  interpWeights[0] = ( frac_x)*( frac_y)*( frac_z);
1586  interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z);
1587  interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z);
1588  interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z);
1589  interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z);
1590  interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z);
1591  interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z);
1592  interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z);
1593 
1594  pix = img->GetPixel(idx) * interpWeights[0];
1595  ItkDoubleImgType::IndexType tmpIdx = idx; tmpIdx[0]++;
1596  pix += img->GetPixel(tmpIdx) * interpWeights[1];
1597  tmpIdx = idx; tmpIdx[1]++;
1598  pix += img->GetPixel(tmpIdx) * interpWeights[2];
1599  tmpIdx = idx; tmpIdx[2]++;
1600  pix += img->GetPixel(tmpIdx) * interpWeights[3];
1601  tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++;
1602  pix += img->GetPixel(tmpIdx) * interpWeights[4];
1603  tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++;
1604  pix += img->GetPixel(tmpIdx) * interpWeights[5];
1605  tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++;
1606  pix += img->GetPixel(tmpIdx) * interpWeights[6];
1607  tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++;
1608  pix += img->GetPixel(tmpIdx) * interpWeights[7];
1609  }
1610 
1611  return pix;
1612  }
1613 
1614  template< class PixelType >
1615  itk::Point<float, 3> TractsToDWIImageFilter< PixelType >::GetItkPoint(double point[3])
1616  {
1617  itk::Point<float, 3> itkPoint;
1618  itkPoint[0] = point[0];
1619  itkPoint[1] = point[1];
1620  itkPoint[2] = point[2];
1621  return itkPoint;
1622  }
1623 
1624  template< class PixelType >
1625  itk::Vector<double, 3> TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3])
1626  {
1627  itk::Vector<double, 3> itkVector;
1628  itkVector[0] = point[0];
1629  itkVector[1] = point[1];
1630  itkVector[2] = point[2];
1631  return itkVector;
1632  }
1633 
1634  template< class PixelType >
1635  vnl_vector_fixed<double, 3> TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3])
1636  {
1637  vnl_vector_fixed<double, 3> vnlVector;
1638  vnlVector[0] = point[0];
1639  vnlVector[1] = point[1];
1640  vnlVector[2] = point[2];
1641  return vnlVector;
1642  }
1643 
1644  template< class PixelType >
1645  vnl_vector_fixed<double, 3> TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector<float,3>& vector)
1646  {
1647  vnl_vector_fixed<double, 3> vnlVector;
1648  vnlVector[0] = vector[0];
1649  vnlVector[1] = vector[1];
1650  vnlVector[2] = vector[2];
1651  return vnlVector;
1652  }
1653 
1654  template< class PixelType >
1655  double TractsToDWIImageFilter< PixelType >::RoundToNearest(double num) {
1656  return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5);
1657  }
1658 
1659  template< class PixelType >
1660  std::string TractsToDWIImageFilter< PixelType >::GetTime()
1661  {
1662  m_TimeProbe.Stop();
1663  unsigned long total = RoundToNearest(m_TimeProbe.GetTotal());
1664  unsigned long hours = total/3600;
1665  unsigned long minutes = (total%3600)/60;
1666  unsigned long seconds = total%60;
1667  std::string out = "";
1668  out.append(boost::lexical_cast<std::string>(hours));
1669  out.append(":");
1670  out.append(boost::lexical_cast<std::string>(minutes));
1671  out.append(":");
1672  out.append(boost::lexical_cast<std::string>(seconds));
1673  m_TimeProbe.Start();
1674  return out;
1675  }
1676 
1677 }
itk::TractsToDWIImageFilter::PrintToLog
void PrintToLog(string m, bool addTime=true, bool linebreak=true, bool stdOut=true)
Definition: itkTractsToDWIImageFilter.cpp:1227
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
itk::TractsToDWIImageFilter::SimulateExtraAxonalSignal
void SimulateExtraAxonalSignal(ItkUcharImgType::IndexType index, double intraAxonalVolume, int g=-1)
Definition: itkTractsToDWIImageFilter.cpp:1377
itk::TractsToDWIImageFilter
Generates artificial diffusion weighted image volume from the input fiberbundle using a generic multi...
Definition: itkTractsToDWIImageFilter.h:38
mitk::IOUtil::GetTempPath
static std::string GetTempPath()
Definition: mitkIOUtil.cpp:372
itk::TractsToDWIImageFilter::NormalizeInsideMask
ItkDoubleImgType::Pointer NormalizeInsideMask(ItkDoubleImgType::Pointer image)
Definition: itkTractsToDWIImageFilter.cpp:343
MITK_ERROR
#define MITK_ERROR
Definition: mitkLogMacros.h:24
mitk::PointSet::New
static Pointer New()
itk::TractsToDWIImageFilter::SimulateKspaceAcquisition
DoubleDwiType::Pointer SimulateKspaceAcquisition(std::vector< DoubleDwiType::Pointer > &images)
Definition: itkTractsToDWIImageFilter.cpp:78
itk::TractsToDWIImageFilter::m_RandGen
itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen
Definition: itkTractsToDWIImageFilter.h:155
itk::DftImageFilter::New
static Pointer New()
itk::KspaceImageFilter::New
static Pointer New()
mitk::DiffusionSignalModel
Abstract class for diffusion signal models.
Definition: mitkDiffusionSignalModel.h:35
itk::TractsToDWIImageFilter::DoubleVectorType
itk::Vector< double, 3 > DoubleVectorType
Definition: itkTractsToDWIImageFilter.h:59
itk::TractsToDWIImageFilter::GetItkPoint
itk::Point< float, 3 > GetItkPoint(double point[3])
Definition: itkTractsToDWIImageFilter.cpp:1615
itk::TractsToDWIImageFilter::InterpolateValue
double InterpolateValue(itk::Point< float, 3 > itkP, ItkDoubleImgType::Pointer img)
Definition: itkTractsToDWIImageFilter.cpp:1544
max
static T max(T x, T y)
Definition: svm.cpp:70
min
static T min(T x, T y)
Definition: svm.cpp:67
itk::ResampleDwiImageFilter
Resample DWI channel by channel.
Definition: itkResampleDwiImageFilter.h:46
PixelType
unsigned short PixelType
Definition: TumorProgressionMapping.cpp:55
itk::TractsToDWIImageFilter::GetVnlVector
vnl_vector_fixed< double, 3 > GetVnlVector(double point[3])
Definition: itkTractsToDWIImageFilter.cpp:1635
itk::TractsToDWIImageFilter::GetItkVector
itk::Vector< double, 3 > GetItkVector(double point[3])
Definition: itkTractsToDWIImageFilter.cpp:1625
images
section MAP_FRAME_Mapper_Settings Mapper settings For the mapping of corrected images
Definition: org.mitk.gui.qt.matchpoint.framereg/documentation/UserManual/Manual.dox:48
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.