ENH: Generate output mask for attenuation statistics

include/itkAttenuationImageFilter.h

@@ -103,8 +103,13 @@ class ITK_TEMPLATE_EXPORT AttenuationImageFilter : public ImageToImageFilter<TIn
   using OutputRegionType = typename TOutputImage::RegionType;
   using OutputPixelType = typename TOutputImage::PixelType;
 
-  itkSetInputMacro(MaskImage, TMaskImage);
-  itkGetInputMacro(MaskImage, TMaskImage);
+  /** Input mask image represents input region for analysis */
+  itkSetInputMacro(InputMaskImage, TMaskImage);
+  itkGetInputMacro(InputMaskImage, TMaskImage);
+
+  /** Output mask image represents output region for analysis
+   *  after padding is applied */
+  itkGetConstMacro(OutputMaskImage, TMaskImage *);
 
   /** Label value indicating which mask pixel values should be included in analysis.
    *  A value of zero indicates that any nonzero pixel should be included.
@@ -248,7 +253,11 @@ class ITK_TEMPLATE_EXPORT AttenuationImageFilter : public ImageToImageFilter<TIn
   ImageRegionSplitterDirection::Pointer m_RegionSplitter = ImageRegionSplitterDirection::New();
 
   /** Cache mask image reference before threaded execution to reduce calls to GetMaskImage() */
-  const MaskImageType * m_ThreadedMaskImage;
+  const MaskImageType * m_ThreadedInputMaskImage;
+
+  /** Output mask image may be eroded via m_PadUpperBounds and m_PadLowerBounds
+   *  along scan line direction */
+  MaskImagePointer m_OutputMaskImage = MaskImageType::New();
 };
 } // end namespace itk
 

include/itkAttenuationImageFilter.hxx

@@ -91,13 +91,13 @@ void
 AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::BeforeThreadedGenerateData()
 {
   // Verify inputs
-  if (this->GetMaskImage() == nullptr)
+  if (this->GetInputMaskImage() == nullptr)
   {
     itkExceptionMacro("Filter requires a mask image for inclusion estimates!");
   }
   else
   {
-    m_ThreadedMaskImage = this->GetMaskImage();
+    m_ThreadedInputMaskImage = this->GetInputMaskImage();
   }
 
   if (this->GetDirection() >= ImageDimension)
@@ -115,6 +115,13 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::BeforeThreadedGen
   // Initialize metric image
   this->GetOutput()->Allocate();
   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);
 }
 
 template <typename TInputImage, typename TOutputImage, typename TMaskImage>
@@ -127,9 +134,9 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::DynamicThreadedGe
     return;
   }
 
-  const MaskImageType *  maskImage = this->GetMaskImage();
   const InputImageType * input = this->GetInput();
   OutputImageType *      output = this->GetOutput();
+  const MaskImageType *  inputMaskImage = this->GetInputMaskImage();
 
   ImageLinearConstIteratorWithIndex<TInputImage> it(input, regionForThread);
   it.SetDirection(m_Direction);
@@ -186,6 +193,9 @@ AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::DynamicThreadedGe
             output->SetPixel(start, estimatedAttenuation);
           }
 
+          // Dynamically generate the output mask with values corresponding to input
+          m_OutputMaskImage->SetPixel(start, inputMaskImage->GetPixel(start));
+
           ++start[m_Direction];
         }
 
@@ -289,7 +299,7 @@ template <typename TInputImage, typename TOutputImage, typename TMaskImage>
 bool
 AttenuationImageFilter<TInputImage, TOutputImage, TMaskImage>::ThreadedIsIncluded(InputIndexType index) const
 {
-  auto maskValue = m_ThreadedMaskImage->GetPixel(index);
+  auto maskValue = m_ThreadedInputMaskImage->GetPixel(index);
   return (m_LabelValue == 0 && maskValue > 0) || (m_LabelValue > 0 && maskValue == m_LabelValue);
 }
 

test/Baseline/itkAttenuationImageFilterTestMaskOutput.mha.sha512

@@ -0,0 +1 @@
+5df3280bf252942f3308e5321b2fdf15cdcd6a61793f089f3eca951edf1203f0ba57aed07e90d19c349171f8790d55dcc44bcf669361a2a8686a7c8e3f5a8d87

test/CMakeLists.txt

