From 46f2c995a5a9c93c531a08626a1e3a2ab8c0f1ee Mon Sep 17 00:00:00 2001 From: Tom Birdsong Date: Wed, 2 Mar 2022 08:45:50 -0500 Subject: [PATCH] ENH: Generate output mask for attenuation statistics --- include/itkAttenuationImageFilter.h | 15 ++++++++++--- include/itkAttenuationImageFilter.hxx | 18 +++++++++++---- ...uationImageFilterTestMaskOutput.mha.sha512 | 1 + test/CMakeLists.txt | 14 ++++++++++++ test/itkAttenuationImageFilterTest.cxx | 22 ++++++++++++++----- 5 files changed, 57 insertions(+), 13 deletions(-) create mode 100644 test/Baseline/itkAttenuationImageFilterTestMaskOutput.mha.sha512 diff --git a/include/itkAttenuationImageFilter.h b/include/itkAttenuationImageFilter.h index 35e49385..ab04ee79 100644 --- a/include/itkAttenuationImageFilter.h +++ b/include/itkAttenuationImageFilter.h @@ -103,8 +103,13 @@ class ITK_TEMPLATE_EXPORT AttenuationImageFilter : public ImageToImageFilter::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::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 @@ -127,9 +134,9 @@ AttenuationImageFilter::DynamicThreadedGe return; } - const MaskImageType * maskImage = this->GetMaskImage(); const InputImageType * input = this->GetInput(); OutputImageType * output = this->GetOutput(); + const MaskImageType * inputMaskImage = this->GetInputMaskImage(); ImageLinearConstIteratorWithIndex it(input, regionForThread); it.SetDirection(m_Direction); @@ -186,6 +193,9 @@ AttenuationImageFilter::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 bool AttenuationImageFilter::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); } diff --git a/test/Baseline/itkAttenuationImageFilterTestMaskOutput.mha.sha512 b/test/Baseline/itkAttenuationImageFilterTestMaskOutput.mha.sha512 new file mode 100644 index 00000000..daba6c58 --- /dev/null +++ b/test/Baseline/itkAttenuationImageFilterTestMaskOutput.mha.sha512 @@ -0,0 +1 @@ +5df3280bf252942f3308e5321b2fdf15cdcd6a61793f089f3eca951edf1203f0ba57aed07e90d19c349171f8790d55dcc44bcf669361a2a8686a7c8e3f5a8d87 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c4fa9a86..6d31aa88 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -56,20 +56,30 @@ 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 @@ -77,10 +87,14 @@ itk_add_test(NAME itkAttenuationImageFilterFixedEstimationDistanceTest --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 ) diff --git a/test/itkAttenuationImageFilterTest.cxx b/test/itkAttenuationImageFilterTest.cxx index f7ed2462..fb8154fa 100644 --- a/test/itkAttenuationImageFilterTest.cxx +++ b/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 "; + std::cerr << " spectraImage maskImage outputImage "; std::cerr << std::endl; return EXIT_FAILURE; } @@ -46,6 +46,12 @@ itkAttenuationImageFilterTest(int argc, char * argv[]) SpectraImageType::Pointer inputImage = itk::ReadImage(std::string(argv[1])); MaskImageType::Pointer maskImage = itk::ReadImage(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; 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; }