ENH: Allow to run AttenuationImageFilter without a mask

This computes attenuation estimates for all pixels in the image.

include/itkAttenuationImageFilter.h

@@ -285,6 +285,8 @@ class ITK_TEMPLATE_EXPORT AttenuationImageFilter : public ImageToImageFilter<TIn
   /** Cache mask image reference before threaded execution to reduce calls to GetMaskImage() */
   mutable const MaskImageType * m_ThreadedInputMaskImage;
 
+  unsigned int m_LastScanlineIndex = 0;
+
   /** Output mask image may be eroded via m_PadUpperBounds and m_PadLowerBounds
    *  along scan line direction */
   MaskImagePointer m_OutputMaskImage = MaskImageType::New();

include/itkAttenuationImageFilter.hxx

@@ -95,14 +95,7 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::VerifyPreconditio
 {
   Superclass::VerifyPreconditions();
 
-  if (this->GetInputMaskImage() == nullptr)
-  {
-    itkExceptionMacro("Filter requires a mask image for inclusion estimates!");
-  }
-  else
-  {
-    m_ThreadedInputMaskImage = this->GetInputMaskImage();
-  }
+  m_ThreadedInputMaskImage = this->GetInputMaskImage();
 
   if (this->GetDirection() >= ImageDimension)
   {
@@ -126,11 +119,26 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::BeforeThreadedGen
   this->GetOutput()->FillBuffer(0.0f);
 
   // Initialize output mask image
-  const MaskImageType * inputMaskImage = this->GetInputMaskImage();
-  m_OutputMaskImage->CopyInformation(inputMaskImage);
-  m_OutputMaskImage->SetRegions(inputMaskImage->GetLargestPossibleRegion());
-  m_OutputMaskImage->Allocate();
-  m_OutputMaskImage->FillBuffer(0.0f);
+  const InputImageType * input = this->GetInput();
+  const MaskImageType *  inputMaskImage = this->GetInputMaskImage();
+
+  auto region = input->GetLargestPossibleRegion();
+  if (inputMaskImage != nullptr)
+  {
+    m_OutputMaskImage->CopyInformation(inputMaskImage);
+    region = inputMaskImage->GetLargestPossibleRegion();
+    m_OutputMaskImage->SetRegions(region);
+    m_OutputMaskImage->Allocate();
+    m_OutputMaskImage->FillBuffer(0);
+  }
+  else
+  {
+    m_OutputMaskImage->CopyInformation(input);
+    m_OutputMaskImage->SetRegions(region);
+    m_OutputMaskImage->Allocate();
+    m_OutputMaskImage->FillBuffer(1);
+  }
+  m_LastScanlineIndex = region.GetIndex(m_Direction) + region.GetSize(m_Direction) - 1;
 
   // Initialize distance weights
   unsigned fourSigma = std::max(m_FixedEstimationDepth, 16u); // An arbitrary default.
@@ -143,7 +151,7 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::BeforeThreadedGen
 
   // Initialize iVars used in ComputeAttenuation()
   float nyquistFrequency = m_SamplingFrequencyMHz / 2;
-  float numComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
+  float numComponents = input->GetNumberOfComponentsPerPixel();
   m_FrequencyDelta = nyquistFrequency / numComponents;
   m_StartComponent = m_FrequencyBandStartMHz / m_FrequencyDelta;
   m_EndComponent = m_FrequencyBandEndMHz / m_FrequencyDelta;
@@ -152,7 +160,7 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::BeforeThreadedGen
     m_EndComponent = numComponents - 1; // Use all components
   }
   m_ConsideredComponents = m_EndComponent - m_StartComponent + 1;
-  m_ScanStepMM = this->GetInput()->GetSpacing()[m_Direction];
+  m_ScanStepMM = input->GetSpacing()[m_Direction];
 }
 
 template <typename TInputImage, typename TOutputImage, typename TMaskImage>
@@ -193,7 +201,7 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::ThreadedGenerateD
       // Advance until an inclusion is found
       InputIndexType index = it.GetIndex();
       bool           inside = ThreadedIsIncluded(index);
-      if (inside)
+      if (inside && index[m_Direction] < m_LastScanlineIndex)
       {
         if (inclusionLength == 0)
         {
@@ -205,9 +213,11 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::ThreadedGenerateD
       {
         inclusionLength = 0; // Prepare for the next one
 
-        // Stay at last pixel in the inclusion
-        end = it.GetIndex();
-        end[m_Direction] -= 1;
+        end = index;
+        if (index[m_Direction] < m_LastScanlineIndex)
+        {
+          end[m_Direction] -= 1; // Stay at last pixel in the inclusion
+        }
 
         // Adjust for pixel padding
         start[m_Direction] += m_PadLowerBounds;
@@ -248,7 +258,7 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::ThreadedGenerateD
               accumulatedWeight[start[m_Direction]] = 0.0f; // reset for next next inclusion segment
 
               // Only set mask for valid estimates
-              if (m_ConsiderNegativeAttenuations || output->GetPixel(start) >= 0.0)
+              if (inputMaskImage != nullptr && (m_ConsiderNegativeAttenuations || output->GetPixel(start) >= 0.0))
               {
                 // Dynamically generate the output mask with values corresponding to input
                 m_OutputMaskImage->SetPixel(start, inputMaskImage->GetPixel(start));
@@ -272,7 +282,7 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::ThreadedGenerateD
             output->SetPixel(target, estimatedAttenuation);
 
             // Only set mask for valid estimates
-            if (m_ConsiderNegativeAttenuations || estimatedAttenuation >= 0.0)
+            if (inputMaskImage != nullptr && (m_ConsiderNegativeAttenuations || estimatedAttenuation >= 0.0))
             {
               m_OutputMaskImage->SetPixel(target, inputMaskImage->GetPixel(target));
             }
@@ -282,7 +292,7 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::ThreadedGenerateD
             itkExceptionMacro(<< "Invalid computation mode: " << m_ComputationMode);
           }
         } // if start<end
-      } // else !inside
+      }   // else !inside
 
       ++it;
     }
