Merge pull request #1630 from ANTsX/SaveCumulativeFields

ENH: Expose forward/inverse cumulative velocity fields.
ANTsX · Dec 2, 2023 · 0ea8e53 · 0ea8e53
2 parents d1dee14 + c0be156
commit 0ea8e53
Show file tree

Hide file tree

Showing 3 changed files with 205 additions and 42 deletions.
diff --git a/Examples/KellyKapowski.cxx b/Examples/KellyKapowski.cxx
@@ -154,7 +154,8 @@ DiReCT(itk::ants::CommandLineParser * parser)
     if (verbose)
     {
       std::cout << "  Grey matter probability image not specified. "
-                << "Creating one from the segmentation image." << std::endl;
+                << "Creating one from the segmentation image using label value " 
+                << direct->GetGrayMatterLabel() << std::endl;
     }
 
     using ThresholderType = itk::BinaryThresholdImageFilter<LabelImageType, LabelImageType>;
@@ -191,8 +192,8 @@ DiReCT(itk::ants::CommandLineParser * parser)
     if (verbose)
     {
       std::cout << "  White matter probability image not specified. "
-                << "Creating one from the segmentation image." << std::endl
-                << std::endl;
+                << "Creating one from the segmentation image using label value " 
+                << direct->GetWhiteMatterLabel() << std::endl;
     }
 
     using ThresholderType = itk::BinaryThresholdImageFilter<LabelImageType, ImageType>;
@@ -366,6 +367,19 @@ DiReCT(itk::ants::CommandLineParser * parser)
     direct->AddObserver(itk::IterationEvent(), observer);
   }
 
+  /**
+   * output
+   */
+  typename itk::ants::CommandLineParser::OptionType::Pointer outputOption = parser->GetOption("output");
+  if (outputOption)
+  {
+    if (outputOption->GetFunction(0)->GetNumberOfParameters() > 1)
+    {
+    direct->SetIncludeCumulativeVelocityFields(true);
+    } 
+  }
+
+
   itk::TimeProbe timer;
   timer.Start();
   try
@@ -391,37 +405,27 @@ DiReCT(itk::ants::CommandLineParser * parser)
   /**
    * output
    */
-  typename itk::ants::CommandLineParser::OptionType::Pointer outputOption = parser->GetOption("output");
 
   if (outputOption && outputOption->GetNumberOfFunctions() > 0)
   {
     if (outputOption->GetFunction(0)->GetNumberOfParameters() == 0)
     {
-      ANTs::WriteImage<ImageType>(direct->GetOutput(), (outputOption->GetFunction(0)->GetName()).c_str());
+      ANTs::WriteImage<ImageType>(direct->GetThicknessImage(), (outputOption->GetFunction(0)->GetName()).c_str());
     }
     else if (outputOption->GetFunction(0)->GetNumberOfParameters() > 0)
     {
-      ANTs::WriteImage<ImageType>(direct->GetOutput(), (outputOption->GetFunction(0)->GetParameter()).c_str());
+      ANTs::WriteImage<ImageType>(direct->GetThicknessImage(), (outputOption->GetFunction(0)->GetParameter(0)).c_str());
       if (outputOption->GetFunction(0)->GetNumberOfParameters() > 1)
       {
-        ANTs::WriteImage<ImageType>(direct->GetOutput(1), (outputOption->GetFunction(0)->GetParameter(1)).c_str());
+        direct->GetForwardCumulativeVelocityField()->Print(std::cout, 3);
+        ANTs::WriteImage<typename DiReCTFilterType::CumulativeVelocityFieldType>(direct->GetForwardCumulativeVelocityField(), 
+                                          (outputOption->GetFunction(0)->GetParameter(1) + "ForwardVelocityField.nii.gz").c_str());
+        ANTs::WriteImage<typename DiReCTFilterType::CumulativeVelocityFieldType>(direct->GetInverseCumulativeVelocityField(), 
+                                          (outputOption->GetFunction(0)->GetParameter(1) + "InverseVelocityField.nii.gz").c_str());
       }
     }
   }
 
-  if (segmentationImageOption->GetFunction(0)->GetNumberOfParameters() == 0)
-  {
-    std::string inputFile = segmentationImageOption->GetFunction(0)->GetName();
-    ReadImage<LabelImageType>(segmentationImage, inputFile.c_str());
-  }
-  else if (segmentationImageOption->GetFunction(0)->GetNumberOfParameters() > 0)
-
-  {}
-  if (outputOption && outputOption->GetNumberOfFunctions() > 1)
-  {
-    ANTs::WriteImage<ImageType>(direct->GetOutput(1), (outputOption->GetFunction(1)->GetName()).c_str());
-  }
-
   return EXIT_SUCCESS;
 }
 
