Diffusion derivatives additions

3_HowToGuides.txt

@@ -982,7 +982,7 @@ References:
 
 /**
 \page Diffusion_Derivatives Miscellaneous: Diffusion Derivatives
-This application extracts various measurements from a Diffusion Weighted MRI scan. The exact measurements comprise i) fractional anisotropy (FA), ii) radial diffusivity (RAD), iii) axial diffusivity (AX), and iv) apparent diffusion coefficient (ADC).
+This application extracts various measurements from a Diffusion Weighted MRI scan. The exact measurements comprise i) fractional anisotropy (FA), ii) radial diffusivity (RAD), iii) axial diffusivity (AX), iv) apparent diffusion coefficient (ADC), and v) averaged b0 image (B0).
 
 This application is also available from the web on the [CBICA Image Processing Portal](https://ipp.cbica.upenn.edu/). Please see the experiment on the portal for details.
 
@@ -994,11 +994,12 @@ This application is also available from the web on the [CBICA Image Processing P
 <b> USAGE:</b>
 	-# Launch the application from "Applications" -> "Diffusion Derivatives".
 	-# Specify the paths to the input DW-MRI image, bval and bvec files with the gradient information, and the mask defining the voxels to extract the measurements.
+	-# If desired, enable registration of outputs and select a fixed image to register to.
 	-# Identify the diffusion measurements to be extracted, and the output directory.
 	-# Press the "Confirm" button.
 	-# The extracted measurements will be saved at the specified location (~5 minutes).
 \verbatim
-${CaPTk_InstallDir}/bin/ DiffusionDerivatives.exe -i C:/inputDWI.nii.gz -m C:/mask.nii.gz -b C:/input.bval -g C:/input.bvec -a 1 -f 1 -r 1 -t 1 -o C:/output
+${CaPTk_InstallDir}/bin/ DiffusionDerivatives.exe -i C:/inputDWI.nii.gz -m C:/mask.nii.gz -b C:/input.bval -g C:/input.bvec -a 1 -f 1 -r 1 -t 1 -z 1 -reg C:/template.nii.gz -o C:/output
 \endverbatim
 
 --------

src/applications/DiffusionDerivatives.cxx

@@ -17,12 +17,15 @@ int main(int argc, char **argv)
   parser.addOptionalParameter("f", "fractional", cbica::Parameter::BOOLEAN, "", "Generate the Fractional Anisotropy image (1=YES, 0=NO, 1 (Default))");
   parser.addOptionalParameter("r", "radial", cbica::Parameter::BOOLEAN, "", "Generate the Radial Diffusivity image (1=YES, 0=NO, 1 (Default))");
   parser.addOptionalParameter("t", "coefficient", cbica::Parameter::BOOLEAN, "", "Generate the Apparent Diffusion Coefficient (1=YES, 0=NO, 1 (Default))");
+  parser.addOptionalParameter("z", "b0", cbica::Parameter::BOOLEAN, "", "Generate the b0 image (1=YES, 0=NO, 1 (Default))");
+
+  parser.addOptionalParameter("reg", "registerTo", cbica::Parameter::FILE, "", "Optional file used as a fixed image for (rigid) registration of the output.");
 
   parser.addRequiredParameter("o", "output", cbica::Parameter::STRING, "", "The output directory.");
   parser.addOptionalParameter("L", "Logger", cbica::Parameter::STRING, "log file which user has write access to", "Full path to log file to store console outputs", "By default, only console output is generated");
   //parser.exampleUsage("");
-  parser.addExampleUsage("-i C:/inputDWI.nii.gz -m C:/mask.nii.gz -b C:/input.bval -g C:/input.bvec -a 1 -f 1 -r 1 -t 1 -o C:/output",
-    "Calculates diffusion derivatives for the specified input image and parameters");
+  parser.addExampleUsage("-i C:/inputDWI.nii.gz -m C:/mask.nii.gz -b C:/input.bval -g C:/input.bvec -a 1 -f 1 -r 1 -t 1 -z 1 -reg C:/registrationTemplate.nii.gz -o C:/output",
+    "Calculates diffusion derivatives and averaged b0 for the specified input image and parameters, then registers (rigidly) to the given registration template.");
   parser.addApplicationDescription("Calculates Diffusion Derivatives for a DWI image for a defined mask, bvec and bval files");
 
   // parameters to get from the command line
@@ -33,9 +36,11 @@ int main(int argc, char **argv)
   bool faPresent = 1;
   bool radPresent = 1;
   bool trPresent = 1;
+  bool bZeroPresent = 1;
+  bool registrationRequested = false;
 
   int tempPosition;
-  std::string inputFileName, inputMaskName, inputBValName, inputBVecName, outputDirectoryName;
+  std::string inputFileName, inputMaskName, inputBValName, inputBVecName, outputDirectoryName, inputRegistrationFile;
 
   if ((argc < 1) || (parser.compareParameter("u", tempPosition)))
   {
@@ -65,6 +70,11 @@ int main(int argc, char **argv)
   {
     inputBVecName = argv[tempPosition + 1];
   }
+  if (parser.compareParameter("reg", tempPosition))
+  {
+      inputRegistrationFile = argv[tempPosition + 1];
+      registrationRequested = true;
+  }
 
 
   if (parser.compareParameter("a", tempPosition))
@@ -83,6 +93,10 @@ int main(int argc, char **argv)
   {
     parser.getParameterValue("t", trPresent);
   }
+  if (parser.compareParameter("z", tempPosition))
+  {
+      parser.getParameterValue("z", bZeroPresent);
+  }
 
   if (parser.compareParameter("o", tempPosition))
   {
@@ -116,6 +130,11 @@ int main(int argc, char **argv)
     std::cout << "The input bvec file does not exist:" << inputBVecName << std::endl;
     return EXIT_FAILURE;
   }
+  if (registrationRequested && !cbica::isFile(inputRegistrationFile))
+  {
+      std::cout << "The registration fixed image file does not exist:" << inputRegistrationFile << std::endl;
+      return EXIT_FAILURE;
+  }
   if (!cbica::directoryExists(outputDirectoryName))
   {
     if (!cbica::createDirectory(outputDirectoryName))
@@ -124,17 +143,112 @@ int main(int argc, char **argv)
   }
   DiffusionDerivatives objDiffusion;
   std::vector<typename ImageTypeFloat3D::Pointer> diffusionDerivatives = objDiffusion.Run(inputFileName, inputMaskName, inputBValName, inputBVecName, outputDirectoryName);
-  std::cout << "Writing measures to the specified output directory.\n";
-
-  //fa,tr, rad , ax
-  if (faPresent== true)
-    cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[0], outputDirectoryName + "/FractionalAnisotropy.nii.gz");
-  if (trPresent == true)
-    cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[1], outputDirectoryName + "/ApparentDiffusionCoefficient.nii.gz");
-  if (radPresent == true)
-    cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[2], outputDirectoryName + "/RadialDiffusivity.nii.gz");
-  if (axPresent == true)
-    cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[3], outputDirectoryName + "/AxialDiffusivity.nii.gz");
+
+  if (registrationRequested)
+  {
+      if (faPresent == true)
+          cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[0], outputDirectoryName + "/unregistered_FractionalAnisotropy.nii.gz");
+      if (trPresent == true)
+          cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[1], outputDirectoryName + "/unregistered_ApparentDiffusionCoefficient.nii.gz");
+      if (radPresent == true)
+          cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[2], outputDirectoryName + "/unregistered_RadialDiffusivity.nii.gz");
+      if (axPresent == true)
+          cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[3], outputDirectoryName + "/unregistered_AxialDiffusivity.nii.gz");
+      if (bZeroPresent == true)
+          cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[4], outputDirectoryName + "/unregistered_b0.nii.gz");
+
+      std::string greedyPathAndDim = cbica::getExecutablePath() + "greedy" +
+#if WIN32
+          ".exe" +
+#endif
+          " -d " + std::to_string(TImageType::ImageDimension);
+      std::string fixedOptions = " -a -ia-image-centers -dof 6 -i " + inputRegistrationFile + " "; // force rigid registration only
+      std::string commandToCall;
+      std::vector<int> returnCodesFromGreedy;
+
+      // TODO (Alex): This needs to be done more dynamically
+      std::string specificParams;
+
+      // Create the intermediate transforms
+      specificParams = outputDirectoryName+"/unregistered_FractionalAnisotropy.nii.gz"+  " -o " + outputDirectoryName+"/affine_FractionalAnisotropy.mat";
+      commandToCall = greedyPathAndDim + fixedOptions + specificParams;
+      returnCodesFromGreedy.push_back(std::system(commandToCall.c_str()));
+      specificParams = outputDirectoryName + "/unregistered_ApparentDiffusionCoefficient.nii.gz" + " -o " + outputDirectoryName + "/affine_ApparentDiffusionCoefficient.mat";
+      commandToCall = greedyPathAndDim + fixedOptions + specificParams;
+      returnCodesFromGreedy.push_back(std::system(commandToCall.c_str()));
+      specificParams = outputDirectoryName + "/unregistered_RadialDiffusivity.nii.gz" + " -o " + outputDirectoryName + "/affine_RadialDiffusivity.mat";
+      commandToCall = greedyPathAndDim + fixedOptions + specificParams;
+      returnCodesFromGreedy.push_back(std::system(commandToCall.c_str()));
+      specificParams = outputDirectoryName + "/unregistered_AxialDiffusivity.nii.gz" + " -o " + outputDirectoryName + "/affine_AxialDiffusivity.mat";
+      commandToCall = greedyPathAndDim + fixedOptions + specificParams;
+      returnCodesFromGreedy.push_back(std::system(commandToCall.c_str()));
+      specificParams = outputDirectoryName + "/unregistered_b0.nii.gz" + " -o " + outputDirectoryName + "/affine_b0.mat";
+      commandToCall = greedyPathAndDim + fixedOptions + specificParams;
+      returnCodesFromGreedy.push_back(std::system(commandToCall.c_str()));
+
+      for (unsigned int i = 0; i < returnCodesFromGreedy.size(); i++)
+      {
+          if (returnCodesFromGreedy[i] != 0)
+          {
+              std::cerr << "Something went wrong when generating transforms with Greedy.\n";
+              return EXIT_FAILURE;
+          }
+      }
+
+      returnCodesFromGreedy.clear();
+      // Now run the transform and image through greedy to produce actual output...
+      std::string fixedOptions2 = " -rf " + inputRegistrationFile + " -ri LINEAR -r ";
+
+      std::string specificParams2 = outputDirectoryName + "/affine_FractionalAnisotropy.mat" +
+          " -rm " + outputDirectoryName + "/unregistered_FractionalAnisotropy.nii.gz" +
+          " " + outputDirectoryName + "/FractionalAnisotropy.nii.gz";
+      commandToCall = greedyPathAndDim + fixedOptions2 + specificParams2;
+      returnCodesFromGreedy.push_back(std::system(commandToCall.c_str()));
+      specificParams2 = outputDirectoryName + "/affine_ApparentDiffusionCoefficient.mat" +
+          " -rm " + outputDirectoryName + "/unregistered_ApparentDiffusionCoefficient.nii.gz" +
+          " " + outputDirectoryName + "/ApparentDiffusionCoefficient.nii.gz";
+      commandToCall = greedyPathAndDim + fixedOptions2 + specificParams2;
+      returnCodesFromGreedy.push_back(std::system(commandToCall.c_str()));
+      specificParams2 = outputDirectoryName + "/affine_RadialDiffusivity.mat" +
+          " -rm " + outputDirectoryName + "/unregistered_RadialDiffusivity.nii.gz" +
+          " " + outputDirectoryName + "/RadialDiffusivity.nii.gz";
+      commandToCall = greedyPathAndDim + fixedOptions2 + specificParams2;
+      returnCodesFromGreedy.push_back(std::system(commandToCall.c_str()));
+      specificParams2 = outputDirectoryName + "/affine_AxialDiffusivity.mat" +
+          " -rm " + outputDirectoryName + "/unregistered_AxialDiffusivity.nii.gz" +
+          " " + outputDirectoryName + "/AxialDiffusivity.nii.gz";
+      commandToCall = greedyPathAndDim + fixedOptions2 + specificParams2;
+      returnCodesFromGreedy.push_back(std::system(commandToCall.c_str()));
+      specificParams2 = outputDirectoryName + "/affine_b0.mat" +
+          " -rm " + outputDirectoryName + "/unregistered_b0.nii.gz" +
+          " " + outputDirectoryName + "/b0.nii.gz";
+      commandToCall = greedyPathAndDim + fixedOptions2 + specificParams2;
+      returnCodesFromGreedy.push_back(std::system(commandToCall.c_str()));
+
+      for (unsigned int i = 0; i < returnCodesFromGreedy.size(); i++)
+      {
+          if (returnCodesFromGreedy[i] != 0)
+          {
+              std::cerr << "Something went wrong when registering final output with Greedy.\n";
+              return EXIT_FAILURE;
+          }
+      }
+
+  }
+  else 
+  {
+      std::cout << "Writing measures to the specified output directory.\n";
+      if (faPresent == true)
+          cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[0], outputDirectoryName + "/FractionalAnisotropy.nii.gz");
+      if (trPresent == true)
+          cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[1], outputDirectoryName + "/ApparentDiffusionCoefficient.nii.gz");
+      if (radPresent == true)
+          cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[2], outputDirectoryName + "/RadialDiffusivity.nii.gz");
+      if (axPresent == true)
+          cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[3], outputDirectoryName + "/AxialDiffusivity.nii.gz");
+      if (bZeroPresent == true)
+          cbica::WriteImage< ImageTypeFloat3D >(diffusionDerivatives[4], outputDirectoryName + "/b0.nii.gz");
+  }
 
   std::cout << "Finished successfully.\n";
 

src/applications/DiffusionDerivatives.h

@@ -30,6 +30,9 @@ See COPYING file or https://www.med.upenn.edu/cbica/captk/license.html
 #include "itkDiffusionTensor3DReconstructionImageFilter.h"
 #include "itkImageRegionIteratorWithIndex.h"
 #include "itkBinaryThresholdImageFilter.h"
+#include "itkNaryAddImageFilter.h"
+#include "itkMaskImageFilter.h"
+#include "itkDivideImageFilter.h"
 #include "cbicaLogging.h"
 #include "CaPTkDefines.h"
 
@@ -410,8 +413,8 @@ std::vector<itk::Image<float, 3>::Pointer> DiffusionDerivatives::Run(std::string
   }
   default:
   {
-	  m_LastError = "Mask is expected to be either a unsigned char or a signed short Image. Please invectigate the image type of your supplied mask.";
-      std::cerr << "Mask is expected to be either a unsigned char or a signed short Image. Please invectigate the image type of your supplied mask" << std::endl;
+	  m_LastError = "Mask is expected to be either a unsigned char or a signed short Image. Please investigate the image type of your supplied mask.";
+      std::cerr << "Mask is expected to be either a unsigned char or a signed short Image. Please investigate the image type of your supplied mask." << std::endl;
     break;
   }
   }
@@ -591,6 +594,8 @@ std::vector<itk::Image<float, 3>::Pointer>  DiffusionDerivatives::dtiRecon(std::
   double bValue = 0;
   std::ifstream bvalIn(bvalFile.c_str());
   std::string line;
+  std::vector<int> indexesWhereBValZero;
+  int currentTimeIndex = 0;
   while (!bvalIn.eof())
   {
     getline(bvalIn, line);
@@ -599,14 +604,18 @@ std::vector<itk::Image<float, 3>::Pointer>  DiffusionDerivatives::dtiRecon(std::
 
     while (ss >> val)
     {
+      if (val == 0) //  mark these indices for b0 extraction
+         indexesWhereBValZero.push_back(currentTimeIndex);
+
       ++NumberOfGradients;
       if (val != 0 && bValue == 0)
         bValue = val;
       else if (val != 0 && bValue != val)
       {
-        std::cerr << "multiple bvalues not allowed" << std::endl;
-		m_LastError = "multiple bvalues not allowed";
+        std::cerr << "multiple unique bvalues are not currently allowed" << std::endl;
+		m_LastError = "multiple unique bvalues are not currently allowed";
       }
+      currentTimeIndex++;
     }
   }
   std::ifstream bvecIn(gradFile.c_str());
@@ -690,6 +699,9 @@ std::vector<itk::Image<float, 3>::Pointer>  DiffusionDerivatives::dtiRecon(std::
   static int writeRadAx = 1;
   static int writeSkew = 0;
   static int writeKurt = 0;
+  static int writeBZero = 1; 
+  // The way this is currently done means the calling code can break depending on if these features ever get enabled.
+  // TODO: Fix this up
 
   typedef float ScalarPixelType;
   typedef itk::DiffusionTensor3D< float > TensorPixelType;
@@ -743,6 +755,9 @@ std::vector<itk::Image<float, 3>::Pointer>  DiffusionDerivatives::dtiRecon(std::
     typename ScalarImageType::Pointer k2Im = ScalarImageType::New();
     typename ScalarImageType::Pointer k3Im = ScalarImageType::New();
 
+    // b0 image
+    typename ScalarImageType::Pointer bZeroIm = ScalarImageType::New();
+
     //Allocate all the images...
     if (writeFA)
       allocateScalarIm<ScalarImageType, TensorImageType>(faIm, tensorIm);
@@ -797,6 +812,51 @@ std::vector<itk::Image<float, 3>::Pointer>  DiffusionDerivatives::dtiRecon(std::
       allocateScalarIm<ScalarImageType, TensorImageType>(k2Im, tensorIm);
       allocateScalarIm<ScalarImageType, TensorImageType>(k3Im, tensorIm);
     }
+    if (writeBZero)
+    {
+      allocateScalarIm<ScalarImageType, TensorImageType>(bZeroIm, tensorIm);
+      typename InputImageType::RegionType inputRegion = img4D->GetLargestPossibleRegion();
+      typename InputImageType::SizeType desiredSize = inputRegion.GetSize();
+      desiredSize[3] = 0; // We'll always collapse along the 4th dimension
+
+      // Add all b0 together first to average the images
+      typedef itk::NaryAddImageFilter<ScalarImageType, ScalarImageType> AddImageFilterType; // NaryAdd makes sure we can loop easily
+      typename AddImageFilterType::Pointer nAdder = AddImageFilterType::New();
+      for (int i = 0; i < indexesWhereBValZero.size(); i++)
+      {
+          typedef itk::ExtractImageFilter<InputImageType, ScalarImageType> ExtractorType;
+          typename ExtractorType::Pointer extractor = ExtractorType::New();
+          typename InputImageType::IndexType index = inputRegion.GetIndex();
+          index[3] = indexesWhereBValZero[i]; // Extract only the current time index
+
+          typename InputImageType::RegionType targetRegion;
+          targetRegion.SetSize(desiredSize);
+          targetRegion.SetIndex(index);
+
+          extractor->SetInput(img4D);
+          extractor->SetExtractionRegion(targetRegion);
+          extractor->SetDirectionCollapseToSubmatrix();
+          extractor->Update();
+
+          nAdder->SetInput(i, extractor->GetOutput());
+      }
+      nAdder->Update();
+      // Now divide the added image by the number of images N
+      typedef itk::DivideImageFilter<ScalarImageType, ScalarImageType, ScalarImageType> DividerType;
+      typename DividerType::Pointer divider = DividerType::New();
+      divider->SetInput(nAdder->GetOutput());
+      divider->SetConstant((double)indexesWhereBValZero.size());
+      divider->Update();
+
+      // Apply mask to the now-averaged b0
+      typedef typename itk::MaskImageFilter<ScalarImageType, ImageMaskType, ScalarImageType> MaskerType;
+      typename MaskerType::Pointer masker = MaskerType::New();
+      masker->SetInput(divider->GetOutput());
+      masker->SetMaskImage(maskReader->GetOutput());
+      masker->Update();
+      bZeroIm = masker->GetOutput();
+
+    }
     //Loop though all the voxels and if compute the needed measures!
     ConstIterType iter(tensorIm, tensorIm->GetLargestPossibleRegion());
     for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
@@ -980,6 +1040,15 @@ std::vector<itk::Image<float, 3>::Pointer>  DiffusionDerivatives::dtiRecon(std::
 	vectorOfDTIScalars.push_back(NULL);
 	vectorOfDTIScalars.push_back(NULL);
 	}
+
+    if (writeBZero)
+    {
+        vectorOfDTIScalars.push_back(bZeroIm);
+    }
+    else
+    {
+        vectorOfDTIScalars.push_back(NULL);
+    }
   }
   catch (itk::ExceptionObject & excp)
   {
@@ -989,4 +1058,4 @@ std::vector<itk::Image<float, 3>::Pointer>  DiffusionDerivatives::dtiRecon(std::
   return vectorOfDTIScalars;
 }
 
-#endif
+#endif