ENH: Add two new (simpler) computation modes

KitwareMedical · Apr 8, 2022 · acc89db · acc89db
1 parent 4ed2a5e
commit acc89db
Show file tree

Hide file tree

Showing 2 changed files with 84 additions and 35 deletions.
diff --git a/include/itkAttenuationImageFilter.h b/include/itkAttenuationImageFilter.h
@@ -50,11 +50,15 @@ namespace itk
  * region over which attenuations should be estimated.
  *
  * AttenuationImageFilter generates a scalar output image with
- * pixel intensities representing attenuation estimates between
- * the given pixel's location in the input image and either a pixel
- * that is a fixed distance away or a pixel at the end of the given
- * mask region, depending on filter settings. Pixels outside of
- * the mask after padding erosion will have a value of zero.
+ * pixel intensities representing attenuation estimates.
+ * Estimates are made in continuous segments of the mask region
+ * along the RF sampling direction.
+ *
+ * There are three modes of computation, see documentation for
+ * ComputationMode for detailed description.
+ * Pixels outside of the mask after padding erosion will have a value of zero.
+ * Pixels inside the mask for which the attenuation has not been computed
+ * will also have a value of zero.
  *
  * \sa MaskedImageToHistogramFilter
  * \sa Spectra1DImageFilter
@@ -113,6 +117,21 @@ class ITK_TEMPLATE_EXPORT AttenuationImageFilter : public ImageToImageFilter<TIn
    *  after padding is applied */
   itkGetConstMacro(OutputMaskImage, TMaskImage *);
 
+  /** Mode of computation (0=all pairs, 1=first and last, 2=fixed distance).
+   * 
+   * 0. (Default) Between all pixel pairs in a segment, producing a weighted
+   *    estimate for each pixel. More distant pixel pairs have higher weights.
+   *    This decreases influence of numerical instability.
+   * 1. Between first and last pixel in a segment.
+   * 2. Between first pixel in a segment, and
+   *    a pixel at a fixed distance away from it.
+   *    This distance is controlled by FixedEstimationDepth parameter.
+   *
+   * Modes 1 and 2 produce an estimate for the midpoint pixel only.
+   */
+  itkSetMacro(ComputationMode, unsigned int);
+  itkGetConstMacro(ComputationMode, unsigned int);
+
   /** Label value indicating which mask pixel values should be included in analysis.
    *  A value of zero indicates that any nonzero pixel should be included.
    */
@@ -239,6 +258,8 @@ class ITK_TEMPLATE_EXPORT AttenuationImageFilter : public ImageToImageFilter<TIn
   ThreadedIsIncluded(InputIndexType index) const;
 
 private:
+  unsigned int m_ComputationMode = 0;
+
   MaskPixelType m_LabelValue = 0;
 
   unsigned int m_Direction = 0;

