PhilipsDWIConverter.cxx

//
// Created by Hui Xie on 12/19/16.
//
#include "PhilipsDWIConverter.h"

PhilipsDWIConverter::PhilipsDWIConverter( DWIDICOMConverterBase::DCMTKFileVector & allHeaders,
                                          DWIConverter::FileNamesContainer &       inputFileNames,
                                          const bool                               useBMatrixGradientDirections )
  : DWIDICOMConverterBase( allHeaders, inputFileNames, useBMatrixGradientDirections )
{}

PhilipsDWIConverter::~PhilipsDWIConverter() {}

void
PhilipsDWIConverter::LoadDicomDirectory()
{
  this->DWIDICOMConverterBase::LoadDicomDirectory();
  if ( !this->m_MultiSliceVolume )
  {
    this->m_NVolume = this->m_NSlice / this->m_SlicesPerVolume;
    this->m_MeasurementFrame = this->m_Volume->GetDirection();
    this->DetermineSliceOrderIS();
    this->SetDirectionsFromSliceOrder();
  }
  // single-frame file handled specially
}

void
PhilipsDWIConverter::ExtractDWIData()
{
  if ( !this->m_MultiSliceVolume )
  {
    // assume volume interleaving
    std::cout << "Number of Slices: " << this->m_NSlice << std::endl;
    std::cout << "Number of Volumes: " << this->m_NVolume << std::endl;
    std::cout << "Number of Slices in each volume: " << this->m_SlicesPerVolume << std::endl;
    // NOTE:  Philips interleaves the directions, so the all gradient directions can be
    // determined in the first "nVolume" slices which represents the first slice from each
    // of the gradient volumes.
    for ( unsigned int k = 0; k < this->m_NVolume; ++k )
    {
      std::string DiffusionDirectionality;
      double      tmpFD;
      bool        useSupplement49Definitions( false );
      if ( ( this->m_Headers[k]->GetElementCSorOB( 0x0018, 0x9075, DiffusionDirectionality, false ) == EXIT_SUCCESS ||
             this->m_Headers[k]->GetElementFD( 0x0018, 0x9087, tmpFD, false ) == EXIT_SUCCESS ) ||
           ( this->m_Headers[k]->GetElementDSorOB( 0x0018, 0x9087, tmpFD, false ) == EXIT_SUCCESS ) )
      {
        useSupplement49Definitions = true;
      }

      bool             B0FieldFound = false;
      double           b = 0.0;
      constexpr double zeroBValueTolerance = 0.1; // Implausibly small value assumed to be b0 images
      {                                           // Fill out the B-Values
        if ( useSupplement49Definitions == true )
        {
          B0FieldFound = this->m_Headers[k]->GetElementFD( 0x0018, 0x9087, b, false ) == EXIT_SUCCESS;
        }
        else
        {
          float floatB;
          if ( this->m_Headers[k]->GetElementFLorOB( 0x2001, 0x1003, floatB, false ) == EXIT_SUCCESS )
          {
            B0FieldFound = true;
          }
          if ( B0FieldFound )
          {
            b = static_cast< double >( floatB );
          }
          std::string tag;
          this->m_Headers[k]->GetElementCSorOB( 0x2001, 0x1004, tag, false );
          if ( StringContains( tag, "I" ) && ( b > zeroBValueTolerance ) )
          {
            DiffusionDirectionality = "ISOTROPIC";
          }
        }
      }

      vnl_vector_fixed< double, 3 > vect3d;
      vect3d.fill( 0.0 );
      { // Try to fill in the vect3d gradient directions
        if ( useSupplement49Definitions == true )
        {
          double doubleArray[3];
          // Use alternate method to get value out of a sequence header (Some Phillips Data).
          if ( this->m_Headers[k]->GetElementFD( 0x0018, 0x9089, 3, doubleArray, false ) == EXIT_FAILURE )
          {
            // Try alternate method.
            // std::cout << "Looking for  0018|9089 in sequence 0018,9076" << std::endl;
            // gdcm::SeqEntry *
            // DiffusionSeqEntry=this->m_Headers[k]->GetSeqEntry(0x0018,0x9076);
            itk::DCMTKSequence DiffusionSeqEntry;
            this->m_Headers[k]->GetElementSQ( 0x0018, 0x9076, DiffusionSeqEntry );
            // const unsigned int
            // n=DiffusionSeqEntry->GetNumberOfSQItems();
            unsigned int n = DiffusionSeqEntry.card();
            if ( n == 0 )
            {
              std::cout << "ERROR:  Sequence entry 0018|9076 has no items." << std::endl;
              throw;
            }
            DiffusionSeqEntry.GetElementFD( 0x0018, 0x9089, 3, doubleArray );
          }
          vect3d[0] = doubleArray[0];
          vect3d[1] = doubleArray[1];
          vect3d[2] = doubleArray[2];
          std::cout << "===== gradient orientations:" << k << " " << this->m_InputFileNames[k] << " (0018,9089) "
                    << " " << vect3d << std::endl;
        }
        else
        {
          float                  tmp[3];
          itk::DCMTKFileReader * hdr = this->m_Headers[k];
          if ( hdr->GetElementFLorOB( 0x2005, 0x10b0, tmp[0], false ) == EXIT_FAILURE ||
               hdr->GetElementFLorOB( 0x2005, 0x10b1, tmp[1], false ) == EXIT_FAILURE ||
               hdr->GetElementFLorOB( 0x2005, 0x10b2, tmp[2], false ) == EXIT_FAILURE )
          {
            // try with new philips tags.
            if ( hdr->GetElementFLorOB( 0x2005, 0x12b0, tmp[0], false ) == EXIT_FAILURE ||
                 hdr->GetElementFLorOB( 0x2005, 0x12b1, tmp[1], false ) == EXIT_FAILURE ||
                 hdr->GetElementFLorOB( 0x2005, 0x12b2, tmp[2], false ) == EXIT_FAILURE )
            {
              // last chance. GetELementFD throws exception on failure
              double tmpD[3];
              hdr->GetElementFD( 0x0018, 0x9089, 3, tmpD );
              tmp[0] = static_cast< float >( tmpD[0] );
              tmp[1] = static_cast< float >( tmpD[1] );
              tmp[2] = static_cast< float >( tmpD[2] );
            }
          }
          vect3d[0] = static_cast< double >( tmp[0] );
          vect3d[1] = static_cast< double >( tmp[1] );
          vect3d[2] = static_cast< double >( tmp[2] );
        }
      }

      if ( StringContains( DiffusionDirectionality, "ISOTROPIC" ) )
      {
        continue;
      }
      else if ( !B0FieldFound || b < zeroBValueTolerance )
      { // Deal with b0 images
        this->m_BValues.push_back( b );
        this->m_DiffusionVectors.push_back( vect3d );
      }
      else if ( StringContains( DiffusionDirectionality, "DIRECTIONAL" ) ||
                ( DiffusionDirectionality == "NONE" ) // Some new Philips data does not specify "DIRECTIONAL"
                || ( DiffusionDirectionality == "" ) )
      { // Deal with gradient direction images
        this->m_BValues.push_back( b );
        this->m_DiffusionVectors.push_back( vect3d );
      }
      else // Have no idea why we'd be here so error out
      {
        std::cout << "ERROR: DiffusionDirectionality was " << DiffusionDirectionality
                  << "  Don't know what to do with that..." << std::endl;
        throw;
      }

      std::cout << "B-value: " << b << "; diffusion direction: " << this->m_DoubleConvert( vect3d[0] ) << ", "
                << this->m_DoubleConvert( vect3d[1] ) << ", " << this->m_DoubleConvert( vect3d[2] ) << std::endl;
    }
  }
  else
  {
    // multi-frame file, everything is inside
    std::map< std::vector< double >, double > gradientDirectionAndBValue;
    std::map< std::string, int >              sliceLocations;
    std::vector< int >                        ignorePhilipsSliceMultiFrame;

    this->m_BValues.clear();
    this->m_DiffusionVectors.clear();

    itk::DCMTKSequence perFrameFunctionalGroup;
    itk::DCMTKSequence innerSeq;
    double             dwbValue;

    this->m_Headers[0]->GetElementSQ( 0x5200, 0x9230, perFrameFunctionalGroup );
    this->m_NSlice = perFrameFunctionalGroup.card();

    // have to determine if volume slices are interleaved
    std::string origins[2];
    for ( unsigned int i = 0; i < this->m_NSlice; ++i )
    {
      itk::DCMTKItem curItem;
      perFrameFunctionalGroup.GetElementItem( i, curItem );

      // index slice locations with string origin
      itk::DCMTKSequence originSeq;
      curItem.GetElementSQ( 0x0020, 0x9113, originSeq );
      std::string originString;
      originSeq.GetElementDS( 0x0020, 0x0032, originString );
      ++sliceLocations[originString];
      // save origin of first 2 slices to compare and see if the
      // volume is interleaved.
      if ( i < 2 )
      {
        origins[i] = originString;
      }

      itk::DCMTKSequence mrDiffusionSeq;
      curItem.GetElementSQ( 0x0018, 0x9117, mrDiffusionSeq );

      std::string dirValue;
      mrDiffusionSeq.GetElementCSorOB( 0x0018, 0x9075, dirValue );

      if ( StringContains( dirValue, "ISO" ) )
      {
        ignorePhilipsSliceMultiFrame.push_back( i );
      }
      else if ( StringContains( dirValue, "NONE" ) )
      {
        std::vector< double > v( 3 );
        v[0] = 0;
        v[1] = 0;
        v[2] = 0;
        unsigned int nOld = gradientDirectionAndBValue.size();
        gradientDirectionAndBValue[v] = 0;
        unsigned int nNew = gradientDirectionAndBValue.size();

        if ( nOld != nNew )
        {
          vnl_vector_fixed< double, 3 > vect3d;
          vect3d.fill( 0 );
          this->m_DiffusionVectors.push_back( vect3d );
          this->m_BValues.push_back( 0 );
        }
      }
      else
      {
        {
          bool preferredExtrationSucceeded = EXIT_FAILURE;
          try
          {
            preferredExtrationSucceeded = mrDiffusionSeq.GetElementDSorOB( 0x0018, 0x9087, dwbValue, false );
          }
          catch ( ... )
          {
            preferredExtrationSucceeded = EXIT_FAILURE;
          }
          // Try alternate method.
          if ( preferredExtrationSucceeded == EXIT_FAILURE )
          {
            mrDiffusionSeq.GetElementFD( 0x0018, 0x9087, dwbValue );
          }
        }
        itk::DCMTKSequence volSeq;
        mrDiffusionSeq.GetElementSQ( 0x0018, 0x9076, volSeq );
        double dwgVal[3];
        {
          bool preferredExtrationSucceeded = EXIT_FAILURE;
          try
          {
            preferredExtrationSucceeded = volSeq.GetElementDSorOB< double >( 0x0018, 0x9089, 3, dwgVal, false );
          }
          catch ( ... )
          {
            preferredExtrationSucceeded = EXIT_FAILURE;
          }
          // Try alternate method.
          if ( preferredExtrationSucceeded == EXIT_FAILURE )
          {
            volSeq.GetElementFD( 0x0018, 0x9089, 3, dwgVal );
          }
        }
        std::vector< double > v( 3 );
        v[0] = dwgVal[0];
        v[1] = dwgVal[1];
        v[2] = dwgVal[2];
        unsigned int nOld = gradientDirectionAndBValue.size();
        gradientDirectionAndBValue[v] = dwbValue;
        unsigned int nNew = gradientDirectionAndBValue.size();

        if ( nOld != nNew )
        {
          vnl_vector_fixed< double, 3 > vect3d;
          vect3d[0] = v[0];
          vect3d[1] = v[1];
          vect3d[2] = v[2];
          // vect3d.normalize();
          this->m_DiffusionVectors.push_back( vect3d );

          this->m_BValues.push_back( dwbValue );
        }
      }
    }
    // update values needed for (possible) de-interleave
    this->m_SlicesPerVolume = sliceLocations.size();


    std::cout << "LPS Matrix: " << std::endl << this->m_Volume->GetDirection() << std::endl;
    std::cout << "Volume Origin: " << std::endl << this->m_Volume->GetOrigin() << std::endl;
    std::cout << "Number of slices per volume: " << this->m_SlicesPerVolume << std::endl;
    std::cout << "Slice matrix size: " << this->GetRows() << " X " << this->GetCols() << std::endl;
    std::cout << "Image resolution: " << this->m_Volume->GetSpacing() << std::endl;

    this->m_MeasurementFrame = this->m_Volume->GetDirection();

    this->m_NVolume = this->m_NSlice / this->m_SlicesPerVolume;
    for ( unsigned int k2 = 0; k2 < this->m_BValues.size(); ++k2 )
    {
      std::cout << k2 << ": direction: " << this->m_DoubleConvert( this->m_DiffusionVectors[k2][0] ) << ", "
                << this->m_DoubleConvert( this->m_DiffusionVectors[k2][1] ) << ", "
                << this->m_DoubleConvert( this->m_DiffusionVectors[k2][2] ) << ", b-value: " << this->m_BValues[k2]
                << std::endl;
    }
    // de-interleave slices if the origins of the first 2 slices
    // are the same.
    if ( origins[0] == origins[1] )
    {
      // interleaved image
      DeInterleaveVolume();
    }
  }
  // deal with trailing isotropic images
  unsigned long trailingVolumes = this->m_NVolume - this->m_DiffusionVectors.size();
  if ( trailingVolumes > 0 )
  {
    std::cout << "# of Volumes " << this->m_NVolume << " # of Diffusion Vectors " << this->m_DiffusionVectors.size()
              << " Removing " << trailingVolumes << " Isotropic volumes." << std::endl;
    using ExtractImageFilterType = itk::ExtractImageFilter< Volume3DUnwrappedType, Volume3DUnwrappedType >;
    ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New();

    Volume3DUnwrappedType::RegionType desiredRegion = this->m_Volume->GetLargestPossibleRegion();
    Volume3DUnwrappedType::SizeType   desiredSize = desiredRegion.GetSize();
    desiredSize[2] -= ( trailingVolumes * this->m_SlicesPerVolume );
    desiredRegion.SetSize( desiredSize );
    extractImageFilter->SetExtractionRegion( desiredRegion );
    extractImageFilter->SetInput( this->m_Volume );
    extractImageFilter->Update();
    this->m_Volume = extractImageFilter->GetOutput();
    this->m_NVolume -= trailingVolumes;
  }
}

