text: /*=========================================================================
   *
   *  Copyright Insight Software Consortium
   *
   *  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 "itkDiffusionTensor3DReconstructionImageFilter.h"
  #include "itkTensorFractionalAnisotropyImageFilter.h"
  #include "itkTensorRelativeAnisotropyImageFilter.h"
  #include "itkNrrdImageIO.h"
  #include "itkImageSeriesReader.h"
  #include "itkImageFileWriter.h"
  #include "itkImageRegionIterator.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;
    using ImageType = itk::VectorImage<unsigned short, 3>;
  
    itk::ImageFileReader<ImageType>::Pointer reader
      = itk::ImageFileReader<ImageType>::New();
  
    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 (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 seperate 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 you 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.
    using DWIIteratorType = itk::ImageRegionConstIterator< ImageType >;
    DWIIteratorType dwiit( img, img->GetBufferedRegion() );
    using IteratorType = itk::ImageRegionIterator< GradientImageType >;
  
    // 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;
    using FAImageType = itk::Image< RealValueType, Dimension >;
    using FAFilterType = itk::TensorFractionalAnisotropyImageFilter<
        TensorReconstructionImageFilterType::OutputImageType, FAImageType >;
  
    FAFilterType::Pointer fractionalAnisotropyFilter = FAFilterType::New();
    fractionalAnisotropyFilter->SetInput( tensorReconstructionFilter->GetOutput() );
  
    // Write the FA image
    //
    using FAWriterType = itk::ImageFileWriter< FAFilterType::OutputImageType >;
    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;
    using RAImageType = itk::Image< RealValueType, Dimension >;
    using RAFilterType = itk::TensorRelativeAnisotropyImageFilter<
        TensorReconstructionImageFilterType::OutputImageType, RAImageType >;
  
    RAFilterType::Pointer relativeAnisotropyFilter = RAFilterType::New();
    relativeAnisotropyFilter->SetInput( tensorReconstructionFilter->GetOutput() );
  
    using RAWriterType = itk::ImageFileWriter< RAFilterType::OutputImageType >;
    RAWriterType::Pointer raWriter = RAWriterType::New();
    raWriter->SetInput( relativeAnisotropyFilter->GetOutput() );
    raWriter->SetFileName(argv[5]);
    raWriter->Update();
  
    return EXIT_SUCCESS;
  }
  