From 75bee914c7393764f9a2e636fcf71e9a4983da16 Mon Sep 17 00:00:00 2001
From: Jose M Pozo <72624361+josempozo@users.noreply.github.com>
Date: Tue, 8 Feb 2022 15:52:10 +0100
Subject: [PATCH] STYLE: Clarifying logic in
 ComputeJacobianWithRespectToPositionInternal

The implementation syntax of the spatial derivative estimation from finite differences have been simplified to make the logic more evident. In addition the robustness for displacement fields with non-zero initial index has been improved.
---
 .../include/itkDisplacementFieldTransform.hxx | 75 ++++++-------------
 1 file changed, 23 insertions(+), 52 deletions(-)

diff --git a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldTransform.hxx b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldTransform.hxx
index 573e646c988..334008b8991 100644
--- a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldTransform.hxx
+++ b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldTransform.hxx
@@ -225,21 +225,12 @@ DisplacementFieldTransform<TParametersValueType, NDimensions>::ComputeJacobianWi
 
   typename DisplacementFieldType::SizeType    size = this->m_DisplacementField->GetLargestPossibleRegion().GetSize();
   typename DisplacementFieldType::IndexType   startingIndex = this->m_DisplacementField->GetLargestPossibleRegion().GetIndex();
+  typename DisplacementFieldType::IndexType   upperIndex = this->m_DisplacementField->GetLargestPossibleRegion().GetUpperIndex();
   typename DisplacementFieldType::SpacingType spacing = this->m_DisplacementField->GetSpacing();
 
-  IndexType ddrindex;
-  IndexType ddlindex;
-  IndexType difIndex[NDimensions][2];
-
-  // Index offset
-  unsigned int posoff = NumericTraits<unsigned int>::OneValue();
-
-  // Apace between indices
+  // Space between indices
   TParametersValueType space = NumericTraits<TParametersValueType>::OneValue();
 
-  // Minimum distance between neighbors
-  IndexValueType mindist = NumericTraits<IndexValueType>::OneValue();
-
   // Flag indicating a valid location for Jacobian calculation
   bool isValidJacobianCalcLocat = true;
 
@@ -248,15 +239,10 @@ DisplacementFieldTransform<TParametersValueType, NDimensions>::ComputeJacobianWi
   dPixSign = doInverseJacobian ? -dPixSign : dPixSign;
   for (unsigned int row = 0; row < NDimensions; ++row)
   {
-    IndexValueType dist = index[row] - startingIndex[row];
-    if (dist < mindist)
-    {
-      isValidJacobianCalcLocat = false;
-    }
-    dist = static_cast<IndexValueType>(size[row]) - dist;
-    if (dist < mindist)
+    if (index[row] <= startingIndex[row] || index[row] >= upperIndex[row])
     {
       isValidJacobianCalcLocat = false;
+      break;
     }
   }
 
@@ -266,51 +252,36 @@ DisplacementFieldTransform<TParametersValueType, NDimensions>::ComputeJacobianWi
     // do manually here
     for (unsigned int row = 0; row < NDimensions; ++row)
     {
-      difIndex[row][0] = index;
-      difIndex[row][1] = index;
-      ddrindex = index;
-      ddlindex = index;
-      if ((int)index[row] < (int)(size[row] - 2))
+      IndexType difIndex[4] = {index, index, index, index};
+      difIndex[0][row] -= 2;
+      difIndex[1][row] -= 1;
+      difIndex[2][row] += 1;
+      difIndex[3][row] += 2;
+      if (difIndex[0][row] < startingIndex[row])
       {
-        difIndex[row][0][row] = index[row] + posoff;
-        ddrindex[row] = index[row] + posoff * 2;
+        difIndex[0][row] = startingIndex[row];
       }
-      if (index[row] > 1)
+      if (difIndex[3][row] > upperIndex[row])
       {
-        difIndex[row][1][row] = index[row] - 1;
-        ddlindex[row] = index[row] - 2;
+        difIndex[3][row] = upperIndex[row];
       }
 
-      OutputVectorType tempPix;
-
-      tempPix = m_DisplacementField->GetPixel(difIndex[row][1]);
-      const auto rpix = m_DisplacementField->TransformLocalVectorToPhysicalVector(tempPix);
-
-      tempPix = m_DisplacementField->GetPixel(difIndex[row][0]);
-      const auto lpix = m_DisplacementField->TransformLocalVectorToPhysicalVector(tempPix);
-
-      tempPix = m_DisplacementField->GetPixel(ddrindex);
-      const auto rrpix = m_DisplacementField->TransformLocalVectorToPhysicalVector(tempPix);
-
-      tempPix = m_DisplacementField->GetPixel(ddlindex);
-      const auto llpix = m_DisplacementField->TransformLocalVectorToPhysicalVector(tempPix);
-
+      OutputVectorType pixDisp[4];
+      for(unsigned int i = 0; i < 4; ++i)
+      {
+        const OutputVectorType tempPix = m_DisplacementField->GetPixel(difIndex[i]);
+        pixDisp[i] = m_DisplacementField->TransformLocalVectorToPhysicalVector(tempPix);
+      }
 
       // 4th order centered difference
-      OutputVectorType dPix = (lpix * 8.0 + llpix - rrpix - rpix * 8.0) * space / (12.0) * dPixSign;
+      OutputVectorType dPix = (pixDisp[0] - pixDisp[1] * 8.0 + pixDisp[2] * 8.0 - pixDisp[3]) / (12.0 * space * spacing[row]) * dPixSign;
 
-      // typename DisplacementFieldType::PixelType dPix=
-      //      ( lpix - rpix )*space/(2.0*h); //2nd order centered difference
+      dPix[row] += 1.0;
       for (unsigned int col = 0; col < NDimensions; ++col)
       {
-        TParametersValueType val = dPix[col] / spacing[row];
-        if (row == col)
-        {
-          val += 1.0;
-        }
-        jacobian(col, row) = val;
+        jacobian(col, row) = dPix[col];
         // Verify it's a real number
-        if (!itk::Math::isfinite(val))
+        if (!itk::Math::isfinite(dPix[col]))
         {
           isValidJacobianCalcLocat = false;
           break;