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