@@ -364,6 +374,10 @@ template <typename TInputImage, typename TOutputImage, typename TMaskImage>
 bool
 AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::ThreadedIsIncluded(InputIndexType index) const
 {
+  if (m_ThreadedInputMaskImage == nullptr)
+  {
+    return true;
+  }
   auto maskValue = m_ThreadedInputMaskImage->GetPixel(index);
   return (m_LabelValue == 0 && maskValue > 0) || (m_LabelValue > 0 && maskValue == m_LabelValue);
 }

test/Baseline/itkAttenuationImageFilterWholeImageTestMask.mha.sha512

@@ -0,0 +1 @@
+2b943f95e921b98f15170ca74b19ccd849e6daa3c8735a4ef09fa9c82e7429c13f9d6435fb5a0759fcc11d68898ac7b19181d76eb35fd54a77114c831fcbe171

test/Baseline/itkAttenuationImageFilterWholeImageTestOutput.mha.sha512

@@ -0,0 +1 @@
+089ccbf2e02e6e30a3165f07508994c4bdd30d6ec7ee2092eae2dc6722264ab02bb8ffb4c338b2d33a0da1091ed22b0556f36e6150b7f5b020a924ddf6cdac1d

test/CMakeLists.txt

@@ -148,6 +148,25 @@ itk_add_test(NAME itkAttenuationImageFilterPairModeTest
     1   # computation mode
     0.0693445 # median attenuation
   )
+itk_add_test(NAME itkAttenuationImageFilterWholeImageTest
+  COMMAND UltrasoundTestDriver
+  --compare
+    DATA{Baseline/itkAttenuationImageFilterWholeImageTestOutput.mha}
+    ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterWholeImageTestOutput.mha
+  --compare
+    DATA{Baseline/itkAttenuationImageFilterWholeImageTestMask.mha}
+    ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterWholeImageTestMask.mha
+  itkAttenuationImageFilterTest
+    DATA{Input/UnfusedRF-a0-spectra-cropped.mhd,UnfusedRF-a0-spectra-cropped.raw}
+    nul
+    ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterWholeImageTestOutput.mha
+    ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterWholeImageTestMask.mha
+    100 # work units
+    0.0 # Distance between pixels sampled for attenuation estimation
+    1   # consider negative attenuations
+    0   # computation mode
+    0.458781 # median attenuation
+  )
 itk_add_test(NAME itkBModeImageFilterTestTiming
   COMMAND UltrasoundTestDriver
   --compare

test/itkAttenuationImageFilterTest.cxx

@@ -48,7 +48,12 @@ itkAttenuationImageFilterTest(int argc, char * argv[])
   using MaskImageType = itk::Image<LabelType, Dimension>;
 
   SpectraImageType::Pointer inputImage = itk::ReadImage<SpectraImageType>(std::string(argv[1]));
-  MaskImageType::Pointer    maskImage = itk::ReadImage<MaskImageType>(std::string(argv[2]));
+  MaskImageType::Pointer    maskImage;
+  std::string               maskFilename = argv[2];
+  if (maskFilename != "nul")
+  {
+    maskImage = itk::ReadImage<MaskImageType>(maskFilename);
+  }
 
   const std::string outputImagePath = argv[3];
   const std::string outputMaskImagePath = (argc > 4 ? argv[4] : "");
@@ -70,9 +75,6 @@ itkAttenuationImageFilterTest(int argc, char * argv[])
   attenuationFilter->SetInput(inputImage);
   ITK_TEST_SET_GET_VALUE(inputImage, attenuationFilter->GetInput());
 
-  // Verify filter does not run without mask
-  ITK_TRY_EXPECT_EXCEPTION(attenuationFilter->Update());
-
   attenuationFilter->SetInputMaskImage(maskImage);
   ITK_TEST_SET_GET_VALUE(maskImage, attenuationFilter->GetInputMaskImage());
 
@@ -133,16 +135,20 @@ itkAttenuationImageFilterTest(int argc, char * argv[])
     itk::WriteImage(attenuationFilter->GetOutputMaskImage(), outputMaskImagePath, true);
   }
 
-  // Discover mask's inside value
-  using MaxType = itk::MinimumMaximumImageCalculator<MaskImageType>;
-  auto maxCalculator = MaxType::New();
-  maxCalculator->SetImage(maskImage);
-  maxCalculator->ComputeMaximum();
-  LabelType maxMaskValue = maxCalculator->GetMaximum();
-  if (maxMaskValue == 0)
+  LabelType maxMaskValue = 1;
+  if (maskFilename != "nul")
   {
-    std::cerr << "The mask is empty!" << std::endl;
-    return EXIT_FAILURE;
+    // Discover mask's inside value
+    using MaxType = itk::MinimumMaximumImageCalculator<MaskImageType>;
+    auto maxCalculator = MaxType::New();
+    maxCalculator->SetImage(maskImage);
+    maxCalculator->ComputeMaximum();
+    LabelType maxMaskValue = maxCalculator->GetMaximum();
+    if (maxMaskValue == 0)
+    {
+      std::cerr << "The mask is empty!" << std::endl;
+      return EXIT_FAILURE;
+    }
   }
 
   // Now boil it down to a single attenuation value