diff --git a/Examples/KellyKapowski.cxx b/Examples/KellyKapowski.cxx
index 9cd53817f..96429c078 100644
--- 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
index 0b226f9ab..dece6eb24 100644
--- 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,6 +323,27 @@ 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;
@@ -319,6 +351,9 @@ class DiReCTImageFilter final : public ImageToImageFilter<TInputImage, TOutputIm
   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;
diff --git a/Utilities/itkDiReCTImageFilter.hxx b/Utilities/itkDiReCTImageFilter.hxx
index cdd9ee086..104887cf6 100644
--- a/Utilities/itkDiReCTImageFilter.hxx
+++ b/Utilities/itkDiReCTImageFilter.hxx
@@ -34,6 +34,7 @@
 #include "itkIterationReporter.h"
 #include "itkMaskedSmoothingImageFilter.h"
 #include "itkMaximumImageFilter.h"
+#include "itkPasteImageFilter.h"
 #include "itkMultiplyByConstantImageFilter.h"
 #include "itkStatisticsImageFilter.h"
 #include "itkVectorLinearInterpolateImageFunction.h"
@@ -63,26 +64,105 @@ DiReCTImageFilter<TInputImage, TOutputImage>::DiReCTImageFilter()
   , m_UseBSplineSmoothing(false)
   , m_UseMaskedSmoothing(false)
   , m_RestrictDeformation(false)
+  , m_IncludeCumulativeVelocityFields(false)
   , m_TimeSmoothingVariance(1.0)
 {
   this->m_ThicknessPriorImage = nullptr;
   this->m_SparseMatrixIndexImage = nullptr;
-  this->SetNumberOfRequiredInputs(3);
-
   this->m_TimePoints.clear();
+
+  this->SetNumberOfRequiredOutputs(3);
+
+  // thickness image
+  this->SetNthOutput(0, this->MakeOutput(0));
+
+  // forward/inverse cumulative velocity fields
+  this->SetNthOutput(1, this->MakeOutput(1));
+  this->SetNthOutput(2, this->MakeOutput(2));
 }
 
 template <typename TInputImage, typename TOutputImage>
 DiReCTImageFilter<TInputImage, TOutputImage>::~DiReCTImageFilter() = default;
 
 template <typename TInputImage, typename TOutputImage>
