ITK  5.2.0
Insight Toolkit
Examples/Filtering/DiffusionTensor3DReconstructionImageFilter.cxx
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
//
// This example shows how to use the
// DiffusionTensor3DReconstructionImageFilter to reconstruct an image of
// tensors from Diffusion weighted images. See the documentation of
// DiffusionTensor3DReconstructionImageFilter,
// TensorRelativeAnisotropyImageFilter, TensorFractionalAnisotropyImageFilter
// first.
//
// The example takes diffusion weighted images in the Nrrd format, writes out
// the gradient and reference images for the users reference and computes and
// writes out the Fractional Anisotropy and Relative Anisotropy images.
//
// The easiest way to get started is to try out the filter on a sample
// dataset.
//
// Acquiring sample datasets:
// 1. Get the DWI datasets from
// ftp://public.kitware.com/pub/namic/DTI/Data/dwi.nhdr
// ftp://public.kitware.com/pub/namic/DTI/Data/dwi.img.gz (gunzip this)
// These datasets contain a reference T1 image and 30 diffusion weighted
// images. See the nrrd header for details such as B value etc.
//
// 2. Run the example with the following args
// dwi.nhdr 80 Tensors.mhd FractionalAnisotropy.mhd
// RelativeAnisotropy.mhd 1
//
// 3. You should find 30 gradient images, 1 reference image, the FA and RA
// images
// in your working directory, which you can fire up in your favourite
// volume browser.
//
// This work is part of the National Alliance for Medical Image
// Computing (NAMIC), funded by the National Institutes of Health
// through the NIH Roadmap for Medical Research, Grant U54 EB005149.
//
// Additional documentation: For details on the Nrrd format for DTI, see
// http://wiki.na-mic.org/Wiki/index.php/NAMIC_Wiki:DTI:Nrrd_format
//
#include "itkNrrdImageIO.h"
#include <iostream>
int
main(int argc, char * argv[])
{
if (argc < 3)
{
std::cerr << "Usage: " << argv[0]
<< " NrrdFileName(.nhdr) threshold(on B0)"
<< " FAImageFileName RelativeAnisotropyFileName "
<< "[ExtractGradientAndReferenceImage from the NRRD file and "
<< "write them as images]" << std::endl;
std::cerr << "\tExample args: xt_dwi.nhdr 80 FA.mhd 1" << std::endl;
return EXIT_FAILURE;
}
constexpr unsigned int Dimension = 3;
unsigned int numberOfImages = 0;
unsigned int numberOfGradientImages = 0;
bool readb0 = false;
double b0 = 0;
using PixelType = unsigned short;
ImageType::Pointer img;
// Set the properties for NrrdReader
reader->SetFileName(argv[1]);
// Read in the nrrd data. The file contains the reference image and the
// gradient images.
try
{
reader->Update();
img = reader->GetOutput();
}
catch (const itk::ExceptionObject & ex)
{
std::cout << ex << std::endl;
return EXIT_FAILURE;
}
// Here we instantiate the DiffusionTensor3DReconstructionImageFilter class.
// The class is templated over the pixel types of the reference, gradient
// and the to be created tensor pixel's precision. (We use double here). It
// takes as input the Reference (B0 image aquired in the absence of
// diffusion sensitizing gradients), 'n' Gradient images and their
// directions and produces as output an image of tensors with pixel-type
// DiffusionTensor3D.
//
using TensorReconstructionImageFilterType = itk::
DiffusionTensor3DReconstructionImageFilter<PixelType, PixelType, double>;
// -------------------------------------------------------------------------
// Parse the Nrrd headers to get the B value and the gradient directions
// used for diffusion weighting.
//
// The Nrrd headers should look like :
// The tags specify the B value and the gradient directions. If gradient
// directions are (0,0,0), it indicates that it is a reference image.
//
// DWMRI_b-value:=800
// DWMRI_gradient_0000:= 0 0 0
// DWMRI_gradient_0001:=-1.000000 0.000000 0.000000
// DWMRI_gradient_0002:=-0.166000 0.986000 0.000000
// DWMRI_gradient_0003:=0.110000 0.664000 0.740000
// ...
//
itk::MetaDataDictionary imgMetaDictionary = img->GetMetaDataDictionary();
std::vector<std::string> imgMetaKeys = imgMetaDictionary.GetKeys();
std::vector<std::string>::const_iterator itKey = imgMetaKeys.begin();
std::string metaString;
TensorReconstructionImageFilterType::GradientDirectionType vect3d;
TensorReconstructionImageFilterType::GradientDirectionContainerType::Pointer
DiffusionVectors = TensorReconstructionImageFilterType::
GradientDirectionContainerType::New();
for (; itKey != imgMetaKeys.end(); ++itKey)
{
double x, y, z;
itk::ExposeMetaData<std::string>(imgMetaDictionary, *itKey, metaString);
if (itKey->find("DWMRI_gradient") != std::string::npos)
{
std::cout << *itKey << " ---> " << metaString << std::endl;
sscanf(metaString.c_str(), "%lf %lf %lf\n", &x, &y, &z);
vect3d[0] = x;
vect3d[1] = y;
vect3d[2] = z;
DiffusionVectors->InsertElement(numberOfImages, vect3d);
++numberOfImages;
// If the direction is 0.0, this is a reference image
if (vect3d[0] == 0.0 && vect3d[1] == 0.0 && vect3d[2] == 0.0)
{
continue;
}
++numberOfGradientImages;
}
else if (itKey->find("DWMRI_b-value") != std::string::npos)
{
std::cout << *itKey << " ---> " << metaString << std::endl;
readb0 = true;
b0 = std::stod(metaString.c_str());
}
}
std::cout << "Number of gradient images: " << numberOfGradientImages
<< " and Number of reference images: "
<< numberOfImages - numberOfGradientImages << std::endl;
if (!readb0)
{
std::cerr << "BValue not specified in header file" << std::endl;
return EXIT_FAILURE;
}
// -------------------------------------------------------------------------
// Extract the Reference and gradient images from the NRRD file
// as separate images.
//
// This is not really necessary, the filter is capable of gobbling the
// entire VectorImage (which contains the reference and the gradient image)
// and chew on it to generate the TensorImage. This is a more intuitive
// formalism since this is what is usually obtained from the Nrrd DWI
// format.
//
// Nevertheless, we go through the "unnecessary pain" of extracting the
// gradient and reference images in separate images and writing them out to
// files, so they can be fired up in your favourite volume viewer.
//
using ReferenceImageType = itk::Image<PixelType, Dimension>;
using GradientImageType = ReferenceImageType;
// A container of smart pointers to images. This container will hold
// the 'numberOfGradientImages' gradient images and the reference image.
//
std::vector<GradientImageType::Pointer> imageContainer;
// iterator to iterate over the DWI Vector image just read in.
DWIIteratorType dwiit(img, img->GetBufferedRegion());
// In this for loop, we will extract the 'n' gradient images + 1 reference
// image from the DWI Vector image.
//
for (unsigned int i = 0; i < numberOfImages; i++)
{
GradientImageType::Pointer image = GradientImageType::New();
image->CopyInformation(img);
image->SetBufferedRegion(img->GetBufferedRegion());
image->SetRequestedRegion(img->GetRequestedRegion());
image->Allocate();
IteratorType it(image, image->GetBufferedRegion());
dwiit.GoToBegin();
it.GoToBegin();
while (!it.IsAtEnd())
{
it.Set(dwiit.Get()[i]);
++it;
++dwiit;
}
imageContainer.push_back(image);
}
// If we need to write out the reference and gradient images to a file..
// Easier viewing them with a viewer than if they were in a single NRRD
// file
if ((argc > 4) && (std::stoi(argv[6])))
{
unsigned int referenceImageIndex = 0;
using GradientWriterType = itk::ImageFileWriter<GradientImageType>;
for (unsigned int i = 0; i < numberOfImages; i++)
{
GradientWriterType::Pointer gradientWriter = GradientWriterType::New();
gradientWriter->SetInput(imageContainer[i]);
char filename[100];
if (DiffusionVectors->ElementAt(i).two_norm() <=
0.0) // this is a reference image
{
std::string fn("ReferenceImage%d.mhd");
sprintf(filename, fn.c_str(), referenceImageIndex);
++referenceImageIndex;
}
else
{
std::string fn("Gradient%d.mhd");
sprintf(filename, fn.c_str(), i);
}
gradientWriter->SetFileName(filename);
gradientWriter->Update();
}
}
// -------------------------------------------------------------------------
TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter =
TensorReconstructionImageFilterType::New();
// The reference and the gradient images are conveniently provided as
// input to the DiffusionTensor3DReconstructionImageFilter as
// filter->SetGradientImage( directionsContainer, nrrdreader->GetOutput()
// );
//
// The output of the nrrdreader is a VectorImage (akin to a multicomponent
// 3D image). The nth component image is treated as a reference image if its
// corresponding gradient direction is (0,0,0). Any number of reference
// images may be specified.
//
// An alternate way to provide the inputs, when you have the reference and
// gradient images in seperate itk::Image< type, 3 > is :
//
// tensorReconstructionFilter->SetReferenceImage( image0 );
// tensorReconstructionFilter->AddGradientImage( direction1, image1 );
// tensorReconstructionFilter->AddGradientImage( direction2, image2 );
//
tensorReconstructionFilter->SetGradientImage(DiffusionVectors,
reader->GetOutput());
// This is necessary until we fix netlib/dsvdc.c
tensorReconstructionFilter->SetNumberOfWorkUnits(1);
tensorReconstructionFilter->SetBValue(b0);
tensorReconstructionFilter->SetThreshold(
static_cast<TensorReconstructionImageFilterType::ReferencePixelType>(
std::stod(argv[2])));
tensorReconstructionFilter->Update();
// -------------------------------------------------------------------------
// Write out the image of tensors. This code snippet goes to show that you
// can use itk::ImageFileWriter to write an image of tensors.
//
using TensorWriterType = itk::ImageFileWriter<
TensorReconstructionImageFilterType::OutputImageType>;
TensorWriterType::Pointer tensorWriter = TensorWriterType::New();
tensorWriter->SetFileName(argv[3]);
tensorWriter->SetInput(tensorReconstructionFilter->GetOutput());
tensorWriter->Update();
// -------------------------------------------------------------------------
// Now that we have the image of tensors, we may use one of the many tensor
// filters in ITK. Below, we use the TensorFractionalAnisotropyImageFilter
// to compute the FA.
//
using TensorPixelType =
TensorReconstructionImageFilterType::TensorPixelType;
using RealValueType = TensorPixelType::RealValueType;
TensorReconstructionImageFilterType::OutputImageType,
FAImageType>;
FAFilterType::Pointer fractionalAnisotropyFilter = FAFilterType::New();
fractionalAnisotropyFilter->SetInput(
tensorReconstructionFilter->GetOutput());
// Write the FA image
//
FAWriterType::Pointer faWriter = FAWriterType::New();
faWriter->SetInput(fractionalAnisotropyFilter->GetOutput());
faWriter->SetFileName(argv[4]);
faWriter->Update();
// Compute and write the Relative Anisotropy
//
using TensorPixelType =
TensorReconstructionImageFilterType::TensorPixelType;
using RealValueType = TensorPixelType::RealValueType;
TensorReconstructionImageFilterType::OutputImageType,
RAImageType>;
RAFilterType::Pointer relativeAnisotropyFilter = RAFilterType::New();
relativeAnisotropyFilter->SetInput(tensorReconstructionFilter->GetOutput());
RAWriterType::Pointer raWriter = RAWriterType::New();
raWriter->SetInput(relativeAnisotropyFilter->GetOutput());
raWriter->SetFileName(argv[5]);
raWriter->Update();
return EXIT_SUCCESS;
}
itkTensorRelativeAnisotropyImageFilter.h
itk::ImageFileReader::New
static Pointer New()
itk::TensorFractionalAnisotropyImageFilter
Computes the Fractional Anisotropy for every pixel of a input tensor image.
Definition: itkTensorFractionalAnisotropyImageFilter.h:75
itkDiffusionTensor3DReconstructionImageFilter.h
itk::VectorImage
Templated n-dimensional vector image class.
Definition: itkImageAlgorithm.h:29
itk::SmartPointer< Self >
itkImageRegionIterator.h
itk::TensorRelativeAnisotropyImageFilter
Computes the Relative Anisotropy for every pixel of a input tensor image.
Definition: itkTensorRelativeAnisotropyImageFilter.h:75
itk::ImageRegionIterator
A multi-dimensional iterator templated over image type that walks a region of pixels.
Definition: itkImageRegionIterator.h:80
itkImageSeriesReader.h
itk::MetaDataDictionary
Provides a mechanism for storing a collection of arbitrary data types.
Definition: itkMetaDataDictionary.h:53
itk::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:88
itkNrrdImageIO.h
itkImageFileWriter.h
itkTensorFractionalAnisotropyImageFilter.h
itk::ImageRegionConstIterator< ImageType >
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itk::MetaDataDictionary::GetKeys
std::vector< std::string > GetKeys() const