diff --git a/include/itkAttenuationImageFilter.hxx b/include/itkAttenuationImageFilter.hxx
@@ -215,47 +215,75 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::ThreadedGenerateD
 
         if (start[m_Direction] < end[m_Direction]) // We need at least a pair of pixels to estimate attenuation
         {
-          // Estimate attenuation for each inclusion pixel
-          // by weighted average of pair-wise attenuations for all pairs
-          while (start[m_Direction] <= end[m_Direction])
+          if (m_ComputationMode == 0)
           {
-            for (IndexValueType k = start[m_Direction] + 1; k <= end[m_Direction]; ++k)
+            // Estimate attenuation for each inclusion pixel
+            // by weighted average of pair-wise attenuations for all pairs
+            while (start[m_Direction] <= end[m_Direction])
             {
-              unsigned pixelDistance = k - start[m_Direction];
-
-              InputIndexType target = start;
-              target[m_Direction] = k;
-              float estimatedAttenuation = ComputeAttenuation(target, start);
-              float weight = 1.0;                      // Weight for this pair's attenuation. 1 for large distances.
-              if (pixelDistance < distanceWeightsSize) // If pixels are close, weight is lower than 1.
+              for (IndexValueType k = start[m_Direction] + 1; k <= end[m_Direction]; ++k)
+              {
+                unsigned pixelDistance = k - start[m_Direction];
+
+                InputIndexType target = start;
+                target[m_Direction] = k;
+                float estimatedAttenuation = ComputeAttenuation(target, start);
+                float weight = 1.0;                      // Weight for this pair's attenuation. 1 for large distances.
+                if (pixelDistance < distanceWeightsSize) // If pixels are close, weight is lower than 1.
+                {
+                  weight = m_DistanceWeights[pixelDistance];
+                }
+
+                // Update this pixel
+                accumulatedWeight[start[m_Direction]] += weight;
+                output->SetPixel(start, estimatedAttenuation * weight + output->GetPixel(start));
+
+                // Update distant pair
+                accumulatedWeight[k] += weight;
+                output->SetPixel(target, estimatedAttenuation * weight + output->GetPixel(target));
+              } // for k
+
+              // Normalize output by accumulated weight
+              output->SetPixel(start, output->GetPixel(start) / accumulatedWeight[start[m_Direction]]);
+              accumulatedWeight[start[m_Direction]] = 0.0f; // reset for next next inclusion segment
+
+              // Possibly eliminate negative attenuations
+              if (!m_ConsiderNegativeAttenuations && output->GetPixel(start) < 0.0)
               {
-                weight = m_DistanceWeights[pixelDistance];
+                output->SetPixel(start, 0.0);
               }
 
-              // Update this pixel
-              accumulatedWeight[start[m_Direction]] += weight;
-              output->SetPixel(start, estimatedAttenuation * weight + output->GetPixel(start));
+              // Dynamically generate the output mask with values corresponding to input
+              m_OutputMaskImage->SetPixel(start, inputMaskImage->GetPixel(start));
 
-              // Update distant pair
-              accumulatedWeight[k] += weight;
-              output->SetPixel(target, estimatedAttenuation * weight + output->GetPixel(target));
-            } // for k
+              ++start[m_Direction];
+            } // while start<=end
+          }
+          else if (m_ComputationMode == 1 || m_ComputationMode == 2)
+          {
+            InputIndexType target = end;
+            IndexValueType fixedEnd = start[m_Direction] + m_FixedEstimationDepth;
+            if (m_ComputationMode == 2 && fixedEnd < end[m_Direction])
+            {
+              target[m_Direction] = fixedEnd;
+            }
+            float estimatedAttenuation = ComputeAttenuation(target, start);
 
-            // Normalize output by accumulated weight
-            output->SetPixel(start, output->GetPixel(start) / accumulatedWeight[start[m_Direction]]);
-            accumulatedWeight[start[m_Direction]] = 0.0f; // reset for next next inclusion segment
+            // Compute mid-point index
+            target[m_Direction] = (start[m_Direction] + target[m_Direction]) / 2;
+            output->SetPixel(target, estimatedAttenuation);
+            m_OutputMaskImage->SetPixel(target, inputMaskImage->GetPixel(target));
 
             // Possibly eliminate negative attenuations
-            if (!m_ConsiderNegativeAttenuations && output->GetPixel(start) < 0.0)
+            if (!m_ConsiderNegativeAttenuations && estimatedAttenuation < 0.0)
             {
-              output->SetPixel(start, 0.0);
+              output->SetPixel(target, 0.0);
             }
-
-            // Dynamically generate the output mask with values corresponding to input
-            m_OutputMaskImage->SetPixel(start, inputMaskImage->GetPixel(start));
-
-            ++start[m_Direction];
-          } // while start<=end
+          }
+          else
+          {
+            itkExceptionMacro(<< "Invalid computation mode: " << m_ComputationMode);
+          }
         } // if start<end
       } // else !inside