void
PhilipsDWIConverter::AddFlagsToDictionary()
{
  // relevant Philips private tags
  DcmDictEntry * PhilipsDictBValue = new DcmDictEntry(
    0x2001, 0x1003, DcmVR( EVR_FL ), "B Value of diffusion weighting", 1, 1, nullptr, true, "dicomtonrrd" );
  itk::DCMTKFileReader::AddDictEntry( PhilipsDictBValue );
  DcmDictEntry * PhilipsDictDiffusionDirection = new DcmDictEntry(
    0x2001, 0x1004, DcmVR( EVR_CS ), "Diffusion Gradient Direction", 1, 1, nullptr, true, "dicomtonrrd" );
  itk::DCMTKFileReader::AddDictEntry( PhilipsDictDiffusionDirection );

  DcmDictEntry * PhilipsDictDiffusionDirectionRL =
    new DcmDictEntry( 0x2005, 0x10b0, DcmVR( EVR_FL ), "Diffusion Direction R/L", 4, 4, nullptr, true, "dicomtonrrd" );
  itk::DCMTKFileReader::AddDictEntry( PhilipsDictDiffusionDirectionRL );
  DcmDictEntry * PhilipsDictDiffusionDirectionAP =
    new DcmDictEntry( 0x2005, 0x10b1, DcmVR( EVR_FL ), "Diffusion Direction A/P", 4, 4, nullptr, true, "dicomtonrrd" );
  itk::DCMTKFileReader::AddDictEntry( PhilipsDictDiffusionDirectionAP );
  DcmDictEntry * PhilipsDictDiffusionDirectionFH =
    new DcmDictEntry( 0x2005, 0x10b2, DcmVR( EVR_FL ), "Diffusion Direction F/H", 4, 4, nullptr, true, "dicomtonrrd" );
  itk::DCMTKFileReader::AddDictEntry( PhilipsDictDiffusionDirectionFH );

  // New data new uses new tags!
  DcmDictEntry * PhilipsDictDiffusionDirectionRLnew =
    new DcmDictEntry( 0x2005, 0x12b0, DcmVR( EVR_FL ), "Diffusion Direction R/L", 4, 4, nullptr, true, "dicomtonrrd" );
  itk::DCMTKFileReader::AddDictEntry( PhilipsDictDiffusionDirectionRLnew );
  DcmDictEntry * PhilipsDictDiffusionDirectionAPnew =
    new DcmDictEntry( 0x2005, 0x12b1, DcmVR( EVR_FL ), "Diffusion Direction A/P", 4, 4, nullptr, true, "dicomtonrrd" );
  itk::DCMTKFileReader::AddDictEntry( PhilipsDictDiffusionDirectionAPnew );
  DcmDictEntry * PhilipsDictDiffusionDirectionFHnew =
    new DcmDictEntry( 0x2005, 0x12b2, DcmVR( EVR_FL ), "Diffusion Direction F/H", 4, 4, nullptr, true, "dicomtonrrd" );
  itk::DCMTKFileReader::AddDictEntry( PhilipsDictDiffusionDirectionFHnew );

  // relevant Philips private tags
}