19 #include "itkImageIOBase.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" 37 #include "gdcmDictEntry.h" 38 #include "gdcmGlobal.h" 39 #include "gdcmElement.h" 40 #include "gdcmPrivateTag.h" 41 #include "gdcmStringFilter.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" 52 #include "itkNrrdImageIO.h" 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 >;
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");
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");
95 template<
typename ImageType,
typename TensorImageType>
96 typename ImageType::Pointer allocateImage(
typename TensorImageType::Pointer tenIm);
98 void GetImageInfo(std::string fName, itk::ImageIOBase::IOPixelType *pixelType, itk::ImageIOBase::IOComponentType *componentType);
100 int GetNumberOfVolumes(VolumeType::Pointer rawVol,
int nVolume,
int nSliceInVolume);
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);
105 std::vector<ImageTypeScalar3D::Pointer> ConvertDWIToScalars(std::string inputDirName, std::string maskFileName);
111 template<
typename ImageType,
typename TensorImageType>
112 typename ImageType::Pointer DTIProcessingManager::allocateImage(
typename TensorImageType::Pointer tenIm)
114 typename ImageType::Pointer outputImage = ImageType::New();
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();
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)
133 typedef itk::VectorImage<TInputPixelType, Dimensions> GradientImageType;
135 std::cout <<
"Converting DWI to vector image\n";
138 typename GradientImageType::Pointer gradIm = inputImage;
172 typename TensorReconstructionImageFilterType::GradientDirectionType vect3d;
173 typename TensorReconstructionImageFilterType::GradientDirectionContainerType::Pointer DiffusionVectors = TensorReconstructionImageFilterType::GradientDirectionContainerType::New();
175 for (
unsigned int index1 = 0; index1 < DiffusionVectorWrite.size(); index1++)
177 vnl_vector_fixed<double, 3> val = DiffusionVectorWrite[index1];
181 DiffusionVectors->InsertElement(index1, vect3d);
195 typename TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New();
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();
209 maskReader->SetFileName(maskFile);
212 maskReader->Update();
214 catch (itk::ExceptionObject & err)
216 std::cerr <<
"ExceptionObject caught !" << std::endl;
217 std::cerr << err << std::endl;
220 thresholder->SetInput(maskReader->GetOutput());
221 thresholder->Update();
222 spatialObjectMask->SetImage(thresholder->GetOutput());
226 tensorReconstructionFilter->SetGradientImage(DiffusionVectors, gradIm);
228 tensorReconstructionFilter->SetBValue(bValue);
229 tensorReconstructionFilter->Update();
232 std::string default_ext =
".nii.gz";
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;
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;
250 std::vector<ScalarImageType::Pointer> vectorOfDTIScalars;
252 typedef itk::ImageRegionConstIteratorWithIndex < TensorImageType > ConstIterType;
256 TensorImageType::Pointer tensorIm = TensorImageType::New();
257 tensorIm = tensorReconstructionFilter->GetOutput();
260 ScalarImageType::Pointer faIm = ScalarImageType::New();
263 ScalarImageType::Pointer trIm = ScalarImageType::New();
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();
274 ScalarImageType::Pointer skIm = ScalarImageType::New();
275 ScalarImageType::Pointer kuIm = ScalarImageType::New();
278 ScalarImageType::Pointer clIm = ScalarImageType::New();
279 ScalarImageType::Pointer cpIm = ScalarImageType::New();
280 ScalarImageType::Pointer csIm = ScalarImageType::New();
283 ScalarImageType::Pointer rdIm = ScalarImageType::New();
284 ScalarImageType::Pointer adIm = ScalarImageType::New();
287 ScalarImageType::Pointer r1Im = ScalarImageType::New();
288 ScalarImageType::Pointer r2Im = ScalarImageType::New();
289 ScalarImageType::Pointer r3Im = ScalarImageType::New();
292 ScalarImageType::Pointer k1Im = ScalarImageType::New();
293 ScalarImageType::Pointer k2Im = ScalarImageType::New();
294 ScalarImageType::Pointer k3Im = ScalarImageType::New();
298 faIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
301 trIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
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);
317 skIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
320 kuIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
324 clIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
325 cpIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
326 csIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
330 rdIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
331 adIm = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
335 r1Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
336 r2Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
337 r3Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
341 k1Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
342 k2Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
343 k3Im = allocateImage<ScalarImageType, TensorImageType>(tensorIm);
346 ConstIterType iter(tensorIm, tensorIm->GetLargestPossibleRegion());
347 for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
349 TensorPixelType tmp = iter.Get();
350 typename TensorPixelType::EigenValuesArrayType lambda;
351 typename TensorPixelType::EigenVectorsMatrixType vMat;
352 tmp.ComputeEigenAnalysis(lambda, vMat);
354 typename ScalarImageType::IndexType index = iter.GetIndex();
357 trIm->SetPixel(index, lambda[0] + lambda[1] + lambda[2]);
361 faIm->SetPixel(index, tmp.GetFractionalAnisotropy());
366 l1Im->SetPixel(index, lambda[2]);
367 l2Im->SetPixel(index, lambda[1]);
368 l3Im->SetPixel(index, lambda[0]);
370 typename VectorImageType::PixelType vec1(3);
371 typename VectorImageType::PixelType vec2(3);
372 typename VectorImageType::PixelType vec3(3);
374 vec1[0] = vMat[2][0];
375 vec1[1] = vMat[2][1];
376 vec1[2] = vMat[2][2];
378 vec2[0] = vMat[1][0];
379 vec2[1] = vMat[1][1];
380 vec2[2] = vMat[1][2];
382 vec3[0] = vMat[0][0];
383 vec3[1] = vMat[0][1];
384 vec3[2] = vMat[0][2];
386 v1Im->SetPixel(index, vec1);
387 v2Im->SetPixel(index, vec2);
388 v3Im->SetPixel(index, vec3);
394 ScalarPixelType m1, m3, l1, l2, l3;
398 m1 = (l1 + l2 + l3) / 3.0;
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));
405 skIm->SetPixel(index, std::pow(m3, static_cast<ScalarPixelType>(1.0 / 3.0)));
409 skIm->SetPixel(index, -1 * std::pow((-1 * m3), static_cast<ScalarPixelType>(1.0 / 3.0)));
414 skIm->SetPixel(index, static_cast<ScalarPixelType>(0));
420 ScalarPixelType m1, m4, l1, l2, l3;
424 m1 = (l1 + l2 + l3) / 3.0;
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)));
432 kuIm->SetPixel(index, static_cast<ScalarPixelType>(0));
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]);
446 clIm->SetPixel(index, 0);
447 cpIm->SetPixel(index, 0);
448 csIm->SetPixel(index, 0);
454 rdIm->SetPixel(index, (lambda[1] + lambda[0]) / 2);
455 adIm->SetPixel(index, lambda[2]);
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;
469 r1 = sqrt(3 * (std::pow(m1, 2) + m2));
470 r2 = sqrt(3 * m2 / 2 / (std::pow(m1, 2) + m2));
472 r3 = (lambda[0] * lambda[1] * lambda[2]) / std::pow(static_cast<double>(sqrt(3 * m2)), 3);
474 r1Im->SetPixel(index, r1);
475 r2Im->SetPixel(index, r2);
476 r3Im->SetPixel(index, r3);
482 ScalarPixelType m1, m2, m3;
483 ScalarPixelType k1, k2, k3;
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;
494 k3 = (lambda[0] * lambda[1] * lambda[2]) / std::pow(static_cast<double>(sqrt(3 * m2)), 3);
497 k1Im->SetPixel(index, k1);
498 k2Im->SetPixel(index, k2);
499 k3Im->SetPixel(index, k3);
503 std::cout <<
"Done Computing Scalars\n";
509 vectorOfDTIScalars.push_back(faIm);
512 vectorOfDTIScalars.push_back(trIm);
516 vectorOfDTIScalars.push_back(rdIm);
517 vectorOfDTIScalars.push_back(adIm);
520 catch (itk::ExceptionObject & excp)
522 std::cerr <<
"Exception caught - " << excp.what() <<
"\n";
524 return vectorOfDTIScalars;
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