-void
-DiReCTImageFilter<TInputImage, TOutputImage>::GenerateData()
+auto
+DiReCTImageFilter<TInputImage, TOutputImage>::MakeOutput(
+  DataObjectPointerArraySizeType idx) -> DataObjectPointer
 {
-  if (this->m_ThicknessPriorImage)
+  if (idx > 0)
   {
-    std::cout << "Using prior thickness image." << std::endl;
+    return CumulativeVelocityFieldType::New().GetPointer();
   }
+  return Superclass::MakeOutput(idx);
+} 
+
+template <typename TInputImage, typename TOutputImage>
+auto
+DiReCTImageFilter<TInputImage, TOutputImage>::GetThicknessImage() -> OutputImageType *
+{
+  return dynamic_cast<OutputImageType *>(this->ProcessObject::GetOutput(0));
+}
+
+template <typename TInputImage, typename TOutputImage>
+auto
+DiReCTImageFilter<TInputImage, TOutputImage>::GetForwardCumulativeVelocityField() -> CumulativeVelocityFieldType *
+{
+  return dynamic_cast<CumulativeVelocityFieldType *>(this->ProcessObject::GetOutput(1));
+}
+
+template <typename TInputImage, typename TOutputImage>
+auto
+DiReCTImageFilter<TInputImage, TOutputImage>::GetInverseCumulativeVelocityField() -> CumulativeVelocityFieldType *
+{
+  return dynamic_cast<CumulativeVelocityFieldType *>(this->ProcessObject::GetOutput(2));
+}
+
+template <typename TInputImage, typename TOutputImage>
+void
+DiReCTImageFilter<TInputImage, TOutputImage>::GenerateOutputInformation()
+{
+  RealImagePointer corticalThicknessImage = this->GetThicknessImage();
+  corticalThicknessImage->CopyInformation(this->GetInput());
+  corticalThicknessImage->SetRegions(this->GetInput()->GetLargestPossibleRegion());
+
+  if (this->m_IncludeCumulativeVelocityFields)
+  {
+    CumulativeVelocityFieldPointer forwardCumulativeVelocityField = this->GetForwardCumulativeVelocityField();
+    CumulativeVelocityFieldPointer inverseCumulativeVelocityField = this->GetInverseCumulativeVelocityField();
+
+    typename CumulativeVelocityFieldType::PointType origin;
+    typename CumulativeVelocityFieldType::DirectionType direction;
+    typename CumulativeVelocityFieldType::SpacingType spacing;
+    typename CumulativeVelocityFieldType::RegionType region;
+    typename CumulativeVelocityFieldType::RegionType::SizeType size;
+    
+    for (unsigned int d = 0; d < ImageDimension; d++ )
+    {
+      origin[d] = this->GetInput()->GetOrigin()[d];
+      spacing[d] = this->GetInput()->GetSpacing()[d];
+      size[d] = this->GetInput()->GetLargestPossibleRegion().GetSize()[d];
+    }
+    origin[ImageDimension] = 0.0;
+    spacing[ImageDimension] = 1.0;
+    direction.SetIdentity();
+    size[ImageDimension] = this->m_NumberOfIntegrationPoints;
+
+    inverseCumulativeVelocityField->SetSpacing(spacing);
+    inverseCumulativeVelocityField->SetOrigin(origin);
+    inverseCumulativeVelocityField->SetDirection(direction);
+    inverseCumulativeVelocityField->SetRegions(size);
+
+    forwardCumulativeVelocityField->SetSpacing(spacing);
+    forwardCumulativeVelocityField->SetOrigin(origin);
+    forwardCumulativeVelocityField->SetDirection(direction);
+    forwardCumulativeVelocityField->SetRegions(size);
+  } 
+}
+
+template <typename TInputImage, typename TOutputImage>
+void
+DiReCTImageFilter<TInputImage, TOutputImage>::GenerateData()
+{
   this->m_CurrentGradientStep = this->m_InitialGradientStep;
 
   // Convert all input direction matrices to identities saving the original
@@ -170,12 +250,19 @@ DiReCTImageFilter<TInputImage, TOutputImage>::GenerateData()
   // Initialize fields and images.
   VectorType zeroVector(0.0);
 
-  RealImagePointer corticalThicknessImage = RealImageType::New();
-  corticalThicknessImage->CopyInformation(segmentationImage);
-  corticalThicknessImage->SetRegions(segmentationImage->GetRequestedRegion());
+  RealImagePointer corticalThicknessImage = this->GetThicknessImage();
+  corticalThicknessImage->SetDirection(segmentationImage->GetDirection());
   corticalThicknessImage->Allocate();
   corticalThicknessImage->FillBuffer(0.0);
 
+  if (this->m_IncludeCumulativeVelocityFields)
+  {
+    CumulativeVelocityFieldPointer forwardCumulativeVelocityField = this->GetForwardCumulativeVelocityField();
+    forwardCumulativeVelocityField->Allocate();
+    CumulativeVelocityFieldPointer inverseCumulativeVelocityField = this->GetInverseCumulativeVelocityField();
+    inverseCumulativeVelocityField->Allocate();
+  }
+
   DisplacementFieldPointer forwardIncrementalField = DisplacementFieldType::New();
   forwardIncrementalField->CopyInformation(segmentationImage);
   forwardIncrementalField->SetRegions(segmentationImage->GetRequestedRegion());
@@ -226,9 +313,9 @@ DiReCTImageFilter<TInputImage, TOutputImage>::GenerateData()
   // Instantiate most of the iterators all in one place
 
   ImageRegionIterator<RealImageType>         ItCorticalThicknessImage(corticalThicknessImage,
-                                                              corticalThicknessImage->GetRequestedRegion());
+                                                                      corticalThicknessImage->GetRequestedRegion());
   ImageRegionConstIterator<RealImageType>    ItGrayMatterProbabilityMap(grayMatterProbabilityImage,
-                                                                     grayMatterProbabilityImage->GetRequestedRegion());
+                                                                        grayMatterProbabilityImage->GetRequestedRegion());
   ImageRegionIterator<RealImageType>         ItHitImage(hitImage, hitImage->GetRequestedRegion());
   ImageRegionIterator<DisplacementFieldType> ItForwardIncrementalField(forwardIncrementalField,
                                                                        forwardIncrementalField->GetRequestedRegion());