@@ -56,31 +56,45 @@ itk_add_test(NAME itkAttenuationImageFilterTest
   --compare
     DATA{Baseline/itkAttenuationImageFilterTestOutputBaseline.mha}
     ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterTestOutput.mha
+  --compare
+    DATA{Baseline/itkAttenuationImageFilterTestMaskOutput.mha}
+    ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterTestMaskOutput.mha
   itkAttenuationImageFilterTest
     DATA{Input/UnfusedRF-a0-spectra-cropped.mhd,UnfusedRF-a0-spectra-cropped.raw}
     DATA{Input/UnfusedRF-a0-spectra-label_1-cropped.mhd,UnfusedRF-a0-spectra-label_1-cropped.raw}
     ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterTestOutput.mha
+    ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterTestMaskOutput.mha
+    1 # work units
+    0.0 # Distance between pixels sampled for attenuation estimation; 0.0 -> inclusion length
   )
 itk_add_test(NAME itkAttenuationImageFilterMultithreadedTest
   COMMAND UltrasoundTestDriver
   --compare
     DATA{Baseline/itkAttenuationImageFilterTestOutputBaseline.mha}
     ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterMultithreadedTestOutput.mha
+  --compare
+    DATA{Baseline/itkAttenuationImageFilterTestMaskOutput.mha}
+    ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterTestMaskOutput.mha
   itkAttenuationImageFilterTest
     DATA{Input/UnfusedRF-a0-spectra-cropped.mhd,UnfusedRF-a0-spectra-cropped.raw}
     DATA{Input/UnfusedRF-a0-spectra-label_1-cropped.mhd,UnfusedRF-a0-spectra-label_1-cropped.raw}
     ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterMultithreadedTestOutput.mha
+    ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterTestMaskOutput.mha
     5 # work units
   )
 itk_add_test(NAME itkAttenuationImageFilterFixedEstimationDistanceTest
   COMMAND UltrasoundTestDriver
   --compare
     DATA{Baseline/itkAttenuationImageFilterTestFixedEstimationDepthBaseline.mha}
     ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterTestFixedEstimationDepthOutput.mha
+  --compare
+    DATA{Baseline/itkAttenuationImageFilterTestMaskOutput.mha}
+    ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterTestMaskOutput.mha
   itkAttenuationImageFilterTest
     DATA{Input/UnfusedRF-a0-spectra-cropped.mhd,UnfusedRF-a0-spectra-cropped.raw}
     DATA{Input/UnfusedRF-a0-spectra-label_1-cropped.mhd,UnfusedRF-a0-spectra-label_1-cropped.raw}
     ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterTestFixedEstimationDepthOutput.mha
+    ${ITK_TEST_OUTPUT_DIR}/itkAttenuationImageFilterTestMaskOutput.mha
     5   # work units
     3.0 # Distance between pixels sampled for attenuation estimation
   )

test/itkAttenuationImageFilterTest.cxx

@@ -30,7 +30,7 @@ itkAttenuationImageFilterTest(int argc, char * argv[])
   if (argc < 3)
   {
     std::cerr << "Usage: " << argv[0];
-    std::cerr << " spectraImage maskImage outputImage <numWorkUnits> <fixedEstimationDepthMM>";
+    std::cerr << " spectraImage maskImage outputImage <outputMaskImage> <numWorkUnits> <fixedEstimationDepthMM>";
     std::cerr << std::endl;
     return EXIT_FAILURE;
   }
@@ -46,6 +46,12 @@ itkAttenuationImageFilterTest(int argc, char * argv[])
   SpectraImageType::Pointer inputImage = itk::ReadImage<SpectraImageType>(std::string(argv[1]));
   MaskImageType::Pointer    maskImage = itk::ReadImage<MaskImageType>(std::string(argv[2]));
 
+  const std::string outputImagePath = argv[3];
+  const std::string outputMaskImagePath = (argc > 4 ? argv[4] : "");
+
+  unsigned int numWorkUnits = (argc > 5 ? std::atoi(argv[5]) : 1);
+  float        fixedEstimationDepthMM = (argc > 6 ? std::atof(argv[6]) : 0.0);
+
   // Initialize the filter
   using AttenuationFilterType = itk::AttenuationImageFilter<SpectraImageType, OutputImageType, MaskImageType>;
   AttenuationFilterType::Pointer attenuationFilter = AttenuationFilterType::New();
@@ -61,8 +67,8 @@ itkAttenuationImageFilterTest(int argc, char * argv[])
   // Verify filter does not run without mask
   ITK_TRY_EXPECT_EXCEPTION(attenuationFilter->Update());
 
-  attenuationFilter->SetMaskImage(maskImage);
-  ITK_TEST_SET_GET_VALUE(maskImage, attenuationFilter->GetMaskImage());
+  attenuationFilter->SetInputMaskImage(maskImage);
+  ITK_TEST_SET_GET_VALUE(maskImage, attenuationFilter->GetInputMaskImage());
 
   // Bad scanline direction produces an error
   attenuationFilter->SetScanDirection(Dimension);
@@ -82,7 +88,6 @@ itkAttenuationImageFilterTest(int argc, char * argv[])
   ITK_TEST_SET_GET_VALUE(3, attenuationFilter->GetFixedEstimationDepth());
 
   // Verify spatial distance is estimated to the nearest pixel center
-  float fixedEstimationDepthMM = (argc > 5 ? std::atof(argv[5]) : 0.0);
   attenuationFilter->SetFixedEstimationDepthMM(fixedEstimationDepthMM);
   ITK_TEST_EXPECT_TRUE(fabs(attenuationFilter->GetFixedEstimationDepthMM() - fixedEstimationDepthMM) <
                        (inputImage->GetSpacing()[0] / 2));
@@ -105,7 +110,6 @@ itkAttenuationImageFilterTest(int argc, char * argv[])
   attenuationFilter->SetPadUpperBoundsMM(0.0);
   ITK_TEST_EXPECT_TRUE(fabs(attenuationFilter->GetPadUpperBoundsMM() - 0.0) < (inputImage->GetSpacing()[0] / 2));
 
-  unsigned int numWorkUnits = (argc > 4 ? std::atoi(argv[4]) : 1);
   attenuationFilter->SetNumberOfWorkUnits(numWorkUnits);
 
   ITK_EXERCISE_BASIC_OBJECT_METHODS(attenuationFilter, AttenuationImageFilter, ImageToImageFilter);
@@ -114,7 +118,13 @@ itkAttenuationImageFilterTest(int argc, char * argv[])
   ITK_TRY_EXPECT_NO_EXCEPTION(attenuationFilter->Update());
 
   // Verify output
-  itk::WriteImage(attenuationFilter->GetOutput(), argv[3]);
+  itk::WriteImage(attenuationFilter->GetOutput(), outputImagePath);
+
+  // Verify mask output after erosion
+  if (outputMaskImagePath != "")
+  {
+    itk::WriteImage(attenuationFilter->GetOutputMaskImage(), outputMaskImagePath);
+  }
 
   return EXIT_SUCCESS;
 }