@@ -642,7 +646,7 @@ KellyKapowskiInitializeCommandLineOptions(itk::ants::CommandLineParser * parser)
     option->SetLongName("output");
     option->SetShortName('o');
     option->SetUsageOption(0, "imageFileName");
-    option->SetUsageOption(1, "[imageFileName,warpedWhiteMatterImageFileName]");
+    option->SetUsageOption(1, "[imageFileName,cumulativeVelocityFieldFileNamePrefix]");
     option->SetDescription(description);
     parser->AddOption(option);
   }

diff --git a/Utilities/itkDiReCTImageFilter.h b/Utilities/itkDiReCTImageFilter.h
@@ -54,6 +54,8 @@ class DiReCTImageFilter final : public ImageToImageFilter<TInputImage, TOutputIm
   using Pointer = SmartPointer<Self>;
   using ConstPointer = SmartPointer<const Self>;
 
+  using DataObjectPointer = DataObject::Pointer;
+
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
 
@@ -77,6 +79,8 @@ class DiReCTImageFilter final : public ImageToImageFilter<TInputImage, TOutputIm
   using VectorType = Vector<RealType, ImageDimension>;
   using DisplacementFieldType = Image<VectorType, ImageDimension>;
   using DisplacementFieldPointer = typename DisplacementFieldType::Pointer;
+  using CumulativeVelocityFieldType = Image<VectorType, ImageDimension+1>;
+  using CumulativeVelocityFieldPointer = typename CumulativeVelocityFieldType::Pointer;
   using VectorValueType = typename VectorType::ValueType;
   using PointType = typename DisplacementFieldType::PointType;
   using SparseMatrixType = vnl_sparse_matrix<RealType>;
@@ -294,6 +298,13 @@ class DiReCTImageFilter final : public ImageToImageFilter<TInputImage, TOutputIm
   itkGetConstMacro(UseBSplineSmoothing, bool);
   itkBooleanMacro(UseBSplineSmoothing);
 
+  /**
+   * Set/Get the option to include the cumulative velocity fields as output.  Default = false.
+   */
+  itkSetMacro(IncludeCumulativeVelocityFields, bool);
+  itkGetConstMacro(IncludeCumulativeVelocityFields, bool);
+  itkBooleanMacro(IncludeCumulativeVelocityFields);
+
   /**
    * Get the number of elapsed iterations.  This is a helper function for
    * reporting observations.
@@ -312,13 +323,37 @@ class DiReCTImageFilter final : public ImageToImageFilter<TInputImage, TOutputIm
    */
   itkGetConstMacro(CurrentConvergenceMeasurement, RealType);
 
+  /**
+   * Get thickness image.
+   */
+  OutputImageType * GetThicknessImage();
+
+  /**
+   * Get forward cumulative velocity
+   */
+  CumulativeVelocityFieldType * GetForwardCumulativeVelocityField();
+
+  /**
+   * Get inverse cumulative velocity
+   */
+  CumulativeVelocityFieldType * GetInverseCumulativeVelocityField();
+
+  /** Standard itk::ProcessObject subclass method. */
+  using DataObjectPointerArraySizeType = ProcessObject::DataObjectPointerArraySizeType;
+  using Superclass::MakeOutput;
+  DataObjectPointer
+  MakeOutput(DataObjectPointerArraySizeType idx) override;
+
 protected:
   DiReCTImageFilter();
   ~DiReCTImageFilter() override;
 
   void
   PrintSelf(std::ostream & os, Indent indent) const override;
 
+  void
+  GenerateOutputInformation() override;
+
   void
   ThreadedGenerateData(const RegionType &, ThreadIdType) override{};
 
@@ -359,12 +394,6 @@ class DiReCTImageFilter final : public ImageToImageFilter<TInputImage, TOutputIm
   RealImagePointer
   WarpImage(const RealImageType *, const DisplacementFieldType *);
 
-  /**
-   * Private function for inverting the deformation field.
-   */
-  void
-  InvertDisplacementField(const DisplacementFieldType *, DisplacementFieldType *);
-
   /**
    * Private function for smoothing the deformation field.
    */
@@ -415,6 +444,7 @@ class DiReCTImageFilter final : public ImageToImageFilter<TInputImage, TOutputIm
   bool                  m_UseBSplineSmoothing;
   bool                  m_UseMaskedSmoothing;
   bool                  m_RestrictDeformation;
+  bool                  m_IncludeCumulativeVelocityFields;
   SparseMatrixType      m_SparseMatrix;
   RealImagePointer      m_SparseMatrixIndexImage;
   std::vector<RealType> m_TimePoints;