@@ -470,6 +557,36 @@ DiReCTImageFilter<TInputImage, TOutputImage>::GenerateData()
 
       inverseField = inverter2->GetOutput();
       inverseField->DisconnectPipeline();
+
+      if (this->m_IncludeCumulativeVelocityFields)
+      {
+        CumulativeVelocityFieldPointer forwardCumulativeVelocityField = this->GetForwardCumulativeVelocityField();
+        CumulativeVelocityFieldPointer inverseCumulativeVelocityField = this->GetInverseCumulativeVelocityField();
+
+        using PasterType = PasteImageFilter<CumulativeVelocityFieldType, DisplacementFieldType>;
+
+        typename CumulativeVelocityFieldType::IndexType destinationIndex;
+        destinationIndex.Fill(NumericTraits<typename IndexType::IndexValueType>::Zero);
+        destinationIndex[ImageDimension] = static_cast<typename IndexType::IndexValueType>(integrationPoint);
+
+        typename PasterType::Pointer paster1 = PasterType::New();
+        paster1->SetSourceImage(integratedField);
+        paster1->SetSourceRegion(integratedField->GetLargestPossibleRegion());
+        paster1->SetDestinationImage(forwardCumulativeVelocityField);
+        paster1->SetDestinationIndex(destinationIndex);
+        paster1->Update();
+
+        this->ProcessObject::SetNthOutput(1, paster1->GetOutput());
+
+        typename PasterType::Pointer paster2 = PasterType::New();
+        paster2->SetSourceImage(inverseField);
+        paster2->SetSourceRegion(inverseField->GetLargestPossibleRegion());
+        paster2->SetDestinationImage(inverseCumulativeVelocityField);
+        paster2->SetDestinationIndex(destinationIndex);
+        paster2->Update();
+
+        this->ProcessObject::SetNthOutput(2, paster2->GetOutput());
+      }
     }
 
     // calculate the size of the solution to allow us to adjust the
@@ -624,12 +741,24 @@ DiReCTImageFilter<TInputImage, TOutputImage>::GenerateData()
 
   // Replace the identity direction with the original direction in the outputs
 
-  RealImagePointer warpedWhiteMatterProbabilityImage = this->WarpImage(whiteMatterProbabilityImage, inverseField);
-  warpedWhiteMatterProbabilityImage->SetDirection(this->GetSegmentationImage()->GetDirection());
   corticalThicknessImage->SetDirection(this->GetSegmentationImage()->GetDirection());
 
-  this->SetNthOutput(0, corticalThicknessImage);
-  this->SetNthOutput(1, warpedWhiteMatterProbabilityImage);
+  if (this->m_IncludeCumulativeVelocityFields)
+  {
+    typename CumulativeVelocityFieldType::DirectionType direction;
+    direction.SetIdentity();
+
+    for (unsigned int d = 0; d < ImageDimension; d++)
+    {
+      for (unsigned int e = 0; e < ImageDimension; e++) 
+      { 
+        direction(d, e) = this->GetSegmentationImage()->GetDirection()(d, e);
+      }
+    }
+   
+    this->GetForwardCumulativeVelocityField()->SetDirection(direction);
+    this->GetInverseCumulativeVelocityField()->SetDirection(direction);
+  }  
 }
 
 template <typename TInputImage, typename TOutputImage>