From 61dfbf3fb8327c6fc6c3d8a14e32e88d382dab9d Mon Sep 17 00:00:00 2001 From: Olivier Commowick Date: Thu, 14 Mar 2019 12:15:16 +0100 Subject: [PATCH 1/3] Add lesion evolution detection tool --- .../validation_tools/CMakeLists.txt | 1 + .../lesion_evolution_detection/CMakeLists.txt | 36 +++ .../animaLesionEvolutionDetection.cxx | 276 ++++++++++++++++++ 3 files changed, 313 insertions(+) create mode 100644 Anima/segmentation/validation_tools/lesion_evolution_detection/CMakeLists.txt create mode 100644 Anima/segmentation/validation_tools/lesion_evolution_detection/animaLesionEvolutionDetection.cxx diff --git a/Anima/segmentation/validation_tools/CMakeLists.txt b/Anima/segmentation/validation_tools/CMakeLists.txt index c28196de5..849cc36ce 100644 --- a/Anima/segmentation/validation_tools/CMakeLists.txt +++ b/Anima/segmentation/validation_tools/CMakeLists.txt @@ -1,3 +1,4 @@ add_subdirectory(detected_components) add_subdirectory(dice_measure) +add_subdirectory(lesion_evolution_detection) add_subdirectory(segmentation_performance_analyzer) diff --git a/Anima/segmentation/validation_tools/lesion_evolution_detection/CMakeLists.txt b/Anima/segmentation/validation_tools/lesion_evolution_detection/CMakeLists.txt new file mode 100644 index 000000000..b3a6004c1 --- /dev/null +++ b/Anima/segmentation/validation_tools/lesion_evolution_detection/CMakeLists.txt @@ -0,0 +1,36 @@ +if(BUILD_TOOLS) + +project(animaLesionEvolutionDetection) + +## ############################################################################# +## List Sources +## ############################################################################# + +list_source_files(${PROJECT_NAME} + ${CMAKE_CURRENT_SOURCE_DIR} + ) + +## ############################################################################# +## add executable +## ############################################################################# + +add_executable(${PROJECT_NAME} + ${${PROJECT_NAME}_CFILES} + ) + + +## ############################################################################# +## Link +## ############################################################################# + +target_link_libraries(${PROJECT_NAME} + ${ITKIO_LIBRARIES} + ) + +## ############################################################################# +## install +## ############################################################################# + +set_exe_install_rules(${PROJECT_NAME}) + +endif() diff --git a/Anima/segmentation/validation_tools/lesion_evolution_detection/animaLesionEvolutionDetection.cxx b/Anima/segmentation/validation_tools/lesion_evolution_detection/animaLesionEvolutionDetection.cxx new file mode 100644 index 000000000..078268b17 --- /dev/null +++ b/Anima/segmentation/validation_tools/lesion_evolution_detection/animaLesionEvolutionDetection.cxx @@ -0,0 +1,276 @@ +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +int main(int argc, char **argv) +{ + TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + + TCLAP::ValueArg refArg("r", "ref", "Reference binary segmentation", true, "", "reference segmentation", cmd); + TCLAP::ValueArg testArg("t", "test", "Test binary segmentation", true, "", "test segmentation", cmd); + TCLAP::ValueArg outArg("o", "out", "output image with labeled lesions", true, "", "output labeled lesions", cmd); + + TCLAP::ValueArg alphaArg("a", "alpha", "Alpha threshold (in between 0 and 1, default: 0)", false, 0, "alpha threshold", cmd); + TCLAP::ValueArg gammaArg("g", "gamma", "Gamma threshold (*100 %, default: 0.5)", false, 0.5, "gamma threshold", cmd); + TCLAP::ValueArg betaArg("b", "beta", "Beta threshold (in between 0 and 1, default: 0.05)", false, 0.05, "beta threshold", cmd); + + TCLAP::ValueArg minVolumeArg("m", "min-vol", "Minimal volume for the component to be considered (default: 3 mm3)", false, 3, "minimal volume", cmd); + TCLAP::ValueArg nbpArg("T","numberofthreads","Number of threads to run on (default : all cores)",false,itk::MultiThreader::GetGlobalDefaultNumberOfThreads(),"number of threads",cmd); + + TCLAP::SwitchArg fullConnectArg("F","full-connect","Use 26-connectivity instead of 6-connectivity",cmd,false); + + try + { + cmd.parse(argc, argv); + } + catch (TCLAP::ArgException& e) + { + std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; + return EXIT_FAILURE; + } + + typedef itk::Image ImageType; + typedef itk::ImageRegionIterator ImageIteratorType; + + ImageType::Pointer refSegmentation = anima::readImage (refArg.getValue()); + ImageType::Pointer testSegmentation = anima::readImage (testArg.getValue()); + + typedef itk::ConnectedComponentImageFilter CCFilterType; + typedef itk::RelabelComponentImageFilter RelabelComponentFilterType; + + CCFilterType::Pointer refCCFilter = CCFilterType::New(); + refCCFilter->SetInput(refSegmentation); + refCCFilter->SetFullyConnected(fullConnectArg.isSet()); + refCCFilter->SetNumberOfThreads(nbpArg.getValue()); + refCCFilter->Update(); + + ImageType::SpacingType spacing = refSegmentation->GetSpacing(); + ImageType::SpacingValueType spacingTot = spacing[0]; + for (unsigned int i = 1;i < 3;++i) + spacingTot *= spacing[i]; + + // Compute minsize in voxels + unsigned int minSizeInVoxel = static_cast (std::ceil(minVolumeArg.getValue() / spacingTot)); + + // Remove too small reference objects + RelabelComponentFilterType::Pointer relabelRefFilter = RelabelComponentFilterType::New(); + relabelRefFilter->SetInput(refCCFilter->GetOutput()); + relabelRefFilter->SetMinimumObjectSize(minSizeInVoxel); + relabelRefFilter->SetNumberOfThreads(nbpArg.getValue()); + relabelRefFilter->Update(); + + // Reference segmentation is now labeled per connected objects + refSegmentation = relabelRefFilter->GetOutput(); + refSegmentation->DisconnectPipeline(); + + CCFilterType::Pointer testCCFilter = CCFilterType::New(); + testCCFilter->SetInput(testSegmentation); + testCCFilter->SetFullyConnected(fullConnectArg.isSet()); + testCCFilter->SetNumberOfThreads(nbpArg.getValue()); + testCCFilter->Update(); + + // Remove too small test objects + RelabelComponentFilterType::Pointer relabelTestFilter = RelabelComponentFilterType::New(); + relabelTestFilter->SetInput(testCCFilter->GetOutput()); + relabelTestFilter->SetMinimumObjectSize(minSizeInVoxel); + relabelTestFilter->SetNumberOfThreads(nbpArg.getValue()); + relabelTestFilter->Update(); + + // Test segmentation is now labeled per connected objects + testSegmentation = relabelTestFilter->GetOutput(); + testSegmentation->DisconnectPipeline(); + + ImageIteratorType testItr(testSegmentation, testSegmentation->GetLargestPossibleRegion()); + + unsigned short maxTestLabel = 0; + while (!testItr.IsAtEnd()) + { + if (testItr.Get() > maxTestLabel) + maxTestLabel = testItr.Get(); + + ++testItr; + } + + ++maxTestLabel; + + ImageIteratorType refItr(refSegmentation, refSegmentation->GetLargestPossibleRegion()); + std::vector labelsOverlap(maxTestLabel,0); + std::vector labelsNonOverlapping(maxTestLabel,0); + std::vector labelsSizes(maxTestLabel,0); + + ImageType::Pointer subImage = ImageType::New(); + subImage->Initialize(); + subImage->SetRegions(testSegmentation->GetLargestPossibleRegion()); + subImage->SetSpacing (testSegmentation->GetSpacing()); + subImage->SetOrigin (testSegmentation->GetOrigin()); + subImage->SetDirection (testSegmentation->GetDirection()); + subImage->Allocate(); + ImageIteratorType subItr(subImage,testSegmentation->GetLargestPossibleRegion()); + + testItr.GoToBegin(); + while (!refItr.IsAtEnd()) + { + if (refItr.Get() != 0) + labelsOverlap[testItr.Get()]++; + else + labelsNonOverlapping[testItr.Get()]++; + + labelsSizes[testItr.Get()]++; + + ++refItr; + ++testItr; + } + + std::cout << "Processing " << maxTestLabel << " lesions in second timepoint..." << std::endl; + for (unsigned int i = 1;i < maxTestLabel;++i) + { + testItr.GoToBegin(); + double ratioNonOverlapOverlap = static_cast (labelsNonOverlapping[i]) / labelsOverlap[i]; + // Shrinking or not enough change -> label = maxTestLabel + 1 + if ((labelsNonOverlapping[i] <= minSizeInVoxel) || (ratioNonOverlapOverlap <= betaArg.getValue())) + { + while (!testItr.IsAtEnd()) + { + if (testItr.Get() == i) + testItr.Set(maxTestLabel + 1); + + ++testItr; + } + + continue; + } + + // New lesion -> put it all as new -> label = maxTestLabel + 2 + double ratioOverlapSizes = static_cast (labelsOverlap[i]) / labelsSizes[i]; + if (ratioOverlapSizes <= alphaArg.getValue()) + { + while (!testItr.IsAtEnd()) + { + if (testItr.Get() == i) + testItr.Set(maxTestLabel + 2); + + ++testItr; + } + + continue; + } + + // Test for growing lesion too large + if (ratioNonOverlapOverlap > gammaArg.getValue()) + { + // New lesion -> label = maxTestLabel + 2 + while (!testItr.IsAtEnd()) + { + if (testItr.Get() == i) + testItr.Set(maxTestLabel + 2); + + ++testItr; + } + + continue; + } + + // We are looking at a growing lesion -> label = maxTestLabel + 3 + subImage->FillBuffer(0); + refItr.GoToBegin(); + subItr.GoToBegin(); + while (!testItr.IsAtEnd()) + { + if ((testItr.Get() == i)&&(refItr.Get() == 0)) + subItr.Set(1); + + ++testItr; + ++refItr; + ++subItr; + } + + CCFilterType::Pointer tmpCCFilter = CCFilterType::New(); + tmpCCFilter->SetInput(subImage); + tmpCCFilter->SetFullyConnected(fullConnectArg.isSet()); + tmpCCFilter->SetNumberOfThreads(nbpArg.getValue()); + tmpCCFilter->Update(); + + // Remove too small test objects + RelabelComponentFilterType::Pointer relabelTmpCCFilter = RelabelComponentFilterType::New(); + relabelTmpCCFilter->SetInput(tmpCCFilter->GetOutput()); + relabelTmpCCFilter->SetMinimumObjectSize(0); + relabelTmpCCFilter->SetNumberOfThreads(nbpArg.getValue()); + relabelTmpCCFilter->Update(); + + ImageIteratorType subCCItr(relabelTmpCCFilter->GetOutput(),testSegmentation->GetLargestPossibleRegion()); + std::vector numVoxelsCCSub(1); + while (!subCCItr.IsAtEnd()) + { + unsigned int currentSize = numVoxelsCCSub.size(); + unsigned int value = subCCItr.Get(); + + if (value >= currentSize) + { + numVoxelsCCSub.resize(value + 1); + for (unsigned int j = currentSize;j <= value;++j) + numVoxelsCCSub[j] = 0; + } + + numVoxelsCCSub[value]++; + ++subCCItr; + } + + unsigned int numSubLabels = numVoxelsCCSub.size() - 1; + double ratioOverlap = static_cast (numVoxelsCCSub[numSubLabels]) / labelsOverlap[i]; + bool okGrowing = (ratioOverlap > betaArg.getValue()); + + refItr.GoToBegin(); + testItr.GoToBegin(); + if (okGrowing) + { + while (!refItr.IsAtEnd()) + { + if (testItr.Get() == i) + { + if (refItr.Get() == 0) + testItr.Set(maxTestLabel + 3); + else + testItr.Set(maxTestLabel + 1); + } + + ++refItr; + ++testItr; + } + } + else + { + // Non growing lesion, setting to + while (!testItr.IsAtEnd()) + { + if (testItr.Get() == i) + testItr.Set(maxTestLabel + 1); + + ++testItr; + } + } + } + + testItr.GoToBegin(); + while (!testItr.IsAtEnd()) + { + unsigned int value = testItr.Get(); + if (value > 0) + testItr.Set(value - maxTestLabel); + + ++testItr; + } + + std::cout << "Writing output to " << outArg.getValue() << std::endl; + + anima::writeImage (outArg.getValue(),testSegmentation); + return EXIT_SUCCESS; +} + From 8d35e7e57f97a0682ff68fc6e1923087f3376180 Mon Sep 17 00:00:00 2001 From: Olivier Commowick Date: Thu, 14 Mar 2019 12:16:25 +0100 Subject: [PATCH 2/3] Add tool to packaging --- superbuild/BinariesPackaging.cmake.in | 1 + 1 file changed, 1 insertion(+) diff --git a/superbuild/BinariesPackaging.cmake.in b/superbuild/BinariesPackaging.cmake.in index 13767ee91..fe55947a4 100644 --- a/superbuild/BinariesPackaging.cmake.in +++ b/superbuild/BinariesPackaging.cmake.in @@ -41,6 +41,7 @@ set(BINARY_TOOLS_LIST animaInfluenceZones animaIsosurface animaKMeansClustering + animaLesionEvolutionDetection animaLinearTransformArithmetic animaLinearTransformToSVF animaLocalPatchCovarianceDistance From bad4e34e87e03631104bd4b8c29ac2104d012023 Mon Sep 17 00:00:00 2001 From: Olivier Commowick Date: Thu, 14 Mar 2019 12:20:32 +0100 Subject: [PATCH 3/3] Update description --- .../animaLesionEvolutionDetection.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Anima/segmentation/validation_tools/lesion_evolution_detection/animaLesionEvolutionDetection.cxx b/Anima/segmentation/validation_tools/lesion_evolution_detection/animaLesionEvolutionDetection.cxx index 078268b17..0c71d502d 100644 --- a/Anima/segmentation/validation_tools/lesion_evolution_detection/animaLesionEvolutionDetection.cxx +++ b/Anima/segmentation/validation_tools/lesion_evolution_detection/animaLesionEvolutionDetection.cxx @@ -12,7 +12,7 @@ int main(int argc, char **argv) { - TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + TCLAP::CmdLine cmd("Computes from two segmentations the connected components that grew (L=3), were already there (L=1), or are new (L=2).\n INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); TCLAP::ValueArg refArg("r", "ref", "Reference binary segmentation", true, "", "reference segmentation", cmd); TCLAP::ValueArg testArg("t", "test", "Test binary segmentation", true, "", "test segmentation", cmd);