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
itkDwiPhantomGenerationFilter.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 #ifndef __itkDwiPhantomGenerationFilter_cpp
17 #define __itkDwiPhantomGenerationFilter_cpp
18 
19 #include <time.h>
20 #include <stdio.h>
21 #include <stdlib.h>
22 
23 #include "itkDwiPhantomGenerationFilter.h"
24 #include <itkImageRegionConstIterator.h>
25 #include <itkImageRegionConstIteratorWithIndex.h>
26 #include <itkImageRegionIterator.h>
27 #include <itkOrientationDistributionFunction.h>
28 
29 #include <vtkSmartPointer.h>
30 #include <vtkPolyData.h>
31 #include <vtkCellArray.h>
32 #include <vtkPoints.h>
33 #include <vtkPolyLine.h>
34 
35 #include <boost/progress.hpp>
36 
37 namespace itk {
38 
39 
40  template< class TOutputScalarType >
41  DwiPhantomGenerationFilter< TOutputScalarType >
42  ::DwiPhantomGenerationFilter()
43  : m_BValue(1000)
44  , m_SignalScale(1000)
45  , m_BaselineImages(0)
46  , m_MaxBaseline(0)
47  , m_MeanBaseline(0)
48  , m_NoiseVariance(0.004)
49  , m_AddNoiseFlag(true) // parameters.m_Misc.m_CheckAddNoiseBox
50  , m_GreyMatterAdc(0.01)
51  , m_SimulateBaseline(true)
52  , m_DefaultBaseline(1000)
53  {
54  this->SetNumberOfRequiredOutputs (1);
55  m_Spacing.Fill(2.5); m_Origin.Fill(0.0);
56  m_DirectionMatrix.SetIdentity();
57  m_ImageRegion.SetSize(0, 10);
58  m_ImageRegion.SetSize(1, 10);
59  m_ImageRegion.SetSize(2, 10);
60 
61  typename OutputImageType::Pointer outImage = OutputImageType::New();
62 
63  outImage->SetSpacing( m_Spacing ); // Set the image spacing
64  outImage->SetOrigin( m_Origin ); // Set the image origin
65  outImage->SetDirection( m_DirectionMatrix ); // Set the image direction
66  outImage->SetLargestPossibleRegion( m_ImageRegion );
67  outImage->SetBufferedRegion( m_ImageRegion );
68  outImage->SetRequestedRegion( m_ImageRegion );
69  outImage->SetVectorLength(QBALL_ODFSIZE);
70  outImage->Allocate();
71  // ITKv4 migration fix : removing OutputImageType::PixelType(0.0)
72  // the conversion is handled internally by the itk::Image
73  typename OutputImageType::PixelType fillValue(0.0);
74  outImage->FillBuffer( fillValue );
75 
76  this->SetNthOutput (0, outImage);
77 
78 
79  }
80 
81  template< class TOutputScalarType >
82  void DwiPhantomGenerationFilter< TOutputScalarType >
83  ::GenerateTensors()
84  {
85  MITK_INFO << "Generating tensors";
86 
87  for (int i=0; i<m_SignalRegions.size(); i++)
88  {
89  float ADC = m_TensorADC.at(i);
90  float FA = m_TensorFA.at(i);
91  itk::DiffusionTensor3D<float> kernel;
92  kernel.Fill(0);
93  float e1=ADC*(1+2*FA/sqrt(3-2*FA*FA));
94  float e2=ADC*(1-FA/sqrt(3-2*FA*FA));
95  float e3=e2;
96  kernel.SetElement(0,e1);
97  kernel.SetElement(3,e2);
98  kernel.SetElement(5,e3);
99 
100  if (m_SimulateBaseline)
101  {
102  double l2 = GetTensorL2Norm(kernel);
103  if (l2>m_MaxBaseline)
104  m_MaxBaseline = l2;
105  }
106 
107  MITK_INFO << "Kernel FA: " << kernel.GetFractionalAnisotropy();
108  vnl_vector_fixed<double, 3> kernelDir; kernelDir[0]=1; kernelDir[1]=0; kernelDir[2]=0;
109 
110  itk::DiffusionTensor3D<float> tensor;
111  vnl_vector_fixed<double, 3> dir = m_TensorDirection.at(i);
112  MITK_INFO << "Tensor direction: " << dir;
113 
114  dir.normalize();
115 
116  vnl_vector_fixed<double, 3> axis = vnl_cross_3d(kernelDir, dir); axis.normalize();
117  vnl_quaternion<double> rotation(axis, acos(dot_product(kernelDir, dir)));
118  rotation.normalize();
119  vnl_matrix_fixed<double, 3, 3> matrix = rotation.rotation_matrix_transpose();
120 
121  vnl_matrix_fixed<double, 3, 3> tensorMatrix;
122  tensorMatrix[0][0] = kernel[0]; tensorMatrix[0][1] = kernel[1]; tensorMatrix[0][2] = kernel[2];
123  tensorMatrix[1][0] = kernel[1]; tensorMatrix[1][1] = kernel[3]; tensorMatrix[1][2] = kernel[4];
124  tensorMatrix[2][0] = kernel[2]; tensorMatrix[2][1] = kernel[4]; tensorMatrix[2][2] = kernel[5];
125 
126  tensorMatrix = matrix.transpose()*tensorMatrix*matrix;
127  tensor[0] = tensorMatrix[0][0]; tensor[1] = tensorMatrix[0][1]; tensor[2] = tensorMatrix[0][2];
128  tensor[3] = tensorMatrix[1][1]; tensor[4] = tensorMatrix[1][2]; tensor[5] = tensorMatrix[2][2];
129 
130  m_TensorList.push_back(tensor);
131  }
132  }
133 
134  template< class TOutputScalarType >
135  void DwiPhantomGenerationFilter< TOutputScalarType >::AddNoise(typename OutputImageType::PixelType& pix)
136  {
137  for( unsigned int i=0; i<m_GradientList.size(); i++)
138  {
139  float signal = pix[i];
140  float val = sqrt(pow(signal + m_SignalScale*m_RandGen->GetNormalVariate(0.0, m_NoiseVariance), 2) + pow(m_SignalScale*m_RandGen->GetNormalVariate(0.0, m_NoiseVariance),2));
141  pix[i] += val;
142  }
143  }
144 
145  template< class TOutputScalarType >
146  typename DwiPhantomGenerationFilter< TOutputScalarType >::OutputImageType::PixelType
147  DwiPhantomGenerationFilter< TOutputScalarType >::SimulateMeasurement(itk::DiffusionTensor3D<float>& T, float weight)
148  {
149  typename OutputImageType::PixelType out;
150  out.SetSize(m_GradientList.size());
151  out.Fill(0);
152 
153  TOutputScalarType s0 = m_DefaultBaseline;
154  if (m_SimulateBaseline)
155  s0 = (GetTensorL2Norm(T)/m_MaxBaseline)*m_SignalScale;
156 
157  for( unsigned int i=0; i<m_GradientList.size(); i++)
158  {
159  GradientType g = m_GradientList[i];
160 
161  if (g.GetNorm()>0.0001)
162  {
163  itk::DiffusionTensor3D<float> S;
164  S[0] = g[0]*g[0];
165  S[1] = g[1]*g[0];
166  S[2] = g[2]*g[0];
167  S[3] = g[1]*g[1];
168  S[4] = g[2]*g[1];
169  S[5] = g[2]*g[2];
170 
171  double D = T[0]*S[0] + T[1]*S[1] + T[2]*S[2] +
172  T[1]*S[1] + T[3]*S[3] + T[4]*S[4] +
173  T[2]*S[2] + T[4]*S[4] + T[5]*S[5];
174 
175  // check for corrupted tensor and generate signal
176  if (D>=0)
177  {
178  D = weight*s0*exp ( -m_BValue * D );
179  out[i] = static_cast<TOutputScalarType>( D );
180  }
181  }
182  else
183  out[i] = s0;
184  }
185 
186  return out;
187  }
188 
189  template< class TOutputScalarType >
190  void DwiPhantomGenerationFilter< TOutputScalarType >
191  ::GenerateData()
192  {
193  if (m_NoiseVariance < 0 && m_AddNoiseFlag)
194  {
195  m_NoiseVariance = 0.001;
196  }
197 
198  if (!m_SimulateBaseline)
199  {
200  MITK_INFO << "Baseline image values are set to default. Noise variance value is treated as SNR!";
201  if (m_NoiseVariance <= 0 && m_AddNoiseFlag)
202  {
203  m_NoiseVariance = 0.0001;
204  }
205  else if (m_NoiseVariance>99 && m_AddNoiseFlag)
206  {
207  m_NoiseVariance = 0;
208  }
209  else
210  {
211  m_NoiseVariance = m_DefaultBaseline/(m_NoiseVariance*m_SignalScale);
212  m_NoiseVariance *= m_NoiseVariance;
213  }
214  }
215 
216  m_RandGen = Statistics::MersenneTwisterRandomVariateGenerator::New();
217  m_RandGen->SetSeed();
218 
219  typename OutputImageType::Pointer outImage = OutputImageType::New();
220  outImage->SetSpacing( m_Spacing );
221  outImage->SetOrigin( m_Origin );
222  outImage->SetDirection( m_DirectionMatrix );
223  outImage->SetLargestPossibleRegion( m_ImageRegion );
224  outImage->SetBufferedRegion( m_ImageRegion );
225  outImage->SetRequestedRegion( m_ImageRegion );
226  outImage->SetVectorLength(m_GradientList.size());
227  outImage->Allocate();
228  typename OutputImageType::PixelType pix;
229  pix.SetSize(m_GradientList.size());
230  pix.Fill(0.0);
231  outImage->FillBuffer(pix);
232  this->SetNthOutput (0, outImage);
233 
234  double minSpacing = m_Spacing[0];
235  if (m_Spacing[1]<minSpacing)
236  minSpacing = m_Spacing[1];
237  if (m_Spacing[2]<minSpacing)
238  minSpacing = m_Spacing[2];
239 
240  m_DirectionImageContainer = ItkDirectionImageContainer::New();
241  for (int i=0; i<m_SignalRegions.size(); i++)
242  {
243  itk::Vector< float, 3 > nullVec; nullVec.Fill(0.0);
244  ItkDirectionImage::Pointer img = ItkDirectionImage::New();
245  img->SetSpacing( m_Spacing );
246  img->SetOrigin( m_Origin );
247  img->SetDirection( m_DirectionMatrix );
248  img->SetRegions( m_ImageRegion );
249  img->Allocate();
250  img->FillBuffer(nullVec);
251  m_DirectionImageContainer->InsertElement(m_DirectionImageContainer->Size(), img);
252  }
253  m_NumDirectionsImage = ItkUcharImgType::New();
254  m_NumDirectionsImage->SetSpacing( m_Spacing );
255  m_NumDirectionsImage->SetOrigin( m_Origin );
256  m_NumDirectionsImage->SetDirection( m_DirectionMatrix );
257  m_NumDirectionsImage->SetRegions( m_ImageRegion );
258  m_NumDirectionsImage->Allocate();
259  m_NumDirectionsImage->FillBuffer(0);
260 
261  m_SNRImage = ItkFloatImgType::New();
262  m_SNRImage->SetSpacing( m_Spacing );
263  m_SNRImage->SetOrigin( m_Origin );
264  m_SNRImage->SetDirection( m_DirectionMatrix );
265  m_SNRImage->SetRegions( m_ImageRegion );
266  m_SNRImage->Allocate();
267  m_SNRImage->FillBuffer(0);
268 
269  vtkSmartPointer<vtkCellArray> m_VtkCellArray = vtkSmartPointer<vtkCellArray>::New();
270  vtkSmartPointer<vtkPoints> m_VtkPoints = vtkSmartPointer<vtkPoints>::New();
271 
272  m_BaselineImages = 0;
273  for( unsigned int i=0; i<m_GradientList.size(); i++)
274  if (m_GradientList[i].GetNorm()<=0.0001)
275  m_BaselineImages++;
276 
277  typedef ImageRegionIterator<OutputImageType> IteratorOutputType;
278  IteratorOutputType it (outImage, m_ImageRegion);
279 
280  // isotropic tensor
281  itk::DiffusionTensor3D<float> isoTensor;
282  isoTensor.Fill(0);
283  float e1 = m_GreyMatterAdc;
284  float e2 = m_GreyMatterAdc;
285  float e3 = m_GreyMatterAdc;
286  isoTensor.SetElement(0,e1);
287  isoTensor.SetElement(3,e2);
288  isoTensor.SetElement(5,e3);
289  m_MaxBaseline = GetTensorL2Norm(isoTensor);
290 
291  GenerateTensors();
292 
293  // simulate measurement
294  m_MeanBaseline = 0;
295  double noiseStdev = sqrt(m_NoiseVariance);
296  while(!it.IsAtEnd())
297  {
298  pix = it.Get();
299  typename OutputImageType::IndexType index = it.GetIndex();
300 
301  int numDirs = 0;
302  for (int i=0; i<m_SignalRegions.size(); i++)
303  {
304  ItkUcharImgType::Pointer region = m_SignalRegions.at(i);
305 
306  if (region->GetPixel(index)!=0)
307  {
308  numDirs++;
309  pix += SimulateMeasurement(m_TensorList[i], m_TensorWeight[i]);
310 
311  // set direction image pixel
312  ItkDirectionImage::Pointer img = m_DirectionImageContainer->GetElement(i);
313  itk::Vector< float, 3 > pixel = img->GetPixel(index);
314  vnl_vector_fixed<double, 3> dir = m_TensorDirection.at(i);
315  dir.normalize();
316  dir *= m_TensorWeight.at(i);
317  pixel.SetElement(0, dir[0]);
318  pixel.SetElement(1, dir[1]);
319  pixel.SetElement(2, dir[2]);
320  img->SetPixel(index, pixel);
321 
322  vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
323  itk::ContinuousIndex<double, 3> center;
324  center[0] = index[0];
325  center[1] = index[1];
326  center[2] = index[2];
327  itk::Point<double> worldCenter;
328  outImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter );
329  itk::Point<double> worldStart;
330  worldStart[0] = worldCenter[0]-dir[0]/2 * minSpacing;
331  worldStart[1] = worldCenter[1]-dir[1]/2 * minSpacing;
332  worldStart[2] = worldCenter[2]-dir[2]/2 * minSpacing;
333  vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer());
334  container->GetPointIds()->InsertNextId(id);
335  itk::Point<double> worldEnd;
336  worldEnd[0] = worldCenter[0]+dir[0]/2 * minSpacing;
337  worldEnd[1] = worldCenter[1]+dir[1]/2 * minSpacing;
338  worldEnd[2] = worldCenter[2]+dir[2]/2 * minSpacing;
339  id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer());
340  container->GetPointIds()->InsertNextId(id);
341  m_VtkCellArray->InsertNextCell(container);
342  }
343  }
344 
345  if (numDirs>1)
346  {
347  for (int i=0; i<m_GradientList.size(); i++) { pix[i] /= numDirs; }
348 
349  }
350  else if (numDirs==0)
351  {
352  if (m_SimulateBaseline) { pix = SimulateMeasurement(isoTensor, 1.0); }
353 
354  else { pix.Fill(0.0); }
355 
356  }
357 
358  m_MeanBaseline += pix[0];
359  it.Set(pix);
360  m_NumDirectionsImage->SetPixel(index, numDirs);
361  if (m_NoiseVariance>0 && m_AddNoiseFlag)
362  {
363  m_SNRImage->SetPixel(index, pix[0]/(noiseStdev*m_SignalScale));
364  }
365  ++it;
366  }
367  m_MeanBaseline /= m_ImageRegion.GetNumberOfPixels();
368  if (m_NoiseVariance>0 && m_AddNoiseFlag)
369  {
370  MITK_INFO << "Mean SNR: " << m_MeanBaseline/(noiseStdev*m_SignalScale);
371  }
372  else
373  {
374  MITK_INFO << "No noise added";
375  }
376 
377  // add rician noise
378  it.GoToBegin();
379  while(!it.IsAtEnd())
380  {
381  pix = it.Get();
382  AddNoise(pix);
383  it.Set(pix);
384  ++it;
385  }
386 
387  // generate fiber bundle
388  vtkSmartPointer<vtkPolyData> directionsPolyData = vtkSmartPointer<vtkPolyData>::New();
389  directionsPolyData->SetPoints(m_VtkPoints);
390  directionsPolyData->SetLines(m_VtkCellArray);
391  m_OutputFiberBundle = mitk::FiberBundle::New(directionsPolyData);
392  }
393  template< class TOutputScalarType >
394  double DwiPhantomGenerationFilter< TOutputScalarType >::GetTensorL2Norm(itk::DiffusionTensor3D<float>& T)
395  {
396  return sqrt(T[0]*T[0] + T[3]*T[3] + T[5]*T[5] + T[1]*T[2]*2.0 + T[2]*T[4]*2.0 + T[1]*T[4]*2.0);
397  }
398 
399 }
400 
401 #endif // __itkDwiPhantomGenerationFilter_cpp
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
MITK_INFO
#define MITK_INFO
Definition: mitkLogMacros.h:22
itk::DwiPhantomGenerationFilter
Generation of synthetic diffusion weighted images using a second rank tensor model.
Definition: itkDwiPhantomGenerationFilter.h:47
rotation
static Matrix3D rotation
Definition: mitkAffineTransformBaseTest.cpp:28
QBALL_ODFSIZE
#define QBALL_ODFSIZE
Definition: mitkQBallImage.h:38
PixelType
unsigned short PixelType
Definition: TumorProgressionMapping.cpp:55
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.