CBICA Toolkit  1.0.0
All Classes Files Functions Variables Typedefs Enumerations Pages
cbicaDTIProcessingManager.h
Go to the documentation of this file.
1 /**
2 \file cbicaDTIProcessingManager.h
3 
4 \brief File that holds the DTIProcessingManager class
5 
6 https://www.med.upenn.edu/sbia/software/// <br>
7 software@cbica.upenn.edu
8 
9 Copyright (c) 2018 University of Pennsylvania. All rights reserved. <br>
10 See COPYING file or https://www.med.upenn.edu/cbica/software-agreement.html
11 
12 */
13 
14 #pragma once
15 
16 #include "string.h"
17 //#include <type_traits>
18 #include <algorithm>
19 #include "itkImageIOBase.h"
20 
21 #include "itkImage.h"
22 #include "itkVectorImage.h"
23 #include "itkImageFileReader.h"
24 #include "itkImageSeriesReader.h"
25 #include "itkImageFileWriter.h"
26 #include "itkCastImageFilter.h"
27 #include "itkGDCMImageIO.h"
28 #include "itkGDCMSeriesFileNames.h"
29 #include "itkExtractImageFilter.h"
30 #include "itkVectorImage.h"
32 #include "itkBinaryThresholdImageFilter.h"
33 #include "itkDiffusionTensor3D.h"
34 #include "gdcmDictEntry.h"
35 #include "gdcmDict.h" // access to dictionary
36 #include "gdcmFile.h" // access to dictionary
37 #include "gdcmDictEntry.h" // access to dictionary
38 #include "gdcmGlobal.h" // access to dictionary
39 #include "gdcmElement.h"
40 #include "gdcmPrivateTag.h"
41 #include "gdcmStringFilter.h"
42 #include "cbicaUtilities.h"
43 #include "itkImageMaskSpatialObject.h"
44 #include "itkBinaryThresholdImageFilter.h"
45 #include "itkNumericTraits.h"
46 #include "itkRescaleIntensityImageFilter.h"
47 #include "itkDiffusionTensor3D.h"
48 #include "itkTensorFractionalAnisotropyImageFilter.h"
49 #include "itkTensorRelativeAnisotropyImageFilter.h"
50 #include "itkDiscreteGaussianImageFilter.h"
51 
52 #include "itkNrrdImageIO.h"
53 
54 //#include "CAPTk.h"
55 
56 const int Dimensions = 3;
57 typedef short PixelValueType;
58 typedef itk::Image<PixelValueType, Dimensions> VolumeType;
59 typedef itk::ImageSeriesReader<VolumeType> ReaderType;
60 typedef itk::GDCMImageIO ImageIOType;
61 typedef itk::GDCMSeriesFileNames InputNamesGeneratorType;
62 typedef itk::MetaDataDictionary DictionaryType;
63 typedef itk::VectorImage<PixelValueType, Dimensions> VectorImageType;
64 using ImageTypeScalar3D = itk::Image< float, 3 >;
65 
66 // relevant GE private tags
67 const gdcm::DictEntry GEDictBValue("0x0043", "0x1039", gdcm::VR::IS, gdcm::VM::VM1, "B Value of diffusion weighting");
68 const gdcm::DictEntry GEDictXGradient("0x0019", "0x10bb", gdcm::VR::DS, gdcm::VM::VM1, "X component of gradient direction");
69 const gdcm::DictEntry GEDictYGradient("0x0019", "0x10bc", gdcm::VR::DS, gdcm::VM::VM1, "Y component of gradient direction");
70 const gdcm::DictEntry GEDictZGradient("0x0019", "0x10bd", gdcm::VR::DS, gdcm::VM::VM1, "Z component of gradient direction");
71 
72 // relevant Siemens private tags
73 const gdcm::DictEntry SiemensMosiacParameters("0x0051", "0x100b", gdcm::VR::IS, gdcm::VM::VM1, "Mosiac Matrix Size");
74 const gdcm::DictEntry SiemensDictNMosiac("0x0019", "0x100a", gdcm::VR::US, gdcm::VM::VM1, "Number of Images In Mosaic");
75 const gdcm::DictEntry SiemensDictBValue("0x0019", "0x100c", gdcm::VR::IS, gdcm::VM::VM1, "B Value of diffusion weighting");
76 const gdcm::DictEntry SiemensDictDiffusionDirection("0x0019", "0x100e", gdcm::VR::FD, gdcm::VM::VM3, "Diffusion Gradient Direction");
77 const gdcm::DictEntry SiemensDictDiffusionMatrix("0x0019", "0x1027", gdcm::VR::FD, gdcm::VM::VM6, "Diffusion Matrix");
78 
79 
80 namespace cbica
81 {
82  /**
83  \class DTIProcessingManager
84 
85  \brief A small description of the class
86 
87  A detailed description.
88  }
89  */
91  {
92  public:
95  template<typename ImageType, typename TensorImageType>
96  typename ImageType::Pointer allocateImage(typename TensorImageType::Pointer tenIm);
97 
98  void GetImageInfo(std::string fName, itk::ImageIOBase::IOPixelType *pixelType, itk::ImageIOBase::IOComponentType *componentType);
99 
100  int GetNumberOfVolumes(VolumeType::Pointer rawVol, int nVolume, int nSliceInVolume);
101 
102  template < typename TInputPixelType, typename TMaskPixelType, typename TOutputTensorCompType>
103  std::vector<ImageTypeScalar3D::Pointer> dtiRecon(VectorImageType::Pointer inputImage, std::string maskFile, int verbose, int inputIsVectorImage, std::vector< vnl_vector_fixed<double, 3> > diffuionVector, double bValue);
104 
105  std::vector<ImageTypeScalar3D::Pointer> ConvertDWIToScalars(std::string inputDirName, std::string maskFileName);
106 
107 
108  };
109 
110 
111  template<typename ImageType, typename TensorImageType>
112  typename ImageType::Pointer DTIProcessingManager::allocateImage(typename TensorImageType::Pointer tenIm)
113  {
114  typename ImageType::Pointer outputImage = ImageType::New();
115 
116  outputImage->SetOrigin(tenIm->GetOrigin());
117  outputImage->SetSpacing(tenIm->GetSpacing());
118  outputImage->SetDirection(tenIm->GetDirection());
119  outputImage->SetLargestPossibleRegion(tenIm->GetLargestPossibleRegion());
120  outputImage->SetRequestedRegion(tenIm->GetRequestedRegion());
121  outputImage->SetBufferedRegion(tenIm->GetBufferedRegion());
122  outputImage->Allocate();
123 
124  return outputImage;
125  }
126  //--------------------------------------------------------------
127  //----------------------------------------------------------
128  template < typename TInputPixelType, typename TMaskPixelType, typename TOutputTensorCompType>
129  std::vector<itk::Image<ImageTypeScalar3D::PixelType, Dimensions>::Pointer> DTIProcessingManager::dtiRecon(VectorImageType::Pointer inputImage, std::string maskFile, int verbose, int inputIsVectorImage, std::vector< vnl_vector_fixed<double, 3> > DiffusionVectorWrite, double bValue)
130  {
131  //bool readb0 = false;
132  //double b0 = 0;
133  typedef itk::VectorImage<TInputPixelType, Dimensions> GradientImageType;
134  //typedef itk::Image<TInputPixelType, 4> InputImageType;
135  std::cout << "Converting DWI to vector image\n";
136 
137  //Set up the gradient image size
138  typename GradientImageType::Pointer gradIm = inputImage;
140  // -------------------------------------------------------------------------
141 
142  //int NumberOfGradients = DiffusionVectorWrite.size();
143  //double bValue = 0;
144  //std::ifstream bvalIn(bvalFile.c_str());
145  //std::string line;
146  //while (!bvalIn.eof())
147  //{
148  // getline(bvalIn, line);
149  // double val;
150  // std::stringstream ss(line);
151 
152  // while (ss >> val)
153  // {
154  // ++NumberOfGradients;
155  // if (val != 0 && bValue == 0)
156  // bValue = val;
157  // else if (val != 0 && bValue != val)
158  // {
159  // std::cerr << "multiple bvalues not allowed" << std::endl;
160  // //return EXIT_FAILURE;
161  // }
162  // }
163  //}
164  //std::ifstream bvecIn(gradFile.c_str());
165  //getline(bvecIn, line);
166  //std::stringstream Xss(line);
167  //getline(bvecIn, line);
168  //std::stringstream Yss(line);
169  //getline(bvecIn, line);
170  //std::stringstream Zss(line);
171 
172  typename TensorReconstructionImageFilterType::GradientDirectionType vect3d;
173  typename TensorReconstructionImageFilterType::GradientDirectionContainerType::Pointer DiffusionVectors = TensorReconstructionImageFilterType::GradientDirectionContainerType::New();
174 
175  for (unsigned int index1 = 0; index1 < DiffusionVectorWrite.size(); index1++)
176  {
177  vnl_vector_fixed<double, 3> val = DiffusionVectorWrite[index1];
178  vect3d[0] = val[0];
179  vect3d[1] = val[1];
180  vect3d[2] = val[2];
181  DiffusionVectors->InsertElement(index1, vect3d);
182  }
183 
184  //int counter = 0;
185  //double x, y, z;
186  //while (Xss >> x)
187  //{
188  // Yss >> y;
189  // Zss >> z;
190  // vect3d[0] = x; vect3d[1] = y; vect3d[2] = z;
191  // DiffusionVectors->InsertElement(counter, vect3d);
192  // ++counter;
193  //}
194  // -------------------------------------------------------------------------
195  typename TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New();
196 
197  typedef itk::Image<TMaskPixelType, Dimensions > ImageMaskType;
198  typedef itk::ImageMaskSpatialObject< Dimensions > MaskType;
199  typename MaskType::Pointer spatialObjectMask = MaskType::New();
200  typedef itk::ImageFileReader<ImageMaskType> MaskReaderType;
201  typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
202  typedef itk::BinaryThresholdImageFilter<ImageMaskType, typename MaskType::ImageType> ThresholderType;
203  typename ThresholderType::Pointer thresholder = ThresholderType::New();
204  thresholder->SetOutsideValue(itk::NumericTraits< MaskType::ImageType::PixelType>::Zero);
205  thresholder->SetInsideValue(itk::NumericTraits< MaskType::ImageType::PixelType>::One);
206  thresholder->SetLowerThreshold(itk::NumericTraits<unsigned short>::One);
207  thresholder->InPlaceOn();
208 
209  maskReader->SetFileName(maskFile);
210  try
211  {
212  maskReader->Update();
213  }
214  catch (itk::ExceptionObject & err)
215  {
216  std::cerr << "ExceptionObject caught !" << std::endl;
217  std::cerr << err << std::endl;
218  //return EXIT_FAILURE;
219  }
220  thresholder->SetInput(maskReader->GetOutput());
221  thresholder->Update();
222  spatialObjectMask->SetImage(thresholder->GetOutput());
223  //tensorReconstructionFilter->SetMaskSpatialObject(spatialObjectMask);
224 
225  //---------------------------------------------------------------------------------
226  tensorReconstructionFilter->SetGradientImage(DiffusionVectors, gradIm);
227  //tensorReconstructionFilter->SetNumberOfThreads(1);
228  tensorReconstructionFilter->SetBValue(bValue);
229  tensorReconstructionFilter->Update();
230 
231  //--------------------------computing scalars----------------------------
232  std::string default_ext = ".nii.gz";
233  // perhaps the following should be replaced with bools
234  static int writeFA = 1;
235  static int writeTR = 1;
236  static int writeEign = 0;
237  static int writeGeo = 0;
238  static int writeGordR = 0;
239  static int writeGordK = 0;
240  static int writeRadAx = 1;
241  static int writeSkew = 0;
242  static int writeKurt = 0;
243 
244  typedef float ScalarPixelType;
245  typedef itk::DiffusionTensor3D< float > TensorPixelType;
246  typedef itk::Image< float, Dimensions > ScalarImageType;
247  typedef itk::VectorImage<float, Dimensions > VectorImageType;
248  typedef itk::Image< TensorPixelType, Dimensions > TensorImageType;
249 
250  std::vector<ScalarImageType::Pointer> vectorOfDTIScalars;
251 
252  typedef itk::ImageRegionConstIteratorWithIndex < TensorImageType > ConstIterType;
253 
254  try
255  {
256  TensorImageType::Pointer tensorIm = TensorImageType::New();
257  tensorIm = tensorReconstructionFilter->GetOutput();
258  // Allocate each image that we will want to make..
259  //FA
260  ScalarImageType::Pointer faIm = ScalarImageType::New();
261 
262  //TR
263  ScalarImageType::Pointer trIm = ScalarImageType::New();
264 
265  //Eigensys
266  ScalarImageType::Pointer l1Im = ScalarImageType::New();
267  ScalarImageType::Pointer l2Im = ScalarImageType::New();
268  ScalarImageType::Pointer l3Im = ScalarImageType::New();
269  VectorImageType::Pointer v1Im = VectorImageType::New();
270  VectorImageType::Pointer v2Im = VectorImageType::New();
271  VectorImageType::Pointer v3Im = VectorImageType::New();
272 
273  //Skewness & Kurtosis
274  ScalarImageType::Pointer skIm = ScalarImageType::New();
275  ScalarImageType::Pointer kuIm = ScalarImageType::New();
276 
277  //Geometric Features
278  ScalarImageType::Pointer clIm = ScalarImageType::New();
279  ScalarImageType::Pointer cpIm = ScalarImageType::New();
280  ScalarImageType::Pointer csIm = ScalarImageType::New();
281 
282  //Radial Axial
283  ScalarImageType::Pointer rdIm = ScalarImageType::New();
284  ScalarImageType::Pointer adIm = ScalarImageType::New();
285 
286  //Gordons R features
287  ScalarImageType::Pointer r1Im = ScalarImageType::New();
288  ScalarImageType::Pointer r2Im = ScalarImageType::New();
289  ScalarImageType::Pointer r3Im = ScalarImageType::New();
290 
291  //Gordons K features
292  ScalarImageType::Pointer k1Im = ScalarImageType::New();
293  ScalarImageType::Pointer k2Im = ScalarImageType::New();
294  ScalarImageType::Pointer k3Im = ScalarImageType::New();
295 
296  //Allocate all the images...
297  if (writeFA)
298  faIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
299 
300  if (writeTR)
301  trIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
302 
303  if (writeEign)
304  {
305  l1Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
306  l2Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
307  l3Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
308  v1Im->SetVectorLength(Dimensions);
309  v2Im->SetVectorLength(Dimensions);
310  v3Im->SetVectorLength(Dimensions);
311  v1Im = allocateImage<VectorImageType, TensorImageType>(tensorIm);
312  v2Im = allocateImage<VectorImageType, TensorImageType>(tensorIm);
313  v3Im = allocateImage<VectorImageType, TensorImageType>(tensorIm);
314  }
315 
316  if (writeSkew)
317  skIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
318 
319  if (writeKurt)
320  kuIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
321 
322  if (writeGeo)
323  {
324  clIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
325  cpIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
326  csIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
327  }
328  if (writeRadAx)
329  {
330  rdIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
331  adIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
332  }
333  if (writeGordR)
334  {
335  r1Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
336  r2Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
337  r3Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
338  }
339  if (writeGordK)
340  {
341  k1Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
342  k2Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
343  k3Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
344  }
345  //Loop though all the voxels and if compute the needed measures!
346  ConstIterType iter(tensorIm, tensorIm->GetLargestPossibleRegion());
347  for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
348  {
349  TensorPixelType tmp = iter.Get();
350  typename TensorPixelType::EigenValuesArrayType lambda;
351  typename TensorPixelType::EigenVectorsMatrixType vMat;
352  tmp.ComputeEigenAnalysis(lambda, vMat);
353 
354  typename ScalarImageType::IndexType index = iter.GetIndex();
355 
356  if (writeTR)
357  trIm->SetPixel(index, lambda[0] + lambda[1] + lambda[2]);
358 
359  if (writeFA)
360  {
361  faIm->SetPixel(index, tmp.GetFractionalAnisotropy());
362  }
363 
364  if (writeEign)
365  {
366  l1Im->SetPixel(index, lambda[2]);
367  l2Im->SetPixel(index, lambda[1]);
368  l3Im->SetPixel(index, lambda[0]);
369 
370  typename VectorImageType::PixelType vec1(3);
371  typename VectorImageType::PixelType vec2(3);
372  typename VectorImageType::PixelType vec3(3);
373 
374  vec1[0] = vMat[2][0];
375  vec1[1] = vMat[2][1];
376  vec1[2] = vMat[2][2];
377 
378  vec2[0] = vMat[1][0];
379  vec2[1] = vMat[1][1];
380  vec2[2] = vMat[1][2];
381 
382  vec3[0] = vMat[0][0];
383  vec3[1] = vMat[0][1];
384  vec3[2] = vMat[0][2];
385 
386  v1Im->SetPixel(index, vec1);
387  v2Im->SetPixel(index, vec2);
388  v3Im->SetPixel(index, vec3);
389 
390  }
391 
392  if (writeSkew)
393  {
394  ScalarPixelType m1, m3, l1, l2, l3;
395  l1 = abs(lambda[0]);
396  l2 = abs(lambda[1]);
397  l3 = abs(lambda[2]);
398  m1 = (l1 + l2 + l3) / 3.0;
399 
400  if (m1 > 0)
401  {
402  m3 = (std::pow(l1 - m1, 3) + std::pow(l2 - m1, 3) + std::pow(l3 - m1, 3)) / (std::pow(l1, 3) + std::pow(l2, 3) + std::pow(l3, 3));
403  if (m3 > 0)
404  {
405  skIm->SetPixel(index, std::pow(m3, static_cast<ScalarPixelType>(1.0 / 3.0)));
406  }
407  else
408  {
409  skIm->SetPixel(index, -1 * std::pow((-1 * m3), static_cast<ScalarPixelType>(1.0 / 3.0)));
410  }
411  }
412  else
413  {
414  skIm->SetPixel(index, static_cast<ScalarPixelType>(0));
415  }
416  }
417 
418  if (writeKurt)
419  {
420  ScalarPixelType m1, m4, l1, l2, l3;
421  l1 = abs(lambda[0]);
422  l2 = abs(lambda[1]);
423  l3 = abs(lambda[2]);
424  m1 = (l1 + l2 + l3) / 3.0;
425  if (m1 > 0)
426  {
427  m4 = (std::pow(l1 - m1, 4) + std::pow(l2 - m1, 4) + std::pow(l3 - m1, 4)) / (std::pow(l1, 4) + std::pow(l2, 4) + std::pow(l3, 4));
428  kuIm->SetPixel(index, std::pow(m4, static_cast<ScalarPixelType>(1.0 / 4.0)));
429  }
430  else
431  {
432  kuIm->SetPixel(index, static_cast<ScalarPixelType>(0));
433  }
434  }
435 
436  if (writeGeo)
437  {
438  if (lambda[2] > 0)
439  {
440  clIm->SetPixel(index, (lambda[2] - lambda[1]) / lambda[2]);
441  cpIm->SetPixel(index, (lambda[1] - lambda[0]) / lambda[2]);
442  csIm->SetPixel(index, lambda[0] / lambda[2]);
443  }
444  else
445  {
446  clIm->SetPixel(index, 0);
447  cpIm->SetPixel(index, 0);
448  csIm->SetPixel(index, 0);
449  }
450  }
451 
452  if (writeRadAx)
453  {
454  rdIm->SetPixel(index, (lambda[1] + lambda[0]) / 2);
455  adIm->SetPixel(index, lambda[2]);
456  }
457 
458  if (writeGordR)
459  {
460  //Compute the moments...
461  ScalarPixelType m1, m2, m3;
462  ScalarPixelType r1, r2, r3;
463  m1 = (lambda[0] + lambda[1] + lambda[2]) / 3.0;
464  m2 = (std::pow(lambda[0] - m1, 2) + std::pow(lambda[1] - m1, 2)
465  + std::pow(lambda[2] - m1, 2)) / 3.0;
466  m3 = (std::pow(lambda[0] - m1, 3) + std::pow(lambda[1] - m1, 3)
467  + std::pow(lambda[2] - m1, 3)) / 3.0;
468 
469  r1 = sqrt(3 * (std::pow(m1, 2) + m2));
470  r2 = sqrt(3 * m2 / 2 / (std::pow(m1, 2) + m2));
471  // r3 = sqrt(2) * m3 / std::pow(static_cast<double>(m2), static_cast<double>(1.5));
472  r3 = (lambda[0] * lambda[1] * lambda[2]) / std::pow(static_cast<double>(sqrt(3 * m2)), 3);
473 
474  r1Im->SetPixel(index, r1);
475  r2Im->SetPixel(index, r2);
476  r3Im->SetPixel(index, r3);
477  }
478 
479  if (writeGordK)
480  {
481  //Compute the moments...
482  ScalarPixelType m1, m2, m3;
483  ScalarPixelType k1, k2, k3;
484 
485  m1 = (lambda[0] + lambda[1] + lambda[2]) / 3.0;
486  m2 = (std::pow(lambda[0] - m1, 2) + std::pow(lambda[1] - m1, 2)
487  + std::pow(lambda[2] - m1, 2)) / 3.0;
488  m3 = (std::pow(lambda[0] - m1, 3) + std::pow(lambda[1] - m1, 3)
489  + std::pow(lambda[2] - m1, 3)) / 3.0;
490 
491  k1 = 3 * m1;
492  k2 = sqrt(3 * m2);
493  // k3 = sqrt(2) * m3 / std::pow(static_cast<double>(m2), static_cast<double>(1.5));
494  k3 = (lambda[0] * lambda[1] * lambda[2]) / std::pow(static_cast<double>(sqrt(3 * m2)), 3);
495 
496 
497  k1Im->SetPixel(index, k1);
498  k2Im->SetPixel(index, k2);
499  k3Im->SetPixel(index, k3);
500  }
501  }
502 
503  std::cout << "Done Computing Scalars\n";
504 
505  //typedef itk::ImageFileWriter< ScalarImageType > ScalarWriterType;
506  //typedef itk::ImageFileWriter< VectorImageType > VectorWriterType;
507 
508  if (writeFA)
509  vectorOfDTIScalars.push_back(faIm);
510 
511  if (writeTR)
512  vectorOfDTIScalars.push_back(trIm);
513 
514  if (writeRadAx)
515  {
516  vectorOfDTIScalars.push_back(rdIm);
517  vectorOfDTIScalars.push_back(adIm);
518  }
519  }
520  catch (itk::ExceptionObject & excp)
521  {
522  std::cerr << "Exception caught - " << excp.what() << "\n";
523  }
524  return vectorOfDTIScalars;
525  }
526  //------------------------------------------------------------------
527 
528 
529 
530 }
Declaration of DiffusionTensor3DReconstructionImageFilter.
This class takes as input one or more reference image (acquired in the absence of diffusion sensitizi...
Definition: itkDiffusionTensor3DReconstructionImageFilter.h:133
Some basic utility functions.
A small description of the class.
Definition: cbicaDTIProcessingManager.h:90