diff --git a/CMake/BuildConfigurations/DiffusionAll.cmake b/CMake/BuildConfigurations/DiffusionAll.cmake index 22fd3320ec..b0b90a6643 100644 --- a/CMake/BuildConfigurations/DiffusionAll.cmake +++ b/CMake/BuildConfigurations/DiffusionAll.cmake @@ -1,39 +1,39 @@ message(STATUS "Configuring MITK Diffusion with all Plugins") # Enable non-optional external dependencies set(MITK_USE_Vigra ON CACHE BOOL "MITK Use Vigra Library" FORCE) set(MITK_USE_HDF5 ON CACHE BOOL "MITK Use HDF5 Library" FORCE) set(MITK_USE_MatchPoint ON CACHE BOOL "" FORCE) set(MITK_USE_DCMTK ON CACHE BOOL "" FORCE) set(MITK_USE_DCMQI OFF CACHE BOOL "" FORCE) set(MITK_USE_PCRE ON CACHE BOOL "" FORCE) set(MITK_USE_SWIG ON CACHE BOOL "" FORCE) set(MITK_USE_Python ON CACHE BOOL "" FORCE) set(MITK_USE_BetData ON CACHE BOOL "" FORCE) # Disable all apps but MITK Diffusion set(MITK_BUILD_ALL_APPS OFF CACHE BOOL "Build all MITK applications" FORCE) set(MITK_BUILD_APP_CoreApp OFF CACHE BOOL "Build the MITK CoreApp" FORCE) set(MITK_BUILD_APP_Diffusion ON CACHE BOOL "Build MITK Diffusion" FORCE) # Activate Diffusion Mini Apps -set(BUILD_DiffusionCoreCmdApps ON CACHE BOOL "Build commandline tools for diffusion" FORCE) set(BUILD_DiffusionFiberProcessingCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber processing" FORCE) set(BUILD_DiffusionFiberfoxCmdApps ON CACHE BOOL "Build commandline tools for diffusion data simulation (Fiberfox)" FORCE) set(BUILD_DiffusionMiscCmdApps ON CACHE BOOL "Build miscellaneous commandline tools for diffusion" FORCE) set(BUILD_DiffusionQuantificationCmdApps ON CACHE BOOL "Build commandline tools for diffusion quantification (IVIM, ADC, ...)" FORCE) set(BUILD_DiffusionTractographyCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber tractography" FORCE) set(BUILD_DiffusionTractographyEvaluationCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber tractography evaluation" FORCE) set(BUILD_DiffusionConnectomicsCmdApps ON CACHE BOOL "Build commandline tools for diffusion connectomics" FORCE) +set(BUILD_DiffusionPythonCmdApps ON CACHE BOOL "Build commandline tools for diffusion with python" FORCE) # Build neither all plugins nor examples set(MITK_BUILD_ALL_PLUGINS OFF CACHE BOOL "Build all MITK plugins" FORCE) set(MITK_BUILD_EXAMPLES OFF CACHE BOOL "Build the MITK examples" FORCE) set(BUILD_TESTING ON CACHE BOOL "Build the MITK tests" FORCE) # Activate in-application help generation set(MITK_DOXYGEN_GENERATE_QCH_FILES ON CACHE BOOL "Use doxygen to generate Qt compressed help files for MITK docs" FORCE) set(BLUEBERRY_USE_QT_HELP ON CACHE BOOL "Enable support for integrating bundle documentation into Qt Help" FORCE) # Disable console window set(MITK_SHOW_CONSOLE_WINDOW ON CACHE BOOL "Use this to enable or disable the console window when starting MITK GUI Applications" FORCE) diff --git a/CMake/BuildConfigurations/DiffusionRelease.cmake b/CMake/BuildConfigurations/DiffusionRelease.cmake index ef55857c2a..7bbb412d7f 100644 --- a/CMake/BuildConfigurations/DiffusionRelease.cmake +++ b/CMake/BuildConfigurations/DiffusionRelease.cmake @@ -1,42 +1,42 @@ message(STATUS "Configuring MITK Diffusion Release Build") # Enable non-optional external dependencies set(MITK_USE_Vigra ON CACHE BOOL "MITK Use Vigra Library" FORCE) set(MITK_USE_HDF5 ON CACHE BOOL "MITK Use HDF5 Library" FORCE) set(MITK_USE_MatchPoint ON CACHE BOOL "" FORCE) set(MITK_USE_DCMTK ON CACHE BOOL "" FORCE) set(MITK_USE_DCMQI OFF CACHE BOOL "" FORCE) set(MITK_USE_PCRE ON CACHE BOOL "" FORCE) set(MITK_USE_SWIG ON CACHE BOOL "" FORCE) set(MITK_USE_Python ON CACHE BOOL "" FORCE) set(MITK_USE_BetData ON CACHE BOOL "" FORCE) # Disable all apps but MITK Diffusion set(MITK_BUILD_ALL_APPS OFF CACHE BOOL "Build all MITK applications" FORCE) set(MITK_BUILD_APP_CoreApp OFF CACHE BOOL "Build the MITK CoreApp" FORCE) set(MITK_BUILD_APP_Workbench OFF CACHE BOOL "Build the MITK Workbench" FORCE) set(MITK_BUILD_APP_Diffusion ON CACHE BOOL "Build MITK Diffusion" FORCE) # Activate Diffusion Mini Apps -set(BUILD_DiffusionCoreCmdApps ON CACHE BOOL "Build commandline tools for diffusion" FORCE) set(BUILD_DiffusionFiberProcessingCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber processing" FORCE) set(BUILD_DiffusionFiberfoxCmdApps ON CACHE BOOL "Build commandline tools for diffusion data simulation (Fiberfox)" FORCE) set(BUILD_DiffusionMiscCmdApps ON CACHE BOOL "Build miscellaneous commandline tools for diffusion" FORCE) set(BUILD_DiffusionQuantificationCmdApps OFF CACHE BOOL "Build commandline tools for diffusion quantification (IVIM, ADC, ...)" FORCE) set(BUILD_DiffusionTractographyCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber tractography" FORCE) set(BUILD_DiffusionTractographyEvaluationCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber tractography evaluation" FORCE) set(BUILD_DiffusionConnectomicsCmdApps ON CACHE BOOL "Build commandline tools for diffusion connectomics" FORCE) +set(BUILD_DiffusionPythonCmdApps ON CACHE BOOL "Build commandline tools for diffusion with python" FORCE) # Build neither all plugins nor examples set(MITK_BUILD_ALL_PLUGINS OFF CACHE BOOL "Build all MITK plugins" FORCE) set(MITK_BUILD_EXAMPLES OFF CACHE BOOL "Build the MITK examples" FORCE) set(BUILD_TESTING OFF CACHE BOOL "Build the MITK tests" FORCE) # Activate in-application help generation set(MITK_DOXYGEN_GENERATE_QCH_FILES ON CACHE BOOL "Use doxygen to generate Qt compressed help files for MITK docs" FORCE) set(BLUEBERRY_USE_QT_HELP ON CACHE BOOL "Enable support for integrating bundle documentation into Qt Help" FORCE) # Enable console window set(MITK_SHOW_CONSOLE_WINDOW ON CACHE BOOL "Use this to enable or disable the console window when starting MITK GUI Applications" FORCE) set(CMAKE_BUILD_TYPE Release CACHE STRING "" FORCE) diff --git a/Modules/DiffusionImaging/CMakeLists.txt b/Modules/DiffusionImaging/CMakeLists.txt index bcd76551e5..0be1046bbc 100644 --- a/Modules/DiffusionImaging/CMakeLists.txt +++ b/Modules/DiffusionImaging/CMakeLists.txt @@ -1,13 +1,14 @@ set( diffusion_module_dirs DiffusionCore FiberTracking Connectomics Quantification DiffusionIO + DiffusionCmdApps ) foreach(diffusion_module_dir ${diffusion_module_dirs}) add_subdirectory(${diffusion_module_dir}) endforeach() configure_file(${CMAKE_CURRENT_SOURCE_DIR}/mitkDiffusionImagingConfigure.h.in ${CMAKE_CURRENT_BINARY_DIR}/DiffusionCore/mitkDiffusionImagingConfigure.h) diff --git a/Modules/DiffusionImaging/Connectomics/CMakeLists.txt b/Modules/DiffusionImaging/Connectomics/CMakeLists.txt index ae7887e933..3099bee377 100644 --- a/Modules/DiffusionImaging/Connectomics/CMakeLists.txt +++ b/Modules/DiffusionImaging/Connectomics/CMakeLists.txt @@ -1,9 +1,8 @@ MITK_CREATE_MODULE( SUBPROJECTS MITK-DTI INCLUDE_DIRS Algorithms Algorithms/BrainParcellation IODataStructures Rendering ${CMAKE_CURRENT_BINARY_DIR} DEPENDS MitkDiffusionCore MitkFiberTracking PACKAGE_DEPENDS VTK|vtkInfovisLayout ) add_subdirectory(Testing) -add_subdirectory(cmdapps) diff --git a/Modules/DiffusionImaging/DiffusionCmdApps/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCmdApps/CMakeLists.txt new file mode 100644 index 0000000000..ffbde72f1d --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCmdApps/CMakeLists.txt @@ -0,0 +1,20 @@ +MITK_CREATE_MODULE( + SUBPROJECTS MITK-DTI + # INCLUDE_DIRS Helpers + DEPENDS MitkDiffusionCore MitkDiffusionIO MitkFiberTracking MitkQuantification MitkConnectomics + PACKAGE_DEPENDS PUBLIC ITK +) + +if(MODULE_IS_ENABLED) + add_subdirectory(Connectomics) + add_subdirectory(Fiberfox) + add_subdirectory(FiberProcessing) + add_subdirectory(Tractography) + add_subdirectory(TractographyEvaluation) + add_subdirectory(Misc) + add_subdirectory(Quantification) + +if(MITK_USE_Python) + add_subdirectory(Python) +endif() +endif() diff --git a/Modules/DiffusionImaging/Connectomics/cmdapps/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/CMakeLists.txt similarity index 85% rename from Modules/DiffusionImaging/Connectomics/cmdapps/CMakeLists.txt rename to Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/CMakeLists.txt index c2050ff407..bcfc777262 100755 --- a/Modules/DiffusionImaging/Connectomics/cmdapps/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/CMakeLists.txt @@ -1,39 +1,39 @@ option(BUILD_DiffusionConnectomicsCmdApps "Build commandline tools for diffusion connectomics" OFF) if(BUILD_DiffusionConnectomicsCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of DiffusionConnectomics cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( DiffusionConnectomicscmdapps - NetworkCreation^^MitkFiberTracking_MitkConnectomics - NetworkStatistics^^MitkConnectomics + NetworkCreation^^ + NetworkStatistics^^ ) foreach(DiffusionConnectomicscmdapp ${DiffusionConnectomicscmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${DiffusionConnectomicscmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} - DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} - PACKAGE_DEPENDS ITK + DEPENDS MitkDiffusionCmdApps ${dependencies_list} + PACKAGE_DEPENDS ) endforeach() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionImaging/Connectomics/cmdapps/NetworkCreation.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkCreation.cpp similarity index 100% rename from Modules/DiffusionImaging/Connectomics/cmdapps/NetworkCreation.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkCreation.cpp diff --git a/Modules/DiffusionImaging/Connectomics/cmdapps/NetworkStatistics.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkStatistics.cpp similarity index 100% rename from Modules/DiffusionImaging/Connectomics/cmdapps/NetworkStatistics.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Connectomics/NetworkStatistics.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/CMakeLists.txt similarity index 66% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/CMakeLists.txt rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/CMakeLists.txt index f3d6510284..fdf23fbbd0 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/CMakeLists.txt @@ -1,49 +1,47 @@ option(BUILD_DiffusionFiberProcessingCmdApps "Build commandline tools for diffusion fiber processing" OFF) if(BUILD_DiffusionFiberProcessingCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionFiberProcessingcmdapps - TractDensity^^MitkFiberTracking - Sift2WeightCopy^^MitkFiberTracking - FiberExtraction^^MitkFiberTracking - FiberExtractionRoi^^MitkFiberTracking - FiberProcessing^^MitkFiberTracking - FitFibersToImage^^MitkFiberTracking - FiberDirectionExtraction^^MitkFiberTracking - FiberJoin^^MitkFiberTracking - FiberClustering^^MitkFiberTracking - GetOverlappingTracts^^MitkFiberTracking - TractDensityFilter^^MitkFiberTracking - TractDistance^^MitkFiberTracking + TractDensity^^ + Sift2WeightCopy^^ + FiberExtraction^^ + FiberExtractionRoi^^ + FiberProcessing^^ + FitFibersToImage^^ + FiberDirectionExtraction^^ + FiberJoin^^ + FiberClustering^^ + TractDensityFilter^^ ) foreach(diffusionFiberProcessingcmdapp ${diffusionFiberProcessingcmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${diffusionFiberProcessingcmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} - DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} - PACKAGE_DEPENDS ITK + DEPENDS MitkDiffusionCmdApps ${dependencies_list} + PACKAGE_DEPENDS ) endforeach() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberClustering.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberClustering.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberDirectionExtraction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberDirectionExtraction.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberDirectionExtraction.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberDirectionExtraction.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtraction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtraction.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtraction.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtraction.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtractionRoi.cpp similarity index 96% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtractionRoi.cpp index 693ed38cda..11be42c254 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberExtractionRoi.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberExtractionRoi.cpp @@ -1,271 +1,276 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; ItkFloatImgType::Pointer LoadItkImage(const std::string& filename) { mitk::Image::Pointer img = mitk::IOUtil::Load(filename); ItkFloatImgType::Pointer itk_image = ItkFloatImgType::New(); mitk::CastToItkImage(img, itk_image); return itk_image; } /*! \brief Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Extraction With ROI Image"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Extract fibers from a tractogram using binary image ROIs"); parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "input tractogram (.fib/.trk/.tck/.dcm)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::String, "Output:", "output root", us::Any(), false); parser.addArgument("rois", "", mitkCommandLineParser::StringList, "ROI images:", "ROI images", us::Any(), false); parser.endGroup(); parser.beginGroup("2. Label based extraction:"); - parser.addArgument("labels", "", mitkCommandLineParser::StringList, "Labels:", "positive means roi image value in labels vector", us::Any()); parser.addArgument("split_labels", "", mitkCommandLineParser::Bool, "Split labels:", "output a separate tractogram for each label-->label tract", false); parser.addArgument("skip_self_connections", "", mitkCommandLineParser::Bool, "Skip self connections:", "ignore streamlines between two identical labels", false); + parser.addArgument("all_labels", "", mitkCommandLineParser::Bool, "All labels:", "use all labels (0 is excluded)", false); + parser.addArgument("labels", "", mitkCommandLineParser::StringList, "Labels:", "positive means roi image value in labels vector", us::Any()); parser.addArgument("start_labels", "", mitkCommandLineParser::StringList, "Start Labels:", "use separate start and end labels instead of one mixed set", us::Any()); parser.addArgument("end_labels", "", mitkCommandLineParser::StringList, "End Labels:", "use separate start and end labels instead of one mixed set", us::Any()); parser.addArgument("paired", "", mitkCommandLineParser::Bool, "Paired:", "start and end label list are paired", false); parser.endGroup(); parser.beginGroup("3. Misc:"); parser.addArgument("both_ends", "", mitkCommandLineParser::Bool, "Both ends:", "Fibers are extracted if both endpoints are located in the ROI.", false); parser.addArgument("overlap_fraction", "", mitkCommandLineParser::Float, "Overlap fraction:", "Extract by overlap, not by endpoints. Extract fibers that overlap to at least the provided factor (0-1) with the ROI.", -1); parser.addArgument("invert", "", mitkCommandLineParser::Bool, "Invert:", "get streamlines not positive for any of the ROI images", false); parser.addArgument("output_negatives", "", mitkCommandLineParser::Bool, "Negatives:", "output negatives", false); parser.addArgument("interpolate", "", mitkCommandLineParser::Bool, "Interpolate:", "interpolate ROI images (only for endpoint based extraction)", false); parser.addArgument("threshold", "", mitkCommandLineParser::Float, "Threshold:", "positive means ROI image value threshold", 0.5); parser.addArgument("min_fibers", "", mitkCommandLineParser::Int, "Min. num. fibers:", "discard positive tracts with less fibers", 0); parser.addArgument("split_rois", "", mitkCommandLineParser::Bool, "Split ROIs:", "output a separate tractogram for each ROI", false); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFib = us::any_cast(parsedArgs["input"]); std::string outFib = us::any_cast(parsedArgs["out"]); mitkCommandLineParser::StringContainerType roi_files = us::any_cast(parsedArgs["rois"]); bool both_ends = false; if (parsedArgs.count("both_ends")) both_ends = us::any_cast(parsedArgs["both_ends"]); bool invert = false; if (parsedArgs.count("invert")) invert = us::any_cast(parsedArgs["invert"]); unsigned int min_fibers = 0; if (parsedArgs.count("min_fibers")) min_fibers = us::any_cast(parsedArgs["min_fibers"]); + bool all_labels = false; + if (parsedArgs.count("all_labels")) + all_labels = us::any_cast(parsedArgs["all_labels"]); + bool split_labels = false; if (parsedArgs.count("split_labels")) split_labels = us::any_cast(parsedArgs["split_labels"]); bool split_rois = false; if (parsedArgs.count("split_rois")) split_rois = us::any_cast(parsedArgs["split_rois"]); bool skip_self_connections = false; if (parsedArgs.count("skip_self_connections")) skip_self_connections = us::any_cast(parsedArgs["skip_self_connections"]); bool output_negatives = false; if (parsedArgs.count("output_negatives")) output_negatives = us::any_cast(parsedArgs["output_negatives"]); float overlap_fraction = -1; if (parsedArgs.count("overlap_fraction")) overlap_fraction = us::any_cast(parsedArgs["overlap_fraction"]); bool any_point = false; if (overlap_fraction>=0) any_point = true; bool interpolate = false; if (parsedArgs.count("interpolate")) interpolate = us::any_cast(parsedArgs["interpolate"]); float threshold = 0.5; if (parsedArgs.count("threshold")) threshold = us::any_cast(parsedArgs["threshold"]); mitkCommandLineParser::StringContainerType labels; if (parsedArgs.count("labels")) labels = us::any_cast(parsedArgs["labels"]); mitkCommandLineParser::StringContainerType start_labels; if (parsedArgs.count("start_labels")) start_labels = us::any_cast(parsedArgs["start_labels"]); mitkCommandLineParser::StringContainerType end_labels; if (parsedArgs.count("end_labels")) end_labels = us::any_cast(parsedArgs["end_labels"]); bool paired = false; if (parsedArgs.count("paired")) paired = us::any_cast(parsedArgs["paired"]); try { // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = mitk::IOUtil::Load(inFib); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect std::vector< ItkFloatImgType::Pointer > roi_images; std::vector< std::string > roi_names; for (std::size_t i=0; i roi_images2; for (auto roi : roi_images) roi_images2.push_back(roi); std::vector< unsigned short > short_labels; for (auto l : labels) short_labels.push_back(boost::lexical_cast(l)); std::vector< unsigned short > short_start_labels; for (auto l : start_labels) short_start_labels.push_back(boost::lexical_cast(l)); std::vector< unsigned short > short_end_labels; for (auto l : end_labels) short_end_labels.push_back(boost::lexical_cast(l)); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(inputTractogram); extractor->SetRoiImages(roi_images2); extractor->SetRoiImageNames(roi_names); extractor->SetOverlapFraction(overlap_fraction); extractor->SetBothEnds(both_ends); extractor->SetInterpolate(interpolate); extractor->SetThreshold(threshold); extractor->SetLabels(short_labels); extractor->SetStartLabels(short_start_labels); extractor->SetEndLabels(short_end_labels); extractor->SetSplitLabels(split_labels); extractor->SetSplitByRoi(split_rois); extractor->SetMinFibersPerTract(min_fibers); extractor->SetSkipSelfConnections(skip_self_connections); extractor->SetPairedStartEndLabels(paired); if (!any_point) extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); - if (short_labels.size()>0 || short_start_labels.size()>0 || short_end_labels.size()>0) + if (all_labels || short_labels.size()>0 || short_start_labels.size()>0 || short_end_labels.size()>0) extractor->SetInputType(itk::FiberExtractionFilter::INPUT::LABEL_MAP); extractor->Update(); std::string ext = itksys::SystemTools::GetFilenameExtension(outFib); if (ext.empty()) ext = ".trk"; outFib = itksys::SystemTools::GetFilenamePath(outFib) + '/' + itksys::SystemTools::GetFilenameWithoutExtension(outFib); if (invert) mitk::IOUtil::Save(extractor->GetNegatives().at(0), outFib + ext); else { int c = 0; std::vector< std::string > positive_labels = extractor->GetPositiveLabels(); for (auto fib : extractor->GetPositives()) { std::string l = positive_labels.at(c); if (l.size()>0) mitk::IOUtil::Save(fib, outFib + "_" + l + ext); else mitk::IOUtil::Save(fib, outFib + ext); ++c; } } if (output_negatives) { invert = !invert; if (invert) mitk::IOUtil::Save(extractor->GetNegatives().at(0), outFib + "_negatives" + ext); else { int c = 0; std::vector< std::string > positive_labels = extractor->GetPositiveLabels(); for (auto fib : extractor->GetPositives()) { std::string l = positive_labels.at(c); if (l.size()>0) mitk::IOUtil::Save(fib, outFib + "_" + l + "_negatives" + ext); else mitk::IOUtil::Save(fib, outFib + "_negatives" + ext); ++c; } } } } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberJoin.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberJoin.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberJoin.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberJoin.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberProcessing.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberProcessing.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberProcessing.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FiberProcessing.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FitFibersToImage.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp similarity index 88% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FitFibersToImage.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp index ff153e017a..ca8e8e9bbe 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FitFibersToImage.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp @@ -1,346 +1,299 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include -#include -#include #include #include #include #include #include #include #include +#include -typedef itksys::SystemTools ist; typedef itk::Point PointType4; typedef itk::Image< float, 4 > PeakImgType; -std::vector< std::string > get_file_list(const std::string& path) -{ - std::vector< std::string > file_list; - itk::Directory::Pointer dir = itk::Directory::New(); - - if (dir->Load(path.c_str())) - { - int n = dir->GetNumberOfFiles(); - for (int r = 0; r < n; r++) - { - const char *filename = dir->GetFile(r); - std::string ext = ist::GetFilenameExtension(filename); - if (ext==".fib" || ext==".trk") - file_list.push_back(path + '/' + filename); - } - } - return file_list; -} - /*! \brief Fits the tractogram to the input peak image by assigning a weight to each fiber (similar to https://doi.org/10.1016/j.neuroimage.2015.06.092). */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fit Fibers To Image"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Assigns a weight to each fiber in order to optimally explain the input peak image"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i1", mitkCommandLineParser::StringList, "Input tractograms:", "input tractograms (.fib, vtk ascii file format)", us::Any(), false); parser.addArgument("", "i2", mitkCommandLineParser::InputFile, "Input image:", "input image", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); parser.addArgument("max_iter", "", mitkCommandLineParser::Int, "Max. iterations:", "maximum number of optimizer iterations", 20); parser.addArgument("bundle_based", "", mitkCommandLineParser::Bool, "Bundle based fit:", "fit one weight per input tractogram/bundle, not for each fiber", false); parser.addArgument("min_g", "", mitkCommandLineParser::Float, "Min. g:", "lower termination threshold for gradient magnitude", 1e-5); parser.addArgument("lambda", "", mitkCommandLineParser::Float, "Lambda:", "modifier for regularization", 0.1); parser.addArgument("save_res", "", mitkCommandLineParser::Bool, "Save Residuals:", "save residual images", false); parser.addArgument("save_weights", "", mitkCommandLineParser::Bool, "Save Weights:", "save fiber weights in a separate text file", false); parser.addArgument("filter_outliers", "", mitkCommandLineParser::Bool, "Filter outliers:", "perform second optimization run with an upper weight bound based on the first weight estimation (99% quantile)", false); parser.addArgument("join_tracts", "", mitkCommandLineParser::Bool, "Join output tracts:", "outout tracts are merged into a single tractogram", false); parser.addArgument("regu", "", mitkCommandLineParser::String, "Regularization:", "MSM, Variance, VoxelVariance (default), Lasso, GroupLasso, GroupVariance, NONE"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType fib_files = us::any_cast(parsedArgs["i1"]); std::string input_image_name = us::any_cast(parsedArgs["i2"]); std::string outRoot = us::any_cast(parsedArgs["o"]); bool single_fib = true; if (parsedArgs.count("bundle_based")) single_fib = !us::any_cast(parsedArgs["bundle_based"]); bool save_residuals = false; if (parsedArgs.count("save_res")) save_residuals = us::any_cast(parsedArgs["save_res"]); bool save_weights = false; if (parsedArgs.count("save_weights")) save_weights = us::any_cast(parsedArgs["save_weights"]); std::string regu = "VoxelVariance"; if (parsedArgs.count("regu")) regu = us::any_cast(parsedArgs["regu"]); bool join_tracts = false; if (parsedArgs.count("join_tracts")) join_tracts = us::any_cast(parsedArgs["join_tracts"]); int max_iter = 20; if (parsedArgs.count("max_iter")) max_iter = us::any_cast(parsedArgs["max_iter"]); float g_tol = 1e-5; if (parsedArgs.count("min_g")) g_tol = us::any_cast(parsedArgs["min_g"]); float lambda = 0.1; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); bool filter_outliers = false; if (parsedArgs.count("filter_outliers")) filter_outliers = us::any_cast(parsedArgs["filter_outliers"]); try { MITK_INFO << "Loading data"; -// std::streambuf *old = cout.rdbuf(); // <-- save -// std::stringstream ss; -// std::cout.rdbuf (ss.rdbuf()); // <-- redirect - std::vector< mitk::FiberBundle::Pointer > input_tracts; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, {}); std::vector< std::string > fib_names; - for (auto item : fib_files) - { - if ( ist::FileIsDirectory(item) ) - { - for ( auto fibFile : get_file_list(item) ) - { - mitk::FiberBundle::Pointer inputTractogram = mitk::IOUtil::Load(fibFile); - if (inputTractogram.IsNull()) - continue; - input_tracts.push_back(inputTractogram); - fib_names.push_back(fibFile); - } - } - else - { - mitk::FiberBundle::Pointer inputTractogram = mitk::IOUtil::Load(item); - if (inputTractogram.IsNull()) - continue; - input_tracts.push_back(inputTractogram); - fib_names.push_back(item); - } - } -// std::cout.rdbuf (old); // <-- restore + auto input_tracts = mitk::DiffusionDataIOHelper::load_fibs(fib_files, &fib_names); itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); mitk::BaseData::Pointer inputData = mitk::IOUtil::Load(input_image_name, &functor)[0].GetPointer(); mitk::Image::Pointer mitk_image = dynamic_cast(inputData.GetPointer()); mitk::PeakImage::Pointer mitk_peak_image = dynamic_cast(inputData.GetPointer()); if (mitk_peak_image.IsNotNull()) { typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitk_peak_image); caster->Update(); mitk::PeakImage::ItkPeakImageType::Pointer peak_image = caster->GetOutput(); fitter->SetPeakImage(peak_image); } else if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(mitk_image)) { fitter->SetDiffImage(mitk::DiffusionPropertyHelper::GetItkVectorImage(mitk_image)); mitk::TensorModel<>* model = new mitk::TensorModel<>(); model->SetBvalue(1000); model->SetDiffusivity1(0.0010); model->SetDiffusivity2(0.00015); model->SetDiffusivity3(0.00015); model->SetGradientList(mitk::DiffusionPropertyHelper::GetGradientContainer(mitk_image)); fitter->SetSignalModel(model); } else if (mitk_image->GetDimension()==3) { itk::FitFibersToImageFilter::DoubleImgType::Pointer scalar_image = itk::FitFibersToImageFilter::DoubleImgType::New(); mitk::CastToItkImage(mitk_image, scalar_image); fitter->SetScalarImage(scalar_image); } else { MITK_INFO << "Input image invalid. Valid options are peak image, 3D scalar image or raw diffusion-weighted image."; return EXIT_FAILURE; } fitter->SetTractograms(input_tracts); fitter->SetFitIndividualFibers(single_fib); fitter->SetMaxIterations(max_iter); fitter->SetGradientTolerance(g_tol); fitter->SetLambda(lambda); fitter->SetFilterOutliers(filter_outliers); if (regu=="MSM") fitter->SetRegularization(VnlCostFunction::REGU::MSM); else if (regu=="Variance") fitter->SetRegularization(VnlCostFunction::REGU::VARIANCE); else if (regu=="Lasso") fitter->SetRegularization(VnlCostFunction::REGU::LASSO); else if (regu=="VoxelVariance") fitter->SetRegularization(VnlCostFunction::REGU::VOXEL_VARIANCE); else if (regu=="GroupLasso") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_LASSO); else if (regu=="GroupVariance") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_VARIANCE); else if (regu=="NONE") fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); if (save_residuals && mitk_peak_image.IsNotNull()) { itk::ImageFileWriter< PeakImgType >::Pointer writer = itk::ImageFileWriter< PeakImgType >::New(); writer->SetInput(fitter->GetFittedImage()); writer->SetFileName(outRoot + "_fitted.nii.gz"); writer->Update(); writer->SetInput(fitter->GetResidualImage()); writer->SetFileName(outRoot + "_residual.nii.gz"); writer->Update(); writer->SetInput(fitter->GetOverexplainedImage()); writer->SetFileName(outRoot + "_overexplained.nii.gz"); writer->Update(); writer->SetInput(fitter->GetUnderexplainedImage()); writer->SetFileName(outRoot + "_underexplained.nii.gz"); writer->Update(); } else if (save_residuals && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(mitk_image)) { { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetFittedImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_fitted_image.nii.gz"); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetResidualImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_residual_image.nii.gz"); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetOverexplainedImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_overexplained_image.nii.gz"); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetUnderexplainedImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_underexplained_image.nii.gz"); } } else if (save_residuals) { itk::ImageFileWriter< itk::FitFibersToImageFilter::DoubleImgType >::Pointer writer = itk::ImageFileWriter< itk::FitFibersToImageFilter::DoubleImgType >::New(); writer->SetInput(fitter->GetFittedImageScalar()); writer->SetFileName(outRoot + "_fitted_image.nii.gz"); writer->Update(); writer->SetInput(fitter->GetResidualImageScalar()); writer->SetFileName(outRoot + "_residual_image.nii.gz"); writer->Update(); writer->SetInput(fitter->GetOverexplainedImageScalar()); writer->SetFileName(outRoot + "_overexplained_image.nii.gz"); writer->Update(); writer->SetInput(fitter->GetUnderexplainedImageScalar()); writer->SetFileName(outRoot + "_underexplained_image.nii.gz"); writer->Update(); } std::vector< mitk::FiberBundle::Pointer > output_tracts = fitter->GetTractograms(); if (!join_tracts) { for (unsigned int bundle=0; bundleGetNumFibers(); ++f) logfile << output_tracts.at(bundle)->GetFiberWeight(f) << "\n"; logfile.close(); } } } else { mitk::FiberBundle::Pointer out = mitk::FiberBundle::New(); out = out->AddBundles(output_tracts); out->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(out, outRoot + "_fitted.fib"); if (save_weights) { ofstream logfile; logfile.open (outRoot + "_weights.txt"); for (int f=0; fGetNumFibers(); ++f) logfile << out->GetFiberWeight(f) << "\n"; logfile.close(); } } } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/Sift2WeightCopy.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/Sift2WeightCopy.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/Sift2WeightCopy.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/Sift2WeightCopy.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/TractDensity.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/TractDensity.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/TractDensity.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/TractDensity.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/TractDensityFilter.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/TractDensityFilter.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/TractDensityFilter.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/FiberProcessing/TractDensityFilter.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/CMakeLists.txt similarity index 86% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/CMakeLists.txt rename to Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/CMakeLists.txt index 74361ab3e7..26dc2d25da 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/CMakeLists.txt @@ -1,39 +1,39 @@ option(BUILD_DiffusionFiberfoxCmdApps "Build commandline tools for diffusion data simulation (Fiberfox)" OFF) if(BUILD_DiffusionFiberfoxCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionFiberfoxcmdapps - Fiberfox^^MitkFiberTracking - FiberfoxOptimization^^MitkFiberTracking + Fiberfox^^ + FiberfoxOptimization^^ ) foreach(diffusionFiberfoxcmdapp ${diffusionFiberfoxcmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${diffusionFiberfoxcmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} - DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} - PACKAGE_DEPENDS ITK + DEPENDS MitkDiffusionCmdApps ${dependencies_list} + PACKAGE_DEPENDS ) endforeach() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/Fiberfox.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/Fiberfox.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/Fiberfox.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/Fiberfox.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/FiberfoxOptimization.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/FiberfoxOptimization.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Fiberfox/FiberfoxOptimization.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Fiberfox/FiberfoxOptimization.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/CMakeLists.txt similarity index 71% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/CMakeLists.txt rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/CMakeLists.txt index 8342f17ffc..d4746a8381 100644 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/CMakeLists.txt @@ -1,37 +1,48 @@ option(BUILD_DiffusionMiscCmdApps "Build miscellaneous commandline tools for diffusion" OFF) if(BUILD_DiffusionMiscCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionMisccmdapps - PeakExtraction^^MitkFiberTracking - FileFormatConverter^^MitkFiberTracking - FlipPeaks^^MitkFiberTracking + PeakExtraction^^ + FileFormatConverter^^ + FlipPeaks^^ + ImageResampler^^ + CopyGeometry^^ + Registration^^ + DiffusionDICOMLoader^^ + ResampleGradients^^ + ShToOdfImage^^ + DImp^^ + DReg^^ + HeadMotionCorrection^^ + DmriDenoising^^MitkImageDenoising + RoundBvalues^^ ) foreach(diffusionMisccmdapp ${diffusionMisccmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${diffusionMisccmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} - DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} - PACKAGE_DEPENDS ITK + DEPENDS MitkDiffusionCmdApps ${dependencies_list} + PACKAGE_DEPENDS ) endforeach() endif() diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/CopyGeometry.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/CopyGeometry.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/CopyGeometry.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/CopyGeometry.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DICOMLoader.cmake b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DICOMLoader.cmake similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/DICOMLoader.cmake rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/DICOMLoader.cmake diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DImp.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DImp.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/DImp.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/DImp.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DReg.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DReg.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/DReg.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/DReg.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionDICOMLoader.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DiffusionDICOMLoader.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionDICOMLoader.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/DiffusionDICOMLoader.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DmriDenoising.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/DmriDenoising.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/DmriDenoising.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/DmriDenoising.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/FileFormatConverter.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/FileFormatConverter.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/FileFormatConverter.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/FileFormatConverter.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/FlipPeaks.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/FlipPeaks.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/FlipPeaks.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/FlipPeaks.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/HeadMotionCorrection.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/HeadMotionCorrection.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/HeadMotionCorrection.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/HeadMotionCorrection.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/ImageResampler.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ImageResampler.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/ImageResampler.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/ImageResampler.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/PeakExtraction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/PeakExtraction.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/PeakExtraction.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/PeakExtraction.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/Registration.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/Registration.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/Registration.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/Registration.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/ResampleGradients.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ResampleGradients.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/ResampleGradients.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/ResampleGradients.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/RoundBvalues.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/RoundBvalues.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/RoundBvalues.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/RoundBvalues.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/ShToOdfImage.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Misc/ShToOdfImage.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/ShToOdfImage.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Misc/ShToOdfImage.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps_python/BrainExtraction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Python/BrainExtraction.cpp similarity index 100% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps_python/BrainExtraction.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Python/BrainExtraction.cpp diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps_python/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCmdApps/Python/CMakeLists.txt similarity index 89% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps_python/CMakeLists.txt rename to Modules/DiffusionImaging/DiffusionCmdApps/Python/CMakeLists.txt index 5cbdd74a2f..8cfbef2cff 100644 --- a/Modules/DiffusionImaging/DiffusionCore/cmdapps_python/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Python/CMakeLists.txt @@ -1,34 +1,34 @@ option(BUILD_DiffusionPythonCmdApps "Build commandline tools for diffusion with python" OFF) if(BUILD_DiffusionPythonCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionpythoncmdapps BrainExtraction^^ ) foreach(diffusionpythoncmdapp ${diffusionpythoncmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${diffusionpythoncmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} - DEPENDS CppMicroServices MitkCore MitkDiffusionCore MitkPython ${dependencies_list} - PACKAGE_DEPENDS ITK + DEPENDS CppMicroServices MitkDiffusionCmdApps MitkPython ${dependencies_list} + PACKAGE_DEPENDS ) endforeach() endif() diff --git a/Modules/DiffusionImaging/Quantification/cmdapps/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/CMakeLists.txt similarity index 92% rename from Modules/DiffusionImaging/Quantification/cmdapps/CMakeLists.txt rename to Modules/DiffusionImaging/DiffusionCmdApps/Quantification/CMakeLists.txt index bd65206b21..d21f29655d 100644 --- a/Modules/DiffusionImaging/Quantification/cmdapps/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/CMakeLists.txt @@ -1,39 +1,39 @@ option(BUILD_DiffusionQuantificationCmdApps "Build commandline tools for diffusion quantification (IVIM, ADC, ...)" OFF) if(BUILD_DiffusionQuantificationCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionQuantificationcmdapps DiffusionIndices^^ QballReconstruction^^ TensorReconstruction^^ TensorDerivedMapsExtraction^^ DiffusionKurtosisFit^^ MultishellMethods^^ ) foreach(diffusionQuantificationcmdapp ${diffusionQuantificationcmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${diffusionQuantificationcmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} - DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} - PACKAGE_DEPENDS ITK + DEPENDS MitkDiffusionCmdApps ${dependencies_list} + PACKAGE_DEPENDS ) endforeach() endif() diff --git a/Modules/DiffusionImaging/Quantification/cmdapps/DiffusionIndices.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/DiffusionIndices.cpp similarity index 100% rename from Modules/DiffusionImaging/Quantification/cmdapps/DiffusionIndices.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Quantification/DiffusionIndices.cpp diff --git a/Modules/DiffusionImaging/Quantification/cmdapps/DiffusionKurtosisFit.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/DiffusionKurtosisFit.cpp similarity index 100% rename from Modules/DiffusionImaging/Quantification/cmdapps/DiffusionKurtosisFit.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Quantification/DiffusionKurtosisFit.cpp diff --git a/Modules/DiffusionImaging/Quantification/cmdapps/MultishellMethods.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/MultishellMethods.cpp similarity index 100% rename from Modules/DiffusionImaging/Quantification/cmdapps/MultishellMethods.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Quantification/MultishellMethods.cpp diff --git a/Modules/DiffusionImaging/Quantification/cmdapps/QballReconstruction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/QballReconstruction.cpp similarity index 100% rename from Modules/DiffusionImaging/Quantification/cmdapps/QballReconstruction.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Quantification/QballReconstruction.cpp diff --git a/Modules/DiffusionImaging/Quantification/cmdapps/TensorDerivedMapsExtraction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/TensorDerivedMapsExtraction.cpp similarity index 100% rename from Modules/DiffusionImaging/Quantification/cmdapps/TensorDerivedMapsExtraction.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Quantification/TensorDerivedMapsExtraction.cpp diff --git a/Modules/DiffusionImaging/Quantification/cmdapps/TensorReconstruction.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Quantification/TensorReconstruction.cpp similarity index 100% rename from Modules/DiffusionImaging/Quantification/cmdapps/TensorReconstruction.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Quantification/TensorReconstruction.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/CMakeLists.txt similarity index 83% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/CMakeLists.txt rename to Modules/DiffusionImaging/DiffusionCmdApps/Tractography/CMakeLists.txt index f8ab7cfa92..a14ed6d992 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/CMakeLists.txt @@ -1,40 +1,40 @@ option(BUILD_DiffusionTractographyCmdApps "Build commandline tools for diffusion fiber tractography" OFF) if(BUILD_DiffusionTractographyCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusiontractographycmdapps - StreamlineTractography^^MitkFiberTracking - GlobalTractography^^MitkFiberTracking - RfTraining^^MitkFiberTracking + StreamlineTractography^^ + GlobalTractography^^ + RfTraining^^ ) foreach(diffusiontractographycmdapp ${diffusiontractographycmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${diffusiontractographycmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} - DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} - PACKAGE_DEPENDS ITK + DEPENDS MitkDiffusionCmdApps ${dependencies_list} + PACKAGE_DEPENDS ) endforeach() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/GlobalTractography.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/GlobalTractography.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/GlobalTractography.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Tractography/GlobalTractography.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/RfTraining.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/RfTraining.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/RfTraining.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Tractography/RfTraining.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/StreamlineTractography.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/Tractography/StreamlineTractography.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp similarity index 61% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp index 885e170bc1..8198a62615 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp @@ -1,694 +1,446 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include -typedef itksys::SystemTools ist; typedef itk::Point PointType4; typedef itk::Image< float, 4 > PeakImgType; typedef itk::Image< unsigned char, 3 > ItkUcharImageType; -std::vector< mitk::FiberBundle::Pointer > CombineTractograms(std::vector< mitk::FiberBundle::Pointer > reference, std::vector< mitk::FiberBundle::Pointer > candidates, int skip=-1) -{ - std::vector< mitk::FiberBundle::Pointer > fib; - for (auto f : reference) - fib.push_back(f); - - int c = 0; - for (auto f : candidates) - { - if (c!=skip) - fib.push_back(f); - ++c; - } - - return fib; -} - -std::vector< std::string > get_file_list(const std::string& path, std::vector< std::string > extensions={".fib", ".trk"}) -{ - std::vector< std::string > file_list; - itk::Directory::Pointer dir = itk::Directory::New(); - - if (dir->Load(path.c_str())) - { - int n = dir->GetNumberOfFiles(); - for (int r = 0; r < n; r++) - { - const char *filename = dir->GetFile(r); - std::string ext = ist::GetFilenameExtension(filename); - for (auto e : extensions) - { - if (ext==e) - { - file_list.push_back(path + '/' + filename); - break; - } - } - } - } - std::sort(file_list.begin(), file_list.end()); - return file_list; -} - /*! \brief Score input candidate tracts using ACP analysis */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Anchor Constrained Plausibility"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription("Score input candidate tracts using ACP analysis"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "a", mitkCommandLineParser::InputFile, "Anchor tractogram:", "anchor tracts in one tractogram file", us::Any()); parser.addArgument("", "p", mitkCommandLineParser::InputFile, "Input peaks:", "input peak image", us::Any(), false); - parser.addArgument("", "c", mitkCommandLineParser::InputDirectory, "Candidates folder:", "folder containing candidate tracts", us::Any(), false); + parser.addArgument("", "c", mitkCommandLineParser::StringList, "Candidates:", "Folder(s) or file list of candidate tracts", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output folder:", "output folder", us::Any(), false); - parser.addArgument("reference_mask_folders", "", mitkCommandLineParser::StringList, "Reference Mask Folder(s):", "Folder(s) containing reference tract masks for accuracy evaluation"); - parser.addArgument("reference_peaks_folders", "", mitkCommandLineParser::StringList, "Reference Peaks Folder(s):", "Folder(s) containing reference peak images for accuracy evaluation"); + parser.addArgument("reference_mask_folders", "", mitkCommandLineParser::StringList, "Reference Mask Folder(s):", "Folder(s) or file list containing reference tract masks for accuracy evaluation"); + parser.addArgument("reference_peaks_folders", "", mitkCommandLineParser::StringList, "Reference Peaks Folder(s):", "Folder(s) or file list containing reference peak images for accuracy evaluation"); parser.addArgument("mask", "", mitkCommandLineParser::InputFile, "Mask image:", "scoring is only performed inside the mask image"); parser.addArgument("greedy_add", "", mitkCommandLineParser::Bool, "Greedy:", "if enabled, the candidate tracts are not jointly fitted to the residual image but one after the other employing a greedy scheme", false); parser.addArgument("lambda", "", mitkCommandLineParser::Float, "Lambda:", "modifier for regularization", 0.1); parser.addArgument("filter_outliers", "", mitkCommandLineParser::Bool, "Filter outliers:", "perform second optimization run with an upper weight bound based on the first weight estimation (99% quantile)", false); parser.addArgument("regu", "", mitkCommandLineParser::String, "Regularization:", "MSM, Variance, VoxelVariance, Lasso, GroupLasso, GroupVariance, NONE (default)"); parser.addArgument("use_num_streamlines", "", mitkCommandLineParser::Bool, "Use number of streamlines as score:", "Don't fit candidates, simply use number of streamlines per candidate as score", false); parser.addArgument("use_weights", "", mitkCommandLineParser::Bool, "Use input weights as score:", "Don't fit candidates, simply use first input streamline weight per candidate as score", false); parser.addArgument("filter_zero_weights", "", mitkCommandLineParser::Bool, "Filter zero-weights", "Remove streamlines with weight 0 from candidates", false); parser.addArgument("flipx", "", mitkCommandLineParser::Bool, "Flip x", "flip along x-axis", false); parser.addArgument("flipy", "", mitkCommandLineParser::Bool, "Flip y", "flip along y-axis", false); parser.addArgument("flipz", "", mitkCommandLineParser::Bool, "Flip z", "flip along z-axis", false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string peak_file_name = us::any_cast(parsedArgs["p"]); - std::string candidate_tract_folder = us::any_cast(parsedArgs["c"]); std::string out_folder = us::any_cast(parsedArgs["o"]); + mitkCommandLineParser::StringContainerType candidate_tract_folders = us::any_cast(parsedArgs["c"]); + if (!out_folder.empty() && out_folder.back() != '/') out_folder += "/"; - if (!candidate_tract_folder.empty() && candidate_tract_folder.back() != '/') - candidate_tract_folder += "/"; bool greedy_add = false; if (parsedArgs.count("greedy_add")) greedy_add = us::any_cast(parsedArgs["greedy_add"]); float lambda = 0.1; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); bool filter_outliers = false; if (parsedArgs.count("filter_outliers")) filter_outliers = us::any_cast(parsedArgs["filter_outliers"]); bool filter_zero_weights = false; if (parsedArgs.count("filter_zero_weights")) filter_zero_weights = us::any_cast(parsedArgs["filter_zero_weights"]); std::string mask_file = ""; if (parsedArgs.count("mask")) mask_file = us::any_cast(parsedArgs["mask"]); mitkCommandLineParser::StringContainerType reference_mask_files_folders; if (parsedArgs.count("reference_mask_folders")) reference_mask_files_folders = us::any_cast(parsedArgs["reference_mask_folders"]); mitkCommandLineParser::StringContainerType reference_peaks_files_folders; if (parsedArgs.count("reference_peaks_folders")) reference_peaks_files_folders = us::any_cast(parsedArgs["reference_peaks_folders"]); std::string regu = "NONE"; if (parsedArgs.count("regu")) regu = us::any_cast(parsedArgs["regu"]); bool use_weights = false; if (parsedArgs.count("use_weights")) use_weights = us::any_cast(parsedArgs["use_weights"]); bool use_num_streamlines = false; if (parsedArgs.count("use_num_streamlines")) use_num_streamlines = us::any_cast(parsedArgs["use_num_streamlines"]); bool flipx = false; if (parsedArgs.count("flipx")) flipx = us::any_cast(parsedArgs["flipx"]); bool flipy = false; if (parsedArgs.count("flipy")) flipy = us::any_cast(parsedArgs["flipy"]); bool flipz = false; if (parsedArgs.count("flipz")) flipz = us::any_cast(parsedArgs["flipz"]); try { itk::TimeProbe clock; clock.Start(); if (!ist::PathExists(out_folder)) { MITK_INFO << "Creating output directory"; ist::MakeDirectory(out_folder); } MITK_INFO << "Loading data"; - std::streambuf *old = cout.rdbuf(); // <-- save - std::stringstream ss; - std::cout.rdbuf (ss.rdbuf()); // <-- redirect - - ofstream logfile; - logfile.open (out_folder + "log.txt"); - logfile << "V3\n"; - - itk::ImageFileWriter< PeakImgType >::Pointer peak_image_writer = itk::ImageFileWriter< PeakImgType >::New(); - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, {}); - mitk::Image::Pointer inputImage = dynamic_cast(mitk::IOUtil::Load(peak_file_name, &functor)[0].GetPointer()); // Load mask file. Fit is only performed inside the mask - itk::FitFibersToImageFilter::UcharImgType::Pointer mask = nullptr; - if (mask_file.compare("")!=0) - { - mitk::Image::Pointer mitk_mask = mitk::IOUtil::Load(mask_file); - mitk::CastToItkImage(mitk_mask, mask); - } + auto mask = mitk::DiffusionDataIOHelper::load_itk_image(mask_file); // Load masks covering the true positives for evaluation purposes - std::vector< itk::FitFibersToImageFilter::UcharImgType::Pointer > reference_masks; std::vector< std::string > anchor_mask_files; - for (auto filename : reference_mask_files_folders) - { - if (itksys::SystemTools::PathExists(filename)) - { - if (!filename.empty() && filename.back() != '/') - filename += "/"; - - auto list = get_file_list(filename, {".nrrd",".nii.gz",".nii"}); - for (auto f : list) - { - MITK_INFO << f; - itk::FitFibersToImageFilter::UcharImgType::Pointer ref_mask = nullptr; - mitk::Image::Pointer ref_mitk_mask = mitk::IOUtil::Load(f); - mitk::CastToItkImage(ref_mitk_mask, ref_mask); - reference_masks.push_back(ref_mask); - anchor_mask_files.push_back(f); - } - } - else if (itksys::SystemTools::FileExists(filename)) - { - anchor_mask_files.push_back(filename); - itk::FitFibersToImageFilter::UcharImgType::Pointer ref_mask = nullptr; - mitk::Image::Pointer ref_mitk_mask = mitk::IOUtil::Load(filename); - mitk::CastToItkImage(ref_mitk_mask, ref_mask); - reference_masks.push_back(ref_mask); - } - } - - typedef mitk::ImageToItk< PeakImgType > CasterType; - std::vector< PeakImgType::Pointer > reference_peaks; - for (auto filename : reference_peaks_files_folders) - { - MITK_INFO << filename; - if (itksys::SystemTools::PathExists(filename)) - { - if (!filename.empty() && filename.back() != '/') - filename += "/"; - - auto list = get_file_list(filename, {".nrrd",".nii.gz",".nii"}); - for (auto f : list) - { - mitk::Image::Pointer ref_mitk_peaks = mitk::IOUtil::Load(f); - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(ref_mitk_peaks); - caster->Update(); - PeakImgType::Pointer peak_image = caster->GetOutput(); - reference_peaks.push_back(peak_image); - } - } - else if (itksys::SystemTools::FileExists(filename)) - { - mitk::Image::Pointer ref_mitk_peaks = mitk::IOUtil::Load(filename); - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(ref_mitk_peaks); - caster->Update(); - PeakImgType::Pointer peak_image = caster->GetOutput(); - reference_peaks.push_back(peak_image); - } - } + auto reference_masks = mitk::DiffusionDataIOHelper::load_itk_images(reference_mask_files_folders, &anchor_mask_files); + auto reference_peaks = mitk::DiffusionDataIOHelper::load_itk_images(reference_peaks_files_folders); // Load peak image - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(inputImage); - caster->Update(); - PeakImgType::Pointer peak_image = caster->GetOutput(); + auto peak_image = mitk::DiffusionDataIOHelper::load_itk_image(peak_file_name); // Load all candidate tracts - std::vector< std::string > candidate_tract_files = get_file_list(candidate_tract_folder); - std::vector< mitk::FiberBundle::Pointer > input_candidates; - for (std::string f : candidate_tract_files) - { - mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(f); - if (fib.IsNull()) - continue; - if (fib->GetNumFibers()<=0) - continue; - input_candidates.push_back(fib); - } - std::cout.rdbuf (old); // <-- restore + std::vector< std::string > candidate_tract_files; + auto input_candidates = mitk::DiffusionDataIOHelper::load_fibs(candidate_tract_folders, &candidate_tract_files); + MITK_INFO << "Loaded " << candidate_tract_files.size() << " candidate tracts."; MITK_INFO << "Loaded " << reference_masks.size() << " reference masks."; MITK_INFO << "Loaded " << reference_peaks.size() << " reference peaks."; if (flipx || flipy || flipz) { itk::FlipPeaksFilter< float >::Pointer flipper = itk::FlipPeaksFilter< float >::New(); flipper->SetInput(peak_image); flipper->SetFlipX(flipx); flipper->SetFlipY(flipy); flipper->SetFlipZ(flipz); flipper->Update(); peak_image = flipper->GetOutput(); } + itk::ImageFileWriter< PeakImgType >::Pointer peak_image_writer = itk::ImageFileWriter< PeakImgType >::New(); + ofstream logfile; + logfile.open (out_folder + "scores.txt"); double rmse = 0.0; int iteration = 0; std::string name = "NOANCHOR"; if (parsedArgs.count("a")) { // Load reference tractogram consisting of all known tracts std::string anchors_file = us::any_cast(parsedArgs["a"]); mitk::FiberBundle::Pointer anchor_tractogram = mitk::IOUtil::Load(anchors_file); if ( !(anchor_tractogram.IsNull() || anchor_tractogram->GetNumFibers()==0) ) { // Fit known tracts to peak image to obtain underexplained image MITK_INFO << "Fit anchor tracts"; itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetTractograms({anchor_tractogram}); fitter->SetLambda(lambda); fitter->SetFilterOutliers(filter_outliers); fitter->SetPeakImage(peak_image); fitter->SetVerbose(true); fitter->SetMaskImage(mask); fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); rmse = fitter->GetRMSE(); vnl_vector rms_diff = fitter->GetRmsDiffPerBundle(); - logfile << "RMS_DIFF: " << setprecision(5) << rms_diff[0] << " " << name << " RMSE: " << rmse << "\n"; name = ist::GetFilenameWithoutExtension(anchors_file); mitk::FiberBundle::Pointer anchor_tracts = fitter->GetTractograms().at(0); anchor_tracts->SetFiberColors(255,255,255); mitk::IOUtil::Save(anchor_tracts, out_folder + boost::lexical_cast((int)(100000*rms_diff[0])) + "_" + name + ".fib"); + logfile << name << " " << setprecision(5) << rms_diff[0] << "\n"; + peak_image = fitter->GetUnderexplainedImage(); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + "Residual_" + name + ".nii.gz"); peak_image_writer->Update(); } } - - struct Overlap { - float v1; - float v2; - int i; - }; if (use_weights || use_num_streamlines) { MITK_INFO << "Using tract weights as scores"; int c = 0; for (auto fib : input_candidates) { int mod = 1; double score = 0; if (use_weights) { score = fib->GetFiberWeight(0); mod = 100000; } else if (use_num_streamlines) score = fib->GetNumFibers(); fib->ColorFibersByOrientation(); std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast((int)(mod*score)) + "_" + bundle_name + ".fib"); - - Overlap classic; classic.v1 = 0; classic.v2 = 0; classic.i = -1; - Overlap directional; directional.v1 = 0; directional.v2 = 0; directional.i = -1; - if(reference_masks.size()==reference_peaks.size()) - { - for (unsigned int i=0; iGetDirectionalOverlap(ref_mask, ref_peak); - if (directional_overlap>directional.v1) - { - directional.v1 = directional_overlap; - directional.v2 = overlap; - directional.i = i; - } - if (overlap>classic.v1) - { - classic.v1 = overlap; - classic.v2 = directional_overlap; - classic.i = i; - } - } - } - else - { - for (unsigned int i=0; iGetOverlap(ref_mask); - if (overlap>classic.v1) - { - classic.v1 = overlap; - classic.i = i; - } - } - } unsigned int num_voxels = 0; { itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New(); masks_filter->SetInputImage(mask); masks_filter->SetBinaryOutput(true); masks_filter->SetFiberBundle(fib); masks_filter->SetUseImageGeometry(true); masks_filter->Update(); num_voxels = masks_filter->GetNumCoveredVoxels(); } double weight_sum = 0; for (int i=0; iGetNumFibers(); i++) weight_sum += fib->GetFiberWeight(i); std::cout.rdbuf (old); // <-- restore - logfile << "RMS_DIFF: " << setprecision(5) << score << " " << bundle_name << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; - if (classic.i>=0) - { - logfile << "Best_overlap: " << setprecision(5) << classic.v1 << " " << classic.v2 << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.at(classic.i)) << "\n"; - if (directional.i>=0) - logfile << "Best_dir_overlap: " << setprecision(5) << directional.v1 << " " << directional.v2 << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.at(directional.i)) << "\n"; - } - else - logfile << "No_overlap\n"; + logfile << bundle_name << " " << setprecision(5) << score << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; ++c; } - } else if (!greedy_add) { MITK_INFO << "Fit candidate tracts"; itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetLambda(lambda); fitter->SetFilterOutliers(filter_outliers); fitter->SetVerbose(true); fitter->SetPeakImage(peak_image); fitter->SetMaskImage(mask); fitter->SetTractograms(input_candidates); fitter->SetFitIndividualFibers(true); if (regu=="MSM") fitter->SetRegularization(VnlCostFunction::REGU::MSM); else if (regu=="Variance") fitter->SetRegularization(VnlCostFunction::REGU::VARIANCE); else if (regu=="Lasso") fitter->SetRegularization(VnlCostFunction::REGU::LASSO); else if (regu=="VoxelVariance") fitter->SetRegularization(VnlCostFunction::REGU::VOXEL_VARIANCE); else if (regu=="GroupLasso") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_LASSO); else if (regu=="GroupVariance") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_VARIANCE); else if (regu=="NONE") fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); vnl_vector rms_diff = fitter->GetRmsDiffPerBundle(); int c = 0; for (auto fib : input_candidates) { std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect if (filter_zero_weights) fib = fib->FilterByWeights(0); mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast((int)(100000*rms_diff[c])) + "_" + bundle_name + ".fib"); - Overlap classic; classic.v1 = 0; classic.v2 = 0; classic.i = -1; - Overlap directional; directional.v1 = 0; directional.v2 = 0; directional.i = -1; - if(reference_masks.size()==reference_peaks.size()) - { - for (unsigned int i=0; iGetDirectionalOverlap(ref_mask, ref_peak); - if (directional_overlap>directional.v1) - { - directional.v1 = directional_overlap; - directional.v2 = overlap; - directional.i = i; - } - if (overlap>classic.v1) - { - classic.v1 = overlap; - classic.v2 = directional_overlap; - classic.i = i; - } - } - } - else - { - for (unsigned int i=0; iGetOverlap(ref_mask); - if (overlap>classic.v1) - { - classic.v1 = overlap; - classic.i = i; - } - } - } - unsigned int num_voxels = 0; { itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New(); masks_filter->SetInputImage(mask); masks_filter->SetBinaryOutput(true); masks_filter->SetFiberBundle(fib); masks_filter->SetUseImageGeometry(true); masks_filter->Update(); num_voxels = masks_filter->GetNumCoveredVoxels(); } double weight_sum = 0; for (int i=0; iGetNumFibers(); i++) weight_sum += fib->GetFiberWeight(i); std::cout.rdbuf (old); // <-- restore - logfile << "RMS_DIFF: " << setprecision(5) << rms_diff[c] << " " << bundle_name << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; - if (classic.i>=0) - { - logfile << "Best_overlap: " << setprecision(5) << classic.v1 << " " << classic.v2 << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.at(classic.i)) << "\n"; - if (directional.i>=0) - logfile << "Best_dir_overlap: " << setprecision(5) << directional.v1 << " " << directional.v2 << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.at(directional.i)) << "\n"; - } - else - logfile << "No_overlap\n"; - + logfile << bundle_name << " " << setprecision(5) << rms_diff[c] << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; ++c; } mitk::FiberBundle::Pointer out_fib = mitk::FiberBundle::New(); out_fib = out_fib->AddBundles(input_candidates); out_fib->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(out_fib, out_folder + "AllCandidates.fib"); peak_image = fitter->GetUnderexplainedImage(); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + "Residual_AllCandidates.nii.gz"); peak_image_writer->Update(); } else { MITK_INFO << "RMSE: " << setprecision(5) << rmse; // fitter->SetPeakImage(peak_image); // Iteratively add candidate bundles in a greedy manner while (!input_candidates.empty()) { double next_rmse = rmse; - double num_peaks = 0; mitk::FiberBundle::Pointer best_candidate = nullptr; PeakImgType::Pointer best_candidate_peak_image = nullptr; for (int i=0; i<(int)input_candidates.size(); ++i) { // WHY NECESSARY AGAIN?? itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetLambda(lambda); fitter->SetFilterOutliers(filter_outliers); fitter->SetVerbose(false); fitter->SetPeakImage(peak_image); fitter->SetMaskImage(mask); // ****************************** fitter->SetTractograms({input_candidates.at(i)}); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect fitter->Update(); std::cout.rdbuf (old); // <-- restore double candidate_rmse = fitter->GetRMSE(); if (candidate_rmseGetNumCoveredDirections(); best_candidate = fitter->GetTractograms().at(0); best_candidate_peak_image = fitter->GetUnderexplainedImage(); } } if (best_candidate.IsNull()) break; // fitter->SetPeakImage(peak_image); peak_image = best_candidate_peak_image; int i=0; std::vector< mitk::FiberBundle::Pointer > remaining_candidates; std::vector< std::string > remaining_candidate_files; for (auto fib : input_candidates) { if (fib!=best_candidate) { remaining_candidates.push_back(fib); remaining_candidate_files.push_back(candidate_tract_files.at(i)); } else name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(i)); ++i; } input_candidates = remaining_candidates; candidate_tract_files = remaining_candidate_files; iteration++; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect // Save winning candidate if (filter_zero_weights) best_candidate = best_candidate->FilterByWeights(0); mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast(iteration) + "_" + name + ".fib"); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + boost::lexical_cast(iteration) + "_" + name + ".nrrd"); peak_image_writer->Update(); - // Calculate best overlap with reference masks for evaluation purposes - float best_overlap = 0; - int best_overlap_index = -1; - i = 0; - for (auto ref_mask : reference_masks) - { - float overlap = best_candidate->GetOverlap(ref_mask); - if (overlap>best_overlap) - { - best_overlap = overlap; - best_overlap_index = i; - } - ++i; - } std::cout.rdbuf (old); // <-- restore - logfile << "RMSE: " << setprecision(5) << rmse << " " << name << " " << num_peaks << "\n"; - if (best_overlap_index>=0) - logfile << "Best_overlap: " << setprecision(5) << best_overlap << " " << ist::GetFilenameWithoutExtension(anchor_mask_files.at(best_overlap_index)) << "\n"; - else - logfile << "No_overlap\n"; +// logfile << name << " " << setprecision(5) << score << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; } } clock.Stop(); int h = clock.GetTotal()/3600; int m = ((int)clock.GetTotal()%3600)/60; int s = (int)clock.GetTotal()%60; MITK_INFO << "Plausibility estimation took " << h << "h, " << m << "m and " << s << "s"; logfile.close(); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CMakeLists.txt similarity index 76% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt rename to Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CMakeLists.txt index 2eb0e28099..6f3e0156a2 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CMakeLists.txt @@ -1,45 +1,48 @@ option(BUILD_DiffusionTractographyEvaluationCmdApps "Build commandline tools for diffusion fiber tractography evaluation" OFF) if(BUILD_DiffusionTractographyEvaluationCmdApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion cmdapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionTractographyEvaluationcmdapps - PeaksAngularError^^MitkFiberTracking - TractometerMetrics^^MitkFiberTracking - MergeOverlappingTracts^^MitkFiberTracking - ExtractSimilarTracts^^MitkFiberTracking - AnchorConstrainedPlausibility^^MitkFiberTracking - CalculateOverlap^^MitkFiberTracking - CheckEpsAndOverlap^^MitkFiberTracking + PeaksAngularError^^ + TractometerMetrics^^ + MergeOverlappingTracts^^ + GetOverlappingTracts^^ + ExtractSimilarTracts^^ + AnchorConstrainedPlausibility^^ + CalculateOverlap^^ + CheckEpsAndOverlap^^ + ReferenceSimilarity^^ + TractDistance^^ # LocalDirectionalFiberPlausibility^^MitkFiberTracking # HAS TO USE NEW PEAK IMAGE FORMAT ) foreach(diffusionTractographyEvaluationcmdapp ${diffusionTractographyEvaluationcmdapps}) # extract cmd app name and dependencies string(REPLACE "^^" "\\;" cmdapp_info ${diffusionTractographyEvaluationcmdapp}) set(cmdapp_info_list ${cmdapp_info}) list(GET cmdapp_info_list 0 appname) list(GET cmdapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitkFunctionCreateCommandLineApp( NAME ${appname} - DEPENDS MitkCore MitkDiffusionCore ${dependencies_list} - PACKAGE_DEPENDS ITK + DEPENDS MitkDiffusionCmdApps ${dependencies_list} + PACKAGE_DEPENDS ) endforeach() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CalculateOverlap.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CalculateOverlap.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CalculateOverlap.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CalculateOverlap.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CheckEpsAndOverlap.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CheckEpsAndOverlap.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/CheckEpsAndOverlap.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/CheckEpsAndOverlap.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ExtractSimilarTracts.cpp similarity index 81% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ExtractSimilarTracts.cpp index edf0a9c4a8..125265b2b3 100644 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/ExtractSimilarTracts.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ExtractSimilarTracts.cpp @@ -1,237 +1,196 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include -#include +#include -typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; -mitk::FiberBundle::Pointer LoadFib(std::string filename) -{ - std::vector fibInfile = mitk::IOUtil::Load(filename); - if( fibInfile.empty() ) - std::cout << "File " << filename << " could not be read!"; - mitk::BaseData::Pointer baseData = fibInfile.at(0); - return dynamic_cast(baseData.GetPointer()); -} - -ItkFloatImgType::Pointer LoadItkImage(const std::string& filename) -{ - mitk::Image::Pointer img = mitk::IOUtil::Load(filename); - ItkFloatImgType::Pointer itkMask = ItkFloatImgType::New(); - mitk::CastToItkImage(img, itkMask); - return itkMask; -} - /*! \brief Spatially cluster fibers */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Extract Similar Tracts"); parser.setCategory("Fiber Tracking Evaluation"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::InputFile, "Input:", "input fiber bundle (.fib, .trk, .tck)", us::Any(), false); parser.addArgument("ref_tracts", "", mitkCommandLineParser::StringList, "Ref. Tracts:", "reference tracts (.fib, .trk, .tck)", us::Any(), false); parser.addArgument("ref_masks", "", mitkCommandLineParser::StringList, "Ref. Masks:", "reference bundle masks", us::Any()); + parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); parser.addArgument("distance", "", mitkCommandLineParser::Int, "Distance:", "", 10); parser.addArgument("metric", "", mitkCommandLineParser::String, "Metric:", "EU_MEAN (default), EU_STD, EU_MAX"); parser.addArgument("subsample", "", mitkCommandLineParser::Float, "Subsampling factor:", "Only use specified fraction of input fibers", 1.0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string in_fib = us::any_cast(parsedArgs["i"]); std::string out_root = us::any_cast(parsedArgs["o"]); mitkCommandLineParser::StringContainerType ref_bundle_files = us::any_cast(parsedArgs["ref_tracts"]); mitkCommandLineParser::StringContainerType ref_mask_files; if (parsedArgs.count("ref_masks")) ref_mask_files = us::any_cast(parsedArgs["ref_masks"]); if (ref_mask_files.size()>0 && ref_mask_files.size()!=ref_bundle_files.size()) { MITK_INFO << "If reference masks are used, there has to be one mask per reference tract."; return EXIT_FAILURE; } int distance = 10; if (parsedArgs.count("distance")) distance = us::any_cast(parsedArgs["distance"]); std::string metric = "EU_MEAN"; if (parsedArgs.count("metric")) metric = us::any_cast(parsedArgs["metric"]); float subsample = 1.0; if (parsedArgs.count("subsample")) subsample = us::any_cast(parsedArgs["subsample"]); try { - mitk::FiberBundle::Pointer fib = LoadFib(in_fib); + mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(in_fib); std::srand(0); if (subsample<1.0) fib = fib->SubsampleFibers(subsample); mitk::FiberBundle::Pointer resampled_fib = fib->GetDeepCopy(); resampled_fib->ResampleToNumPoints(12); - std::vector< mitk::FiberBundle::Pointer > ref_fibs; - std::vector< ItkFloatImgType::Pointer > ref_masks; - for (std::size_t i=0; i(ref_mask_files); std::vector< float > distances; distances.push_back(distance); mitk::FiberBundle::Pointer extracted = mitk::FiberBundle::New(nullptr); unsigned int c = 0; for (auto ref_fib : ref_fibs) { MITK_INFO << "Extracting " << ist::GetFilenameName(ref_bundle_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect try { itk::TractClusteringFilter::Pointer segmenter = itk::TractClusteringFilter::New(); // calculate centroids from reference bundle { itk::TractClusteringFilter::Pointer clusterer = itk::TractClusteringFilter::New(); clusterer->SetDistances({10,20,30}); clusterer->SetTractogram(ref_fib); clusterer->SetMetrics({new mitk::ClusteringMetricEuclideanStd()}); clusterer->SetMergeDuplicateThreshold(0.0); clusterer->Update(); std::vector tracts = clusterer->GetOutCentroids(); ref_fib = mitk::FiberBundle::New(nullptr); ref_fib = ref_fib->AddBundles(tracts); mitk::IOUtil::Save(ref_fib, out_root + "centroids_" + ist::GetFilenameName(ref_bundle_files.at(c))); segmenter->SetInCentroids(ref_fib); } // segment tract segmenter->SetFilterMask(ref_masks.at(c)); segmenter->SetOverlapThreshold(0.8); segmenter->SetDistances(distances); segmenter->SetTractogram(resampled_fib); segmenter->SetMergeDuplicateThreshold(0.0); segmenter->SetDoResampling(false); if (metric=="EU_MEAN") segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanMean()}); else if (metric=="EU_STD") segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanStd()}); else if (metric=="EU_MAX") segmenter->SetMetrics({new mitk::ClusteringMetricEuclideanMax()}); segmenter->Update(); std::vector< std::vector< long > > clusters = segmenter->GetOutFiberIndices(); if (clusters.size()>0) { vtkSmartPointer weights = vtkSmartPointer::New(); mitk::FiberBundle::Pointer result = mitk::FiberBundle::New(nullptr); std::vector< mitk::FiberBundle::Pointer > result_fibs; for (unsigned int cluster_index=0; cluster_indexGeneratePolyDataByIds(clusters.at(cluster_index), weights))); result = result->AddBundles(result_fibs); extracted = extracted->AddBundle(result); mitk::IOUtil::Save(result, out_root + "extracted_" + ist::GetFilenameName(ref_bundle_files.at(c))); fib = mitk::FiberBundle::New(fib->GeneratePolyDataByIds(clusters.back(), weights)); resampled_fib = mitk::FiberBundle::New(resampled_fib->GeneratePolyDataByIds(clusters.back(), weights)); } } catch(itk::ExceptionObject& excpt) { MITK_INFO << "Exception while processing " << ist::GetFilenameName(ref_bundle_files.at(c)); MITK_INFO << excpt.GetDescription(); } catch(std::exception& excpt) { MITK_INFO << "Exception while processing " << ist::GetFilenameName(ref_bundle_files.at(c)); MITK_INFO << excpt.what(); } std::cout.rdbuf (old); // <-- restore if (fib->GetNumFibers()==0) break; ++c; } MITK_INFO << "Extracted streamlines: " << extracted->GetNumFibers(); mitk::IOUtil::Save(extracted, out_root + "extracted_streamlines.trk"); MITK_INFO << "Residual streamlines: " << fib->GetNumFibers(); mitk::IOUtil::Save(fib, out_root + "residual_streamlines.trk"); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/GetOverlappingTracts.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/GetOverlappingTracts.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/LocalDirectionalFiberPlausibility.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/LocalDirectionalFiberPlausibility.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/LocalDirectionalFiberPlausibility.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/LocalDirectionalFiberPlausibility.cpp diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeOverlappingTracts.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/MergeOverlappingTracts.cpp similarity index 85% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeOverlappingTracts.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/MergeOverlappingTracts.cpp index a6afab75a0..f2b12676e4 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/MergeOverlappingTracts.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/MergeOverlappingTracts.cpp @@ -1,249 +1,215 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include -#include -#include #include #include #include #include #include #include #include #include #include +#include -typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUIntImgType; -std::vector< std::string > get_file_list(const std::string& path, std::vector< std::string > extensions={".fib", ".trk"}) -{ - std::vector< std::string > file_list; - itk::Directory::Pointer dir = itk::Directory::New(); - - if (dir->Load(path.c_str())) - { - int n = dir->GetNumberOfFiles(); - for (int r = 0; r < n; r++) - { - const char *filename = dir->GetFile(r); - std::string ext = ist::GetFilenameExtension(filename); - for (auto e : extensions) - { - if (ext==e) - { - file_list.push_back(path + '/' + filename); - break; - } - } - } - } - return file_list; -} /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Merge Overlapping Tracts"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::InputFile, "Input Folder:", "input folder", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output Folder:", "output folder", us::Any(), false); parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap threshold:", "Tracts with overlap larger than this threshold are merged", 0.8, false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string input_folder = us::any_cast(parsedArgs["i"]); std::string out_folder = us::any_cast(parsedArgs["o"]); float overlap = 0.8; if (parsedArgs.count("overlap")) overlap = us::any_cast(parsedArgs["overlap"]); try { if (!ist::PathExists(out_folder)) ist::MakeDirectory(out_folder); - std::vector< std::string > fib_files = get_file_list(input_folder, {".fib", ".trk", ".tck"}); - if (fib_files.empty()) - return EXIT_FAILURE; + std::vector< mitk::FiberBundle::Pointer > fibs = mitk::DiffusionDataIOHelper::load_fibs({input_folder}); + std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect - std::vector< mitk::FiberBundle::Pointer > fibs; - for (std::string f : fib_files) - { - mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(f); - fibs.push_back(fib); - } - mitk::FiberBundle::Pointer combined = mitk::FiberBundle::New(); combined = combined->AddBundles(fibs); itk::TractsToFiberEndingsImageFilter< ItkFloatImgType >::Pointer endings = itk::TractsToFiberEndingsImageFilter< ItkFloatImgType >::New(); endings->SetFiberBundle(combined); endings->SetUpsamplingFactor(0.25); endings->Update(); ItkFloatImgType::Pointer ref_image = endings->GetOutput(); std::cout.rdbuf (old); // <-- restore for (int its = 0; its<3; its++) { std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect std::vector< ItkFloatImgType::Pointer > mask_images; for (auto fib : fibs) { itk::TractDensityImageFilter< ItkFloatImgType >::Pointer masks = itk::TractDensityImageFilter< ItkFloatImgType >::New(); masks->SetInputImage(ref_image); masks->SetBinaryOutput(true); masks->SetFiberBundle(fib); masks->SetUseImageGeometry(true); masks->Update(); mask_images.push_back(masks->GetOutput()); } int r=0; vnl_matrix< int > mat; mat.set_size(mask_images.size(), mask_images.size()); mat.fill(0); for (auto m1 : mask_images) { float max_overlap = overlap; int c = 0; for (auto m2 : mask_images) { if (c<=r) { ++c; continue; } itk::ImageRegionConstIterator it1(m1, m1->GetLargestPossibleRegion()); itk::ImageRegionConstIterator it2(m2, m2->GetLargestPossibleRegion()); unsigned int c1 = 0; unsigned int c2 = 0; unsigned int intersect = 0; while( !it1.IsAtEnd() ) { if( it1.Get()>0 && it2.Get()>0) ++intersect; if(it1.Get()>0) ++c1; if(it2.Get()>0) ++c2; ++it1; ++it2; } if ( (float)intersect/c1>max_overlap ) { max_overlap = (float)intersect/c1; mat.put(r,c, 1); } if ( (float)intersect/c2>max_overlap ) { max_overlap = (float)intersect/c2; mat.put(r,c, 1); } ++c; } ++r; } std::vector< mitk::FiberBundle::Pointer > out_fibs; std::vector< bool > used; for (unsigned int i=0; i0) { fib = fib->AddBundle(fibs.at(c)); used[c] = true; } } out_fibs.push_back(fib); } std::cout.rdbuf (old); // <-- restore MITK_INFO << fibs.size() << " --> " << out_fibs.size(); if (fibs.size()==out_fibs.size()) break; fibs = out_fibs; } int c = 0; for (auto fib : fibs) { std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::IOUtil::Save(fib, out_folder + "/bundle_" + boost::lexical_cast(c) + ".trk"); std::cout.rdbuf (old); // <-- restore ++c; } } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/PeaksAngularError.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/PeaksAngularError.cpp similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/PeaksAngularError.cpp rename to Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/PeaksAngularError.cpp diff --git a/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp new file mode 100644 index 0000000000..28f376323b --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp @@ -0,0 +1,149 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +typedef itk::Image< unsigned char, 3 > ItkUcharImageType; +typedef itk::Image< float, 4 > ItkPeakImgType; + +int main(int argc, char* argv[]) +{ + mitkCommandLineParser parser; + + parser.setTitle("Reference Similarity"); + parser.setCategory("Fiber Tracking Evaluation"); + parser.setContributor("MIC"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("", "i", mitkCommandLineParser::StringList, "Input Tracts:", "input tracts folder", us::Any(), false); + parser.addArgument("reference_tracts", "", mitkCommandLineParser::StringList, "", "", us::Any(), false); + parser.addArgument("reference_masks", "", mitkCommandLineParser::StringList, "", "", us::Any(), false); + parser.addArgument("reference_peaks", "", mitkCommandLineParser::StringList, "", "", us::Any(), false); + + parser.addArgument("", "o", mitkCommandLineParser::String, "", "", us::Any(), false); + + parser.addArgument("fiber_points", "", mitkCommandLineParser::Int, "Fiber points:", "", 20); + + std::map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + std::string out_folder = us::any_cast(parsedArgs["o"]); + mitkCommandLineParser::StringContainerType input_tract_files = us::any_cast(parsedArgs["i"]); + mitkCommandLineParser::StringContainerType reference_tract_files = us::any_cast(parsedArgs["reference_tracts"]); + mitkCommandLineParser::StringContainerType reference_mask_files = us::any_cast(parsedArgs["reference_masks"]); + mitkCommandLineParser::StringContainerType reference_peak_files = us::any_cast(parsedArgs["reference_peaks"]); + + int fiber_points = 20; + if (parsedArgs.count("fiber_points")) + fiber_points = us::any_cast(parsedArgs["fiber_points"]); + + try + { + std::vector input_tract_names; + std::vector ref_tract_names; + std::vector< mitk::FiberBundle::Pointer > input_tracts = mitk::DiffusionDataIOHelper::load_fibs(input_tract_files, &input_tract_names); + std::vector< mitk::FiberBundle::Pointer > reference_tracts = mitk::DiffusionDataIOHelper::load_fibs(reference_tract_files, &ref_tract_names); + std::vector< ItkUcharImageType::Pointer > reference_masks = mitk::DiffusionDataIOHelper::load_itk_images(reference_mask_files); + std::vector< ItkPeakImgType::Pointer > reference_peaks = mitk::DiffusionDataIOHelper::load_itk_images(reference_peak_files); + + MITK_INFO << "Calculating distances"; + itk::TractDistanceFilter::Pointer distance_calculator = itk::TractDistanceFilter::New(); + distance_calculator->SetNumPoints(fiber_points); + distance_calculator->SetTracts1(input_tracts); + distance_calculator->SetTracts2(reference_tracts); + distance_calculator->SetMetrics({new mitk::ClusteringMetricEuclideanMean()}); + distance_calculator->Update(); + auto distances = distance_calculator->GetAllDistances(); + + vnl_matrix voxel_overlap; voxel_overlap.set_size(input_tracts.size(), reference_tracts.size()); + vnl_matrix dir_overlap; dir_overlap.set_size(input_tracts.size(), reference_tracts.size()); + + MITK_INFO << "Calculating overlap"; + boost::progress_display disp(input_tracts.size()*reference_tracts.size()); + int r=0; + for (auto fib : input_tracts) + { + int c=0; + for (auto ref_mask : reference_masks) + { + ++disp; + std::streambuf *old = cout.rdbuf(); // <-- save + std::stringstream ss; + std::cout.rdbuf (ss.rdbuf()); // <-- redirect + float overlap = 0; + float directional_overlap = 0; + std::tie(directional_overlap, overlap) = fib->GetDirectionalOverlap(ref_mask, reference_peaks.at(c)); + voxel_overlap[r][c] = overlap; + dir_overlap[r][c] = directional_overlap; + std::cout.rdbuf (old); // <-- restore + ++c; + } + ++r; + } + + ofstream logfile; + logfile.open(out_folder + "ref_tract_names.txt"); + for (unsigned int i=0; i #include #include #include #include #include #include #include -#include -#include +#include -typedef itksys::SystemTools ist; - -std::vector< std::string > get_file_list(const std::string& path, std::vector< std::string > extensions={".fib", ".trk"}) -{ - std::vector< std::string > file_list; - itk::Directory::Pointer dir = itk::Directory::New(); - - if (dir->Load(path.c_str())) - { - int n = dir->GetNumberOfFiles(); - for (int r = 0; r < n; r++) - { - const char *filename = dir->GetFile(r); - std::string ext = ist::GetFilenameExtension(filename); - for (auto e : extensions) - { - if (ext==e) - { - file_list.push_back(path + '/' + filename); - break; - } - } - } - } - std::sort(file_list.begin(), file_list.end()); - return file_list; -} int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Tract Distance"); parser.setCategory("Fiber Processing"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i1", mitkCommandLineParser::String, "Input Folder 1:", "input tracts folder 1", us::Any(), false); parser.addArgument("", "i2", mitkCommandLineParser::String, "Input Folder 2:", "input tracts folder 2", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output logfile", us::Any(), false); parser.addArgument("fiber_points", "", mitkCommandLineParser::Int, "Fiber points:", "", 12); parser.addArgument("metrics", "", mitkCommandLineParser::StringList, "Metrics:", "EU_MEAN (default), EU_STD, EU_MAX"); parser.addArgument("metric_weights", "", mitkCommandLineParser::StringList, "Metric weights:", "add one float weight for each used metric"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string t1_folder = us::any_cast(parsedArgs["i1"]); std::string t2_folder = us::any_cast(parsedArgs["i2"]); std::string out_file = us::any_cast(parsedArgs["o"]); int fiber_points = 12; if (parsedArgs.count("fiber_points")) fiber_points = us::any_cast(parsedArgs["fiber_points"]); std::vector< std::string > metric_strings = {"EU_MEAN"}; if (parsedArgs.count("metrics")) metric_strings = us::any_cast(parsedArgs["metrics"]); std::vector< std::string > metric_weights = {"1.0"}; if (parsedArgs.count("metric_weights")) metric_weights = us::any_cast(parsedArgs["metric_weights"]); if (metric_strings.size()!=metric_weights.size()) { MITK_INFO << "Each metric needs an associated metric weight!"; return EXIT_FAILURE; } try { - std::streambuf *old = cout.rdbuf(); // <-- save - std::stringstream ss; - std::cout.rdbuf (ss.rdbuf()); // <-- redirect - std::vector< mitk::FiberBundle::Pointer > tractograms1; - std::vector< mitk::FiberBundle::Pointer > tractograms2; - - auto t1_files = get_file_list(t1_folder, {".fib",".trk",".tck"}); - for (auto f : t1_files) - { - mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(f); - tractograms1.push_back(fib); - } - auto t2_files = get_file_list(t2_folder, {".fib",".trk",".tck"}); - for (auto f : t2_files) - { - mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(f); - tractograms2.push_back(fib); - } + std::vector t1_files; + std::vector< mitk::FiberBundle::Pointer > tractograms1 = mitk::DiffusionDataIOHelper::load_fibs({t1_folder}, &t1_files); + std::vector t2_files; + std::vector< mitk::FiberBundle::Pointer > tractograms2 = mitk::DiffusionDataIOHelper::load_fibs({t2_folder}, &t2_files); - std::cout.rdbuf (old); // <-- restore MITK_INFO << "Loaded " << tractograms1.size() << " source tractograms."; MITK_INFO << "Loaded " << tractograms2.size() << " target tractograms."; itk::TractDistanceFilter::Pointer distance_calculator = itk::TractDistanceFilter::New(); distance_calculator->SetNumPoints(fiber_points); distance_calculator->SetTracts1(tractograms1); distance_calculator->SetTracts2(tractograms2); std::vector< mitk::ClusteringMetric* > metrics; int mc = 0; for (auto m : metric_strings) { float w = boost::lexical_cast(metric_weights.at(mc)); MITK_INFO << "Metric: " << m << " (w=" << w << ")"; if (m=="EU_MEAN") metrics.push_back({new mitk::ClusteringMetricEuclideanMean()}); else if (m=="EU_STD") metrics.push_back({new mitk::ClusteringMetricEuclideanStd()}); else if (m=="EU_MAX") metrics.push_back({new mitk::ClusteringMetricEuclideanMax()}); metrics.back()->SetScale(w); mc++; } if (metrics.empty()) { MITK_INFO << "No metric selected!"; return EXIT_FAILURE; } distance_calculator->SetMetrics(metrics); distance_calculator->Update(); MITK_INFO << "Distances:"; - auto distances = distance_calculator->GetDistances(); - auto indices = distance_calculator->GetIndices(); + auto distances = distance_calculator->GetMinDistances(); + auto indices = distance_calculator->GetMinIndices(); ofstream logfile; logfile.open (out_file); for (unsigned int i=0; i mitk::DiffusionDataIOHelper::get_file_list(const std::string& path, const std::vector extensions) +{ + std::vector< std::string > file_list; + itk::Directory::Pointer dir = itk::Directory::New(); + + if (dir->Load(path.c_str())) + { + int n = dir->GetNumberOfFiles(); + for (int r = 0; r < n; r++) + { + const char *filename = dir->GetFile(r); + std::string ext = ist::GetFilenameExtension(filename); + for (auto e : extensions) + { + if (ext==e) + { + file_list.push_back(path + '/' + filename); + break; + } + } + } + } + std::sort(file_list.begin(), file_list.end()); + return file_list; +} + +std::vector< mitk::FiberBundle::Pointer > mitk::DiffusionDataIOHelper::load_fibs(const std::vector files, std::vector* filenames) +{ + std::streambuf *old = cout.rdbuf(); // <-- save + std::stringstream ss; + std::cout.rdbuf (ss.rdbuf()); // <-- redirect + std::vector< mitk::FiberBundle::Pointer > out; + for (auto f : files) + { + if (itksys::SystemTools::PathExists(f)) + { + if (!f.empty() && f.back() != '/') + f += "/"; + + auto list = get_file_list(f, {".fib",".trk",".tck"}); + for (auto file : list) + { + mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(file); + out.push_back(fib); + if (filenames!=nullptr) + filenames->push_back(file); + } + } + else if (itksys::SystemTools::FileExists(f)) + { + mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(f); + out.push_back(fib); + if (filenames!=nullptr) + filenames->push_back(f); + } + } + std::cout.rdbuf (old); // <-- restore + MITK_INFO << "Loaded " << out.size() << " tractograms"; + return out; +} diff --git a/Modules/DiffusionImaging/DiffusionCmdApps/mitkDiffusionDataIOHelper.h b/Modules/DiffusionImaging/DiffusionCmdApps/mitkDiffusionDataIOHelper.h new file mode 100644 index 0000000000..4923527fc4 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCmdApps/mitkDiffusionDataIOHelper.h @@ -0,0 +1,116 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef __mitkDiffusionDataIOHelper_h_ +#define __mitkDiffusionDataIOHelper_h_ + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +typedef itksys::SystemTools ist; + +namespace mitk{ + +class MITKDIFFUSIONCMDAPPS_EXPORT DiffusionDataIOHelper +{ +public: + + static std::vector< std::string > get_file_list(const std::string& path, const std::vector< std::string > extensions={".fib", ".trk"}); + + static std::vector< mitk::FiberBundle::Pointer > load_fibs(const std::vector files, std::vector* filenames=nullptr); + + template< class TYPE > + static typename TYPE::Pointer load_itk_image(const std::string file) + { + if (file.compare("")!=0 && itksys::SystemTools::FileExists(file)) + { + mitk::Image::Pointer image = mitk::IOUtil::Load(file); + + typedef mitk::ImageToItk< TYPE > CasterType; + typename CasterType::Pointer caster = CasterType::New(); + caster->SetInput(image); + caster->Update(); + typename TYPE::Pointer itk_image = caster->GetOutput(); + return itk_image; + } + return nullptr; + } + + template< class TYPE > + static std::vector< typename TYPE::Pointer > load_itk_images(const std::vector files, std::vector* filenames=nullptr) + { + std::streambuf *old = cout.rdbuf(); // <-- save + std::stringstream ss; + std::cout.rdbuf (ss.rdbuf()); // <-- redirect + std::vector< typename TYPE::Pointer > out; + for (auto f : files) + { + if (itksys::SystemTools::PathExists(f)) + { + if (!f.empty() && f.back() != '/') + f += "/"; + + auto list = get_file_list(f, {".nrrd",".nii.gz",".nii"}); + for (auto file : list) + { + mitk::Image::Pointer image = mitk::IOUtil::Load(file); + + typedef mitk::ImageToItk< TYPE > CasterType; + typename CasterType::Pointer caster = CasterType::New(); + caster->SetInput(image); + caster->Update(); + typename TYPE::Pointer itk_image = caster->GetOutput(); + + out.push_back(itk_image); + if (filenames!=nullptr) + filenames->push_back(file); + } + } + else if (itksys::SystemTools::FileExists(f)) + { + mitk::Image::Pointer image = mitk::IOUtil::Load(f); + + typedef mitk::ImageToItk< TYPE > CasterType; + typename CasterType::Pointer caster = CasterType::New(); + caster->SetInput(image); + caster->Update(); + typename TYPE::Pointer itk_image = caster->GetOutput(); + + out.push_back(itk_image); + if (filenames!=nullptr) + filenames->push_back(f); + } + } + + std::cout.rdbuf (old); // <-- restore + MITK_INFO << "Loaded " << out.size() << " images"; + return out; + } + +}; + +} + +#endif //__mitkDiffusionDataIOHelper_h_ + diff --git a/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt index 1b9cedf87a..c2b6ac476b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt @@ -1,26 +1,24 @@ # With apple gcc 4.2.1 the following waring leads to an build error if boost is enabled if(APPLE) mitkFunctionCheckCAndCXXCompilerFlags("-Wno-error=empty-body" CMAKE_C_FLAGS CMAKE_CXX_FLAGS) endif() MITK_CREATE_MODULE( SUBPROJECTS MITK-DTI INCLUDE_DIRS include/ include/Algorithms include/Algorithms/Reconstruction include/Algorithms/Registration include/Algorithms/Reconstruction/MultishellProcessing include/Algorithms/Reconstruction/FittingFunctions include/DicomImport include/IODataStructures/DiffusionWeightedImages include/IODataStructures/Properties include/IODataStructures/OdfImages include/IODataStructures/TensorImages include/IODataStructures include/Rendering ${CMAKE_CURRENT_BINARY_DIR} DEPENDS MitkMapperExt MitkPlanarFigure MitkImageExtraction MitkSceneSerializationBase MitkDICOMReader MitkMatchPointRegistration PACKAGE_DEPENDS PUBLIC ITK|ITKTestKernel+ITKRegistrationCommon+ITKMetricsv4+ITKRegistrationMethodsv4+ITKDistanceMap+ITKLabelVoting+ITKVTK PUBLIC VTK|vtkFiltersProgrammable ) if(MSVC) mitkFunctionCheckCAndCXXCompilerFlags("/wd4005" CMAKE_C_FLAGS CMAKE_CXX_FLAGS) endif() add_subdirectory(Testing) -add_subdirectory(cmdapps) add_subdirectory(autoload/IO) if(MITK_USE_Python) MITK_INSTALL(FILES PythonRequirements.txt) -add_subdirectory(cmdapps_python) endif() diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt deleted file mode 100644 index 19e45bfb5f..0000000000 --- a/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt +++ /dev/null @@ -1,44 +0,0 @@ -option(BUILD_DiffusionCoreCmdApps "Build commandline tools for diffusion" OFF) - -if(BUILD_DiffusionCoreCmdApps OR MITK_BUILD_ALL_APPS) - - # needed include directories - include_directories( - ${CMAKE_CURRENT_SOURCE_DIR} - ${CMAKE_CURRENT_BINARY_DIR} - ) - - # list of diffusion cmdapps - # if an app requires additional dependencies - # they are added after a "^^" and separated by "_" - set( diffusioncorecmdapps - ImageResampler^^ - CopyGeometry^^ - Registration^^ - DiffusionDICOMLoader^^ - ResampleGradients^^ - ShToOdfImage^^ - DImp^^ - DReg^^ - HeadMotionCorrection^^ - DmriDenoising^^ - RoundBvalues^^ - ) - - foreach(diffusioncorecmdapp ${diffusioncorecmdapps}) - # extract cmd app name and dependencies - string(REPLACE "^^" "\\;" cmdapp_info ${diffusioncorecmdapp}) - set(cmdapp_info_list ${cmdapp_info}) - list(GET cmdapp_info_list 0 appname) - list(GET cmdapp_info_list 1 raw_dependencies) - string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") - set(dependencies_list ${dependencies}) - - mitkFunctionCreateCommandLineApp( - NAME ${appname} - DEPENDS MitkCore MitkDiffusionCore MitkImageDenoising ${dependencies_list} - PACKAGE_DEPENDS ITK - ) - endforeach() - -endif() diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp index 767c342fe2..18bb94c24b 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberExtractionFilter.cpp @@ -1,502 +1,504 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkFiberExtractionFilter_cpp #define __itkFiberExtractionFilter_cpp #include "itkFiberExtractionFilter.h" #define _USE_MATH_DEFINES #include #include #include #include namespace itk{ template< class PixelType > FiberExtractionFilter< PixelType >::FiberExtractionFilter() : m_DontResampleFibers(false) , m_Mode(MODE::OVERLAP) , m_InputType(INPUT::SCALAR_MAP) , m_BothEnds(true) , m_OverlapFraction(0.8) , m_NoNegatives(false) , m_NoPositives(false) , m_Interpolate(false) , m_Threshold(0.5) - , m_Labels({1}) + , m_Labels() , m_SkipSelfConnections(false) , m_OnlySelfConnections(false) , m_SplitByRoi(false) , m_SplitLabels(false) , m_MinFibersPerTract(0) , m_PairedStartEndLabels(false) { m_Interpolator = itk::LinearInterpolateImageFunction< itk::Image< PixelType, 3 >, float >::New(); } template< class PixelType > FiberExtractionFilter< PixelType >::~FiberExtractionFilter() { } template< class PixelType > void FiberExtractionFilter< PixelType >::SetRoiImages(const std::vector &rois) { for (auto roi : rois) { if (roi==nullptr) { MITK_INFO << "ROI image is NULL!"; return; } } m_RoiImages = rois; } template< class PixelType > void FiberExtractionFilter< PixelType >::SetRoiImageNames(const std::vector< std::string > roi_names) { m_RoiImageNames = roi_names; } template< class PixelType > mitk::FiberBundle::Pointer FiberExtractionFilter< PixelType >::CreateFib(std::vector< long >& ids) { vtkSmartPointer weights = vtkSmartPointer::New(); vtkSmartPointer pTmp = m_InputFiberBundle->GeneratePolyDataByIds(ids, weights); mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); fib->SetFiberWeights(weights); return fib; } template< class PixelType > bool FiberExtractionFilter< PixelType >::IsPositive(const itk::Point& itkP) { if( m_InputType == INPUT::SCALAR_MAP ) return mitk::imv::IsInsideMask(itkP, m_Interpolate, m_Interpolator, m_Threshold); else if( m_InputType == INPUT::LABEL_MAP ) { auto val = mitk::imv::GetImageValue(itkP, m_Interpolate, m_Interpolator); for (auto l : m_Labels) if (l==val) return true; } else mitkThrow() << "No valid input type selected!"; return false; } template< class PixelType > std::vector< std::string > FiberExtractionFilter< PixelType >::GetPositiveLabels() const { return m_PositiveLabels; } template< class PixelType > void FiberExtractionFilter< PixelType >::ExtractOverlap(mitk::FiberBundle::Pointer fib) { MITK_INFO << "Extracting fibers (min. overlap " << m_OverlapFraction << ")"; vtkSmartPointer polydata = fib->GetFiberPolyData(); std::vector< std::vector< long > > positive_ids; // one ID vector per ROI positive_ids.resize(m_RoiImages.size()); std::vector< long > negative_ids; // fibers not overlapping with ANY mask boost::progress_display disp(m_InputFiberBundle->GetNumFibers()); for (int i=0; iGetNumFibers(); i++) { ++disp; vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float best_ol = 0; int best_ol_idx = -1; for (unsigned int m=0; mSetInputImage(roi); PixelType inside = 0; PixelType outside = 0; for (int j=0; j startVertex = mitk::imv::GetItkPoint(points->GetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; roi->TransformPhysicalPointToIndex(startVertex, startIndex); roi->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); itk::Point endVertex = mitk::imv::GetItkPoint(points->GetPoint(j + 1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; roi->TransformPhysicalPointToIndex(endVertex, endIndex); roi->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(roi->GetSpacing(), startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > segment : segments) { if (!roi->GetLargestPossibleRegion().IsInside(segment.first)) continue; if (roi->GetPixel(segment.first)>=m_Threshold) inside += segment.second; else outside += segment.second; } } float overlap = (float)inside/(inside+outside); if (overlap > best_ol && overlap >= m_OverlapFraction) { best_ol_idx = m; best_ol = overlap; } } if (best_ol_idx<0) negative_ids.push_back(i); else positive_ids.at(best_ol_idx).push_back(i); } if (!m_NoNegatives) m_Negatives.push_back(CreateFib(negative_ids)); if (!m_NoPositives) { for (unsigned int i=0; i(i); if (positive_ids.at(i).size()>=m_MinFibersPerTract) { m_Positives.push_back(CreateFib(positive_ids.at(i))); m_PositiveLabels.push_back(name); } } if (!m_SplitByRoi) { mitk::FiberBundle::Pointer output = mitk::FiberBundle::New(nullptr); output = output->AddBundles(m_Positives); m_Positives.clear(); m_Positives.push_back(output); m_PositiveLabels.clear(); m_PositiveLabels.push_back(""); } } } template< class PixelType > void FiberExtractionFilter< PixelType >::ExtractEndpoints(mitk::FiberBundle::Pointer fib) { MITK_INFO << "Extracting fibers (endpoints in mask)"; vtkSmartPointer polydata = fib->GetFiberPolyData(); std::vector< std::vector< long > > positive_ids; // one ID vector per ROI positive_ids.resize(m_RoiImages.size()); std::vector< long > negative_ids; // fibers not overlapping with ANY mask boost::progress_display disp(m_InputFiberBundle->GetNumFibers()); for (int i=0; iGetNumFibers(); i++) { ++disp; vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); bool positive = false; if (numPoints>1) for (unsigned int m=0; mSetInputImage(roi); int inside = 0; // check first fiber point { double* p = points->GetPoint(0); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; if ( IsPositive(itkP) ) inside++; } // check second fiber point { double* p = points->GetPoint(numPoints-1); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; if ( IsPositive(itkP) ) inside++; } if (inside==2 || (inside==1 && !m_BothEnds)) { positive = true; positive_ids[m].push_back(i); } } if (!positive) negative_ids.push_back(i); } if (!m_NoNegatives) m_Negatives.push_back(CreateFib(negative_ids)); if (!m_NoPositives) { for (unsigned int i=0; i(i); if (positive_ids.at(i).size()>=m_MinFibersPerTract) { m_Positives.push_back(CreateFib(positive_ids.at(i))); m_PositiveLabels.push_back(name); } } if (!m_SplitByRoi) { mitk::FiberBundle::Pointer output = mitk::FiberBundle::New(nullptr); output = output->AddBundles(m_Positives); m_Positives.clear(); m_Positives.push_back(output); m_PositiveLabels.clear(); m_PositiveLabels.push_back(""); } } } template< class PixelType > void FiberExtractionFilter< PixelType >::ExtractLabels(mitk::FiberBundle::Pointer fib) { MITK_INFO << "Extracting fibers by labels"; vtkSmartPointer polydata = fib->GetFiberPolyData(); std::map< std::string, std::vector< long > > positive_ids; std::vector< long > negative_ids; // fibers not overlapping with ANY label boost::progress_display disp(m_InputFiberBundle->GetNumFibers()); for (int i=0; iGetNumFibers(); i++) { ++disp; vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); bool positive = false; if (numPoints>1) for (unsigned int m=0; mSetInputImage(roi); int inside = 0; double* p1 = points->GetPoint(0); itk::Point itkP1; itkP1[0] = p1[0]; itkP1[1] = p1[1]; itkP1[2] = p1[2]; short label1 = mitk::imv::GetImageValue(itkP1, m_Interpolate, m_Interpolator); double* p2 = points->GetPoint(numPoints-1); itk::Point itkP2; itkP2[0] = p2[0]; itkP2[1] = p2[1]; itkP2[2] = p2[2]; short label2 = mitk::imv::GetImageValue(itkP2, m_Interpolate, m_Interpolator); if (!m_Labels.empty()) // extract fibers from all pairwise label combinations { for (auto l : m_Labels) { if (l==label1) inside++; if (l==label2) inside++; if (inside==2) break; } } - else // extract fibers between start and end labels + else if (!m_StartLabels.empty() || !m_EndLabels.empty()) // extract fibers between start and end labels { - m_BothEnds = true; // if we have start and end labels it does not make sense to not use both endpoints + if (!m_StartLabels.empty() && !m_EndLabels.empty()) + m_BothEnds = true; // if we have start and end labels it does not make sense to not use both endpoints if (m_PairedStartEndLabels) { if (m_StartLabels.size()!=m_EndLabels.size()) mitkThrow() << "Start and end label lists must have same size if paired labels are used"; for (unsigned int ii=0; ii(label1) + "-" + boost::lexical_cast(label2); else key = boost::lexical_cast(label2) + "-" + boost::lexical_cast(label1); } if (m_SplitByRoi) { if (m(m) + "_" + key; } if (m_BothEnds) { if ( (inside==2 && (!m_SkipSelfConnections || label1!=label2)) || (inside==2 && m_OnlySelfConnections && label1==label2) ) { positive = true; if ( positive_ids.count(key)==0 ) positive_ids.insert( std::pair< std::string, std::vector< long > >( key, {i}) ); else positive_ids[key].push_back(i); } } else { if ( (inside>=1 && (!m_SkipSelfConnections || label1!=label2)) || (inside==2 && m_OnlySelfConnections && label1==label2) ) { positive = true; if ( positive_ids.count(key)==0 ) positive_ids.insert( std::pair< std::string, std::vector< long > >( key, {i}) ); else positive_ids[key].push_back(i); } } } if (!positive) negative_ids.push_back(i); } if (!m_NoNegatives) m_Negatives.push_back(CreateFib(negative_ids)); if (!m_NoPositives) { for (auto label : positive_ids) { if (label.second.size()>=m_MinFibersPerTract) { m_Positives.push_back(CreateFib(label.second)); m_PositiveLabels.push_back(label.first); } } - -// if (!m_SplitLabels) -// { -// mitk::FiberBundle::Pointer output = mitk::FiberBundle::New(nullptr); -// output = output->AddBundles(m_Positives); -// m_Positives.clear(); -// m_Positives.push_back(output); -// m_PositiveLabels.clear(); -// } } } template< class PixelType > void FiberExtractionFilter< PixelType >::SetLabels(const std::vector &Labels) { m_Labels = Labels; } template< class PixelType > void FiberExtractionFilter< PixelType >::SetStartLabels(const std::vector &Labels) { m_StartLabels = Labels; } template< class PixelType > void FiberExtractionFilter< PixelType >::SetEndLabels(const std::vector &Labels) { m_EndLabels = Labels; } template< class PixelType > std::vector FiberExtractionFilter< PixelType >::GetNegatives() const { return m_Negatives; } template< class PixelType > std::vector FiberExtractionFilter< PixelType >::GetPositives() const { return m_Positives; } template< class PixelType > void FiberExtractionFilter< PixelType >::GenerateData() { mitk::FiberBundle::Pointer fib = m_InputFiberBundle; if (fib->GetNumFibers()<=0) { MITK_INFO << "No fibers in tractogram!"; return; } if (m_Mode == MODE::OVERLAP) ExtractOverlap(fib); else if (m_Mode == MODE::ENDPOINTS) { if (m_InputType==INPUT::LABEL_MAP) ExtractLabels(fib); else ExtractEndpoints(fib); } } } #endif // __itkFiberExtractionFilter_cpp diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractDistanceFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractDistanceFilter.cpp index 6245c6bb4e..a86705a7d3 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractDistanceFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractDistanceFilter.cpp @@ -1,181 +1,189 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "itkTractDistanceFilter.h" #define _USE_MATH_DEFINES #include #include namespace itk{ TractDistanceFilter::TractDistanceFilter() : m_NumPoints(12) , disp(0) { } TractDistanceFilter::~TractDistanceFilter() { for (auto m : m_Metrics) delete m; } -std::vector TractDistanceFilter::GetIndices() const +vnl_matrix TractDistanceFilter::GetAllDistances() const { - return m_Indices; + return m_AllDistances; } -std::vector TractDistanceFilter::GetDistances() const +vnl_vector TractDistanceFilter::GetMinIndices() const { - return m_Distances; + return m_MinIndices; +} + +vnl_vector TractDistanceFilter::GetMinDistances() const +{ + return m_MinDistances; } void TractDistanceFilter::SetTracts2(const std::vector &Tracts2) { m_Tracts2 = Tracts2; } void TractDistanceFilter::SetTracts1(const std::vector &Tracts1) { m_Tracts1 = Tracts1; } void TractDistanceFilter::SetMetrics(const std::vector &Metrics) { m_Metrics = Metrics; } std::vector > TractDistanceFilter::ResampleFibers(mitk::FiberBundle::Pointer tractogram) { std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::FiberBundle::Pointer temp_fib = tractogram->GetDeepCopy(); temp_fib->ResampleToNumPoints(m_NumPoints); std::vector< vnl_matrix > out_fib; for (int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = temp_fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vnl_matrix streamline; streamline.set_size(3, m_NumPoints); streamline.fill(0.0); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< float, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; streamline.set_column(j, candV); } out_fib.push_back(streamline); } std::cout.rdbuf (old); // <-- restore return out_fib; } float TractDistanceFilter::calc_distance(const std::vector >& T1, const std::vector >& T2) { float distance = 0; for (auto metric : m_Metrics) { unsigned int r = 0; for (auto f1 : T1) { unsigned int c = 0; float min_d = 99999; for (auto f2 : T2) { #pragma omp critical ++disp; bool flipped = false; float d = metric->CalculateDistance(f1, f2, flipped); if (d < min_d) min_d = d; ++c; } distance += min_d; ++r; } } distance /= (T1.size() * m_Metrics.size()); return distance; } void TractDistanceFilter::GenerateData() { if (m_Metrics.empty()) { mitkThrow() << "No metric selected!"; return; } - m_Indices.clear(); - m_Distances.clear(); + m_MinIndices.clear(); + m_MinDistances.clear(); std::vector>> T1; std::vector>> T2; unsigned int num_fibs1 = 0; unsigned int num_fibs2 = 0; for (auto t : m_Tracts1) { T1.push_back(ResampleFibers(t)); num_fibs1 += T1.back().size(); } for (auto t : m_Tracts2) { T2.push_back(ResampleFibers(t)); num_fibs2 += T2.back().size(); } disp.restart(m_Metrics.size() * num_fibs1 * num_fibs2); - m_Indices.resize(T1.size()); - m_Distances.resize(T1.size()); + m_AllDistances.set_size(T1.size(), T2.size()); + + m_MinIndices.set_size(T1.size()); + m_MinDistances.set_size(T1.size()); + #pragma omp parallel for for (unsigned int i=0; i #include #include #include #include #include #include #include #include #include #include namespace itk{ /** * \brief */ class TractDistanceFilter : public ProcessObject { public: typedef TractDistanceFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef itk::Image< float, 3 > FloatImageType; typedef itk::Image< unsigned char, 3 > UcharImageType; itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkTypeMacro( TractDistanceFilter, ProcessObject ) itkSetMacro(NumPoints, unsigned int) ///< Fibers are resampled to the specified number of points. If scalar maps are used, a larger number of points is recommended. virtual void Update() override{ this->GenerateData(); } void SetMetrics(const std::vector &Metrics); void SetTracts1(const std::vector &Tracts1); void SetTracts2(const std::vector &Tracts2); - std::vector GetDistances() const; + vnl_vector GetMinDistances() const; - std::vector GetIndices() const; + vnl_vector GetMinIndices() const; + + vnl_matrix GetAllDistances() const; protected: void GenerateData() override; std::vector< vnl_matrix > ResampleFibers(FiberBundle::Pointer tractogram); float calc_distance(const std::vector > &T1, const std::vector > &T2); TractDistanceFilter(); virtual ~TractDistanceFilter(); - std::vector< FiberBundle::Pointer > m_Tracts1; - std::vector< FiberBundle::Pointer > m_Tracts2; - unsigned int m_NumPoints; - std::vector< mitk::ClusteringMetric* > m_Metrics; - std::vector< float > m_Distances; - std::vector< int > m_Indices; - boost::progress_display disp; + std::vector m_Tracts1; + std::vector m_Tracts2; + unsigned int m_NumPoints; + std::vector m_Metrics; + vnl_vector m_MinDistances; + vnl_vector m_MinIndices; + vnl_matrix m_AllDistances; + boost::progress_display disp; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkTractDistanceFilter.cpp" #endif #endif diff --git a/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt index 36f9917b3c..c0be517ec4 100644 --- a/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt +++ b/Modules/DiffusionImaging/FiberTracking/CMakeLists.txt @@ -1,52 +1,47 @@ set(_module_deps MitkDiffusionCore MitkGraphAlgorithms MitkCLVigraRandomForest) mitk_check_module_dependencies( MODULES ${_module_deps} MISSING_DEPENDENCIES_VAR _missing_deps ) if(NOT _missing_deps) set(lut_url http://mitk.org/download/data/FibertrackingLUT.tar.gz) set(lut_tarball ${CMAKE_CURRENT_BINARY_DIR}/FibertrackingLUT.tar.gz) file(DOWNLOAD ${lut_url} ${lut_tarball} EXPECTED_MD5 38ecb6d4a826c9ebb0f4965eb9aeee44 TIMEOUT 60 STATUS status SHOW_PROGRESS ) list(GET status 0 status_code) list(GET status 1 status_msg) if(NOT status_code EQUAL 0) message(SEND_ERROR "${status_msg} (error code ${status_code})") endif() file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/Resources) execute_process(COMMAND ${CMAKE_COMMAND} -E tar xzf ../FibertrackingLUT.tar.gz WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/Resources RESULT_VARIABLE result ERROR_VARIABLE err_msg) if(result) message(SEND_ERROR "Unpacking FibertrackingLUT.tar.gz failed: ${err_msg}") endif() endif() MITK_CREATE_MODULE( SUBPROJECTS MITK-DTI INCLUDE_DIRS Fiberfox Fiberfox/SignalModels Fiberfox/Sequences Algorithms Algorithms/TrackingHandlers Algorithms/ClusteringMetrics Algorithms/GibbsTracking Algorithms/StochasticTracking IODataStructures IODataStructures/FiberBundle IODataStructures/PlanarFigureComposite Rendering ${CMAKE_CURRENT_BINARY_DIR} DEPENDS ${_module_deps} PACKAGE_DEPENDS PUBLIC ITK|ITKFFT ITK|ITKDiffusionTensorImage Vigra HDF5 ) if(MODULE_IS_ENABLED) - add_subdirectory(cmdapps/Fiberfox) - add_subdirectory(cmdapps/FiberProcessing) - add_subdirectory(cmdapps/Tractography) - add_subdirectory(cmdapps/TractographyEvaluation) - add_subdirectory(cmdapps/Misc) add_subdirectory(Testing) endif() diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/PythonTest.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/PythonTest.cpp deleted file mode 100755 index b5c21e3952..0000000000 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Misc/PythonTest.cpp +++ /dev/null @@ -1,80 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center, -Division of Medical and Biological Informatics. -All rights reserved. - -This software is distributed WITHOUT ANY WARRANTY; without -even the implied warranty of MERCHANTABILITY or FITNESS FOR -A PARTICULAR PURPOSE. - -See LICENSE.txt or http://www.mitk.org for details. - -===================================================================*/ - -#include -#include - -#include -#include -#include - -#include -#include - -#include -#include -#include -#include - -typedef itk::Image< float, 3 > ItkImageType; - -/*! -\brief -*/ -int main(int argc, char* argv[]) -{ - mitkCommandLineParser parser; - parser.setArgumentPrefix("--", "-"); - parser.addArgument("image", "i", mitkCommandLineParser::InputFile, "Input image", "sh coefficient image", us::Any(), false); - parser.addArgument("script", "s", mitkCommandLineParser::InputFile, "", "", us::Any(), false); - parser.addArgument("params", "p", mitkCommandLineParser::InputFile, "", "", us::Any(), false); - parser.addArgument("out_file", "o", mitkCommandLineParser::OutputDirectory, "Output image", "output image", us::Any(), false); - - parser.setCategory("TEST"); - parser.setTitle("TEST"); - parser.setDescription("TEST"); - parser.setContributor("MIC"); - - std::map parsedArgs = parser.parseArguments(argc, argv); - if (parsedArgs.size()==0) - return EXIT_FAILURE; - - std::string image_file = us::any_cast(parsedArgs["image"]); - std::string script = us::any_cast(parsedArgs["script"]); - std::string params = us::any_cast(parsedArgs["params"]); - std::string out_file = us::any_cast(parsedArgs["out_file"]); - - mitk::Image::Pointer mitk_image = mitk::IOUtil::Load(image_file); - us::ModuleContext* context = us::GetModuleContext(); - - us::ServiceReference m_PythonServiceRef = context->GetServiceReference(); - mitk::IPythonService* m_PythonService = dynamic_cast ( context->GetService(m_PythonServiceRef) ); - - mitk::IPythonService::ForceLoadModule(); - - m_PythonService->Execute("import SimpleITK as sitk"); - m_PythonService->Execute("import SimpleITK._SimpleITK as _SimpleITK"); - m_PythonService->Execute("import numpy"); - - m_PythonService->CopyToPythonAsSimpleItkImage( mitk_image, "myvar"); - m_PythonService->Execute("myparams=\""+params+"\""); - m_PythonService->ExecuteScript(script); - - mitk::Image::Pointer out_seg = m_PythonService->CopySimpleItkImageFromPython("out_image"); - mitk::IOUtil::Save(out_seg, out_file); - - return EXIT_FAILURE; -} diff --git a/Modules/DiffusionImaging/MiniApps/ExtractAllGradients.cpp b/Modules/DiffusionImaging/MiniApps/ExtractAllGradients.cpp deleted file mode 100644 index 848d857171..0000000000 --- a/Modules/DiffusionImaging/MiniApps/ExtractAllGradients.cpp +++ /dev/null @@ -1,86 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center, -Division of Medical and Biological Informatics. -All rights reserved. - -This software is distributed WITHOUT ANY WARRANTY; without -even the implied warranty of MERCHANTABILITY or FITNESS FOR -A PARTICULAR PURPOSE. - -See LICENSE.txt or http://www.mitk.org for details. - -===================================================================*/ - -#include -#include -#include -#include "mitkCommandLineParser.h" -#include - -using namespace mitk; -using namespace std; - -/*! -\brief Copies transformation matrix of one image to another -*/ -int main(int argc, char* argv[]) -{ - typedef itk::VectorImage< short, 3 > ItkDwiType; - - mitkCommandLineParser parser; - - parser.setTitle("Extract all gradients"); - parser.setCategory("Preprocessing Tools"); - parser.setDescription("Extract all gradients from an diffusion image"); - parser.setContributor("MBI"); - - parser.setArgumentPrefix("--", "-"); - parser.addArgument("in", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(), false); - //parser.addArgument("extension", "e", mitkCommandLineParser::String, "File Extension:", "Extension of the output file", us::Any(), false); - //parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output:", "output image", us::Any(), false); - - map parsedArgs = parser.parseArguments(argc, argv); - if (parsedArgs.size() == 0) - return EXIT_FAILURE; - - MITK_INFO << "Extract parameter"; - // mandatory arguments - /*string inputName = us::any_cast(parsedArgs["in"]); - string extensionName = us::any_cast(parsedArgs["extension"]); - string ouputName = us::any_cast(parsedArgs["out"]);*/ - string in = us::any_cast(parsedArgs["in"]); - string inputName = "E:\\Kollektive\\R02-Lebertumore-Diffusion\\01-Extrahierte-Daten\\" + in + "\\" + in + "-DWI.dwi"; - string extensionName = ".nrrd"; - string ouputName = "E:\\Kollektive\\R02-Lebertumore-Diffusion\\01-Extrahierte-Daten\\" + in + "\\" + in + "-"; - - MITK_INFO << "Load Image: "; - mitk::Image::Pointer image = mitk::IOUtil::Load(inputName); - - //bool isDiffusionImage(mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(image)); - //if (!isDiffusionImage) - //{ - // MITK_INFO << "Input file is not of type diffusion image"; - // return; - //} - - ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); - mitk::CastToItkImage(image, itkVectorImagePointer); - - unsigned int channel = 0; - for (unsigned int channel = 0; channel < itkVectorImagePointer->GetVectorLength(); ++channel) - { - itk::ExtractDwiChannelFilter< short >::Pointer filter = itk::ExtractDwiChannelFilter< short >::New(); - filter->SetInput(itkVectorImagePointer); - filter->SetChannelIndex(channel); - filter->Update(); - - mitk::Image::Pointer newImage = mitk::Image::New(); - newImage->InitializeByItk(filter->GetOutput()); - newImage->SetImportChannel(filter->GetOutput()->GetBufferPointer()); - mitk::IOUtil::SaveImage(newImage, ouputName + to_string(channel) + extensionName); - } - return EXIT_SUCCESS; -} diff --git a/Modules/DiffusionImaging/Quantification/CMakeLists.txt b/Modules/DiffusionImaging/Quantification/CMakeLists.txt index 4f465c02c4..d3b6e845cc 100644 --- a/Modules/DiffusionImaging/Quantification/CMakeLists.txt +++ b/Modules/DiffusionImaging/Quantification/CMakeLists.txt @@ -1,15 +1,14 @@ # With apple gcc 4.2.1 the following waring leads to an build error if boost is enabled if(APPLE) mitkFunctionCheckCAndCXXCompilerFlags("-Wno-error=empty-body" CMAKE_C_FLAGS CMAKE_CXX_FLAGS) endif() #DiffusionImaging/Quantification MITK_CREATE_MODULE( SUBPROJECTS MITK-DTI INCLUDE_DIRS Algorithms IODataStructures IODataStructures/TbssImages Rendering ${CMAKE_CURRENT_BINARY_DIR} DEPENDS MitkDiffusionCore MitkFiberTracking MitkGraphAlgorithms PACKAGE_DEPENDS ITK|ITKThresholding ) add_subdirectory(Testing) -add_subdirectory(cmdapps) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp index 86635c518d..6098c451b7 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp @@ -1,1614 +1,1733 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkFiberProcessingView.h" -// Qt #include -// MITK #include #include #include #include #include #include #include #include #include #include #include "usModuleRegistry.h" #include #include "mitkNodePredicateDataType.h" #include #include #include #include -// ITK #include #include #include #include #include #include #include #include +#include +#include + const std::string QmitkFiberProcessingView::VIEW_ID = "org.mitk.views.fiberprocessing"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace mitk; QmitkFiberProcessingView::QmitkFiberProcessingView() : QmitkAbstractView() , m_Controls( 0 ) , m_CircleCounter(0) , m_PolygonCounter(0) , m_UpsamplingFactor(1) { } // Destructor QmitkFiberProcessingView::~QmitkFiberProcessingView() { RemoveObservers(); } void QmitkFiberProcessingView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberProcessingViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( OnDrawCircle() ) ); connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( OnDrawPolygon() ) ); connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(GenerateAndComposite()) ); connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) ); connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) ); connect(m_Controls->m_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) ); connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) ); connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) ); connect(m_Controls->m_CopyBundle, SIGNAL(clicked()), this, SLOT(CopyBundles()) ); connect(m_Controls->m_ExtractFibersButton, SIGNAL(clicked()), this, SLOT(Extract())); connect(m_Controls->m_RemoveButton, SIGNAL(clicked()), this, SLOT(Remove())); connect(m_Controls->m_ModifyButton, SIGNAL(clicked()), this, SLOT(Modify())); connect(m_Controls->m_ExtractionMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_RemovalMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ModificationMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ExtractionBoxMask, SIGNAL(currentIndexChanged(int)), this, SLOT(OnMaskExtractionChanged())); m_Controls->m_ColorMapBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isOdf = mitk::NodePredicateDataType::New("OdfImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isOdf); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); m_Controls->m_ColorMapBox->SetPredicate(finalPredicate); - - m_Controls->label_17->setVisible(false); - m_Controls->m_FiberExtractionFractionBox->setVisible(false); } UpdateGui(); + OnMaskExtractionChanged(); } void QmitkFiberProcessingView::OnMaskExtractionChanged() { + m_Controls->m_FiberExtractionFractionLabel->setVisible(false); + m_Controls->m_FiberExtractionFractionBox->setVisible(false); + + m_Controls->m_FiberExtractionThresholdLabel->setVisible(false); + m_Controls->m_FiberExtractionThresholdBox->setVisible(false); + + m_Controls->m_InterpolateRoiBox->setVisible(false); + m_Controls->m_BothEnds->setVisible(false); + + m_Controls->m_LabelsBox->setVisible(false); + m_Controls->m_LabelsLabel->setVisible(false); + if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2 || m_Controls->m_ExtractionBoxMask->currentIndex() == 3) { - m_Controls->label_17->setVisible(true); + m_Controls->m_FiberExtractionFractionLabel->setVisible(true); m_Controls->m_FiberExtractionFractionBox->setVisible(true); - m_Controls->m_BothEnds->setVisible(false); + + m_Controls->m_FiberExtractionThresholdLabel->setVisible(true); + m_Controls->m_FiberExtractionThresholdBox->setVisible(true); + + m_Controls->m_InterpolateRoiBox->setVisible(true); } - else + else if (m_Controls->m_ExtractionBoxMask->currentIndex() == 0 || m_Controls->m_ExtractionBoxMask->currentIndex() == 1) { - m_Controls->label_17->setVisible(false); - m_Controls->m_FiberExtractionFractionBox->setVisible(false); if (m_Controls->m_ExtractionBoxMask->currentIndex() != 3) m_Controls->m_BothEnds->setVisible(true); + m_Controls->m_InterpolateRoiBox->setVisible(true); + + m_Controls->m_FiberExtractionThresholdLabel->setVisible(true); + m_Controls->m_FiberExtractionThresholdBox->setVisible(true); + } + else if (m_Controls->m_ExtractionBoxMask->currentIndex() == 4) + { + m_Controls->m_BothEnds->setVisible(true); + + m_Controls->m_LabelsBox->setVisible(true); + m_Controls->m_LabelsLabel->setVisible(true); } } void QmitkFiberProcessingView::SetFocus() { m_Controls->toolBoxx->setFocus(); } void QmitkFiberProcessingView::Modify() { switch (m_Controls->m_ModificationMethodBox->currentIndex()) { case 0: { ResampleSelectedBundlesSpline(); break; } case 1: { ResampleSelectedBundlesLinear(); break; } case 2: { CompressSelectedBundles(); break; } case 3: { DoImageColorCoding(); break; } case 4: { MirrorFibers(); break; } case 5: { WeightFibers(); break; } case 6: { DoCurvatureColorCoding(); break; } case 7: { DoWeightColorCoding(); break; } case 8: { DoLengthColorCoding(); break; } } } void QmitkFiberProcessingView::WeightFibers() { float weight = this->m_Controls->m_BundleWeightBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->SetFiberWeights(weight); } } void QmitkFiberProcessingView::Remove() { switch (m_Controls->m_RemovalMethodBox->currentIndex()) { case 0: { RemoveDir(); break; } case 1: { PruneBundle(); break; } case 2: { ApplyCurvatureThreshold(); break; } case 3: { RemoveWithMask(false); break; } case 4: { RemoveWithMask(true); break; } case 5: { ApplyWeightThreshold(); break; } + case 6: + { + ApplyDensityThreshold(); + break; + } } } void QmitkFiberProcessingView::Extract() { switch (m_Controls->m_ExtractionMethodBox->currentIndex()) { case 0: { ExtractWithPlanarFigure(); break; } case 1: { switch (m_Controls->m_ExtractionBoxMask->currentIndex()) { { case 0: - ExtractWithMask(true, false); + ExtractWithMask(true, false, false); break; } { case 1: - ExtractWithMask(true, true); + ExtractWithMask(true, true, false); break; } { case 2: - ExtractWithMask(false, false); + ExtractWithMask(false, false, false); break; } { case 3: - ExtractWithMask(false, true); + ExtractWithMask(false, true, false); + break; + } + { + case 4: + ExtractWithMask(true, false, true); break; } } break; } } } void QmitkFiberProcessingView::PruneBundle() { int minLength = this->m_Controls->m_PruneFibersMinBox->value(); int maxLength = this->m_Controls->m_PruneFibersMaxBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); if (!fib->RemoveShortFibers(minLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); else if (!fib->RemoveLongFibers(maxLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyWeightThreshold() { float thr = this->m_Controls->m_WeightThresholdBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->FilterByWeights(thr); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByFiberWeights(false, true); node->SetData(newFib); } else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } +void QmitkFiberProcessingView::ApplyDensityThreshold() +{ + float thr = this->m_Controls->m_DensityThresholdBox->value(); + float ol = this->m_Controls->m_DensityOverlapBox->value(); + std::vector< DataNode::Pointer > nodes = m_SelectedFB; + for (auto node : nodes) + { + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + + itk::TractDensityImageFilter< ItkFloatImageType >::Pointer generator = itk::TractDensityImageFilter< ItkFloatImageType >::New(); + generator->SetFiberBundle(fib); + generator->SetBinaryOutput(false); + generator->SetOutputAbsoluteValues(false); + generator->Update(); + + itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); + extractor->SetRoiImages({generator->GetOutput()}); + extractor->SetInputFiberBundle(fib); + extractor->SetOverlapFraction(ol); + extractor->SetInterpolate(true); + extractor->SetThreshold(thr); + extractor->SetNoNegatives(true); + extractor->Update(); + + if (extractor->GetPositives().empty()) + { + QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); + continue; + } + + mitk::FiberBundle::Pointer newFib = extractor->GetPositives().at(0); + if (newFib->GetNumFibers()>0) + node->SetData(newFib); + else + QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); + } + RenderingManager::GetInstance()->RequestUpdateAll(); +} + void QmitkFiberProcessingView::ApplyCurvatureThreshold() { int angle = this->m_Controls->m_CurvSpinBox->value(); int dist = this->m_Controls->m_CurvDistanceSpinBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); itk::FiberCurvatureFilter::Pointer filter = itk::FiberCurvatureFilter::New(); filter->SetInputFiberBundle(fib); filter->SetAngularDeviation(angle); filter->SetDistance(dist); filter->SetRemoveFibers(m_Controls->m_RemoveCurvedFibersBox->isChecked()); filter->Update(); mitk::FiberBundle::Pointer newFib = filter->GetOutputFiberBundle(); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByOrientation(); node->SetData(newFib); } else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveDir() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); vnl_vector_fixed dir; dir[0] = m_Controls->m_ExtractDirX->value(); dir[1] = m_Controls->m_ExtractDirY->value(); dir[2] = m_Controls->m_ExtractDirZ->value(); fib->RemoveDir(dir,cos((float)m_Controls->m_ExtractAngle->value()*itk::Math::pi/180)); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveWithMask(bool removeInside) { if (m_RoiImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_RoiImageNode->GetData()); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(mitkMask, mask); mitk::FiberBundle::Pointer newFib = fib->RemoveFibersOutside(mask, removeInside); if (newFib->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } node->SetData(newFib); } RenderingManager::GetInstance()->RequestUpdateAll(); } -void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert) +void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert, bool labelmap) { if (m_RoiImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_RoiImageNode->GetData()); for (auto node : m_SelectedFB) { + std::string roi_name = m_RoiImageNode->GetName(); mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); - QString name(node->GetName().c_str()); - ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(mitkMask, mask); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(fib); extractor->SetRoiImages({mask}); + extractor->SetRoiImageNames({roi_name}); extractor->SetThreshold(m_Controls->m_FiberExtractionThresholdBox->value()); extractor->SetOverlapFraction(m_Controls->m_FiberExtractionFractionBox->value()); extractor->SetBothEnds(m_Controls->m_BothEnds->isChecked()); extractor->SetInterpolate(m_Controls->m_InterpolateRoiBox->isChecked()); + extractor->SetMinFibersPerTract(m_Controls->m_MinExtractedFibersBox->value()); + extractor->SetSplitByRoi(true); if (invert) extractor->SetNoPositives(true); else extractor->SetNoNegatives(true); + + if (labelmap) + { + std::string labels_string = m_Controls->m_LabelsBox->text().toStdString(); + if (labels_string!="ALL") + { + std::vector strs; + boost::split(strs,labels_string,boost::is_any_of(" ,;\t")); + + std::vector< unsigned short > labels_vector; + for (auto v : strs) + { + try{ + unsigned short l = boost::lexical_cast(v); + labels_vector.push_back(l); + } + catch(...) + { + + } + } + extractor->SetLabels(labels_vector); + } + + extractor->SetInterpolate(false); + extractor->SetInputType(itk::FiberExtractionFilter::INPUT::LABEL_MAP); + extractor->SetSplitLabels(true); + onlyEnds = true; + } + if (onlyEnds) extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); extractor->Update(); - mitk::FiberBundle::Pointer newFib; + std::vector< mitk::FiberBundle::Pointer > newFibs; if (invert) - newFib = extractor->GetNegatives().at(0); + newFibs = extractor->GetNegatives(); else - newFib = extractor->GetPositives().at(0); + newFibs = extractor->GetPositives(); - if (newFib.IsNull() || newFib->GetNumFibers()<=0) + if (newFibs.empty()) { - QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); + QMessageBox::information(nullptr, "No output generated:", "No fibers could be extracted."); continue; } - DataNode::Pointer newNode = DataNode::New(); - newNode->SetData(newFib); + auto labels = extractor->GetPositiveLabels(); + MITK_INFO << labels.size(); - if (invert) + for (unsigned int i=0; iSetData(fib); + std::string name = roi_name; + if (iSetName(name.toStdString()); - GetDataStorage()->Add(newNode); + newNode->SetName(name); + GetDataStorage()->Add(newNode, node); + } node->SetVisibility(false); } } void QmitkFiberProcessingView::GenerateRoiImage() { if (m_SelectedPF.empty()) return; mitk::BaseGeometry::Pointer geometry; if (!m_SelectedFB.empty()) { mitk::FiberBundle::Pointer fib = dynamic_cast(m_SelectedFB.front()->GetData()); geometry = fib->GetGeometry(); } else if (m_SelectedImage) geometry = m_SelectedImage->GetGeometry(); else return; itk::Vector spacing = geometry->GetSpacing(); spacing /= m_UpsamplingFactor; mitk::Point3D newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); newOrigin[0] += bounds.GetElement(0); newOrigin[1] += bounds.GetElement(2); newOrigin[2] += bounds.GetElement(4); itk::Matrix direction; itk::ImageRegion<3> imageRegion; for (int i=0; i<3; i++) for (int j=0; j<3; j++) direction[j][i] = geometry->GetMatrixColumn(i)[j]/spacing[j]; imageRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); imageRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); imageRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); m_PlanarFigureImage = ItkUCharImageType::New(); m_PlanarFigureImage->SetSpacing( spacing ); // Set the image spacing m_PlanarFigureImage->SetOrigin( newOrigin ); // Set the image origin m_PlanarFigureImage->SetDirection( direction ); // Set the image direction m_PlanarFigureImage->SetRegions( imageRegion ); m_PlanarFigureImage->Allocate(); m_PlanarFigureImage->FillBuffer( 0 ); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); std::string name = m_SelectedPF.at(0)->GetName(); WritePfToImage(m_SelectedPF.at(0), tmpImage); for (unsigned int i=1; iGetName(); WritePfToImage(m_SelectedPF.at(i), tmpImage); } DataNode::Pointer node = DataNode::New(); tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); node->SetData(tmpImage); node->SetName(name); this->GetDataStorage()->Add(node); } void QmitkFiberProcessingView::WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image) { if (dynamic_cast(node->GetData())) { m_PlanarFigure = dynamic_cast(node->GetData()); AccessFixedDimensionByItk_2( image, InternalReorientImagePlane, 3, m_PlanarFigure->GetGeometry(), -1); AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 3, 2, node->GetName() ); } else if (dynamic_cast(node->GetData())) { DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node); for (unsigned int i=0; iSize(); i++) { WritePfToImage(children->at(i), image); } } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< float, VImageDimension > FloatImageType; typedef itk::ResampleImageFilter ResamplerType; typename ResamplerType::Pointer resampler = ResamplerType::New(); mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); float upsamp = m_UpsamplingFactor; float gausssigma = 0.5; // Spacing typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); spacing[0] = image->GetSpacing()[0] / upsamp; spacing[1] = image->GetSpacing()[1] / upsamp; spacing[2] = image->GetSpacing()[2]; resampler->SetOutputSpacing( spacing ); // Size typename ResamplerType::SizeType size; size[0] = planegeo->GetExtentInMM(0) / spacing[0]; size[1] = planegeo->GetExtentInMM(1) / spacing[1]; size[2] = 1; resampler->SetSize( size ); // Origin typename mitk::Point3D orig = planegeo->GetOrigin(); typename mitk::Point3D corrorig; planegeo3D->WorldToIndex(orig,corrorig); corrorig[0] += 0.5/upsamp; corrorig[1] += 0.5/upsamp; corrorig[2] += 0; planegeo3D->IndexToWorld(corrorig,corrorig); resampler->SetOutputOrigin(corrorig ); // Direction typename ResamplerType::DirectionType direction; typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); for(unsigned int c=0; cSetOutputDirection( direction ); // Gaussian interpolation if(gausssigma != 0) { double sigma[3]; for( unsigned int d = 0; d < 3; d++ ) sigma[d] = gausssigma * image->GetSpacing()[d]; double alpha = 2.0; typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New(); interpolator->SetInputImage( image ); interpolator->SetParameters( sigma, alpha ); resampler->SetInterpolator( interpolator ); } else { typedef typename itk::LinearInterpolateImageFunction InterpolatorType; typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); } resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); resampler->Update(); if(additionalIndex < 0) { this->m_InternalImage = mitk::Image::New(); this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string ) { typedef itk::Image< TPixel, VImageDimension > ImageType; // Generate mask image as new image with same header as input image and // initialize with "1". ItkUCharImageType::Pointer newMaskImage = ItkUCharImageType::New(); newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); newMaskImage->Allocate(); newMaskImage->FillBuffer( 1 ); // Generate VTK polygon from (closed) PlanarFigure polyline // (The polyline points are shifted by -0.5 in z-direction to make sure // that the extrusion filter, which afterwards elevates all points by +0.5 // in z-direction, creates a 3D object which is cut by the the plane z=0) const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const BaseGeometry *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); vtkPolyData *polyline = vtkPolyData::New(); polyline->Allocate( 1, 1 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // Create VTK polydata object of polyline contour vtkPoints *points = vtkPoints::New(); PlanarFigure::PolyLineType::const_iterator it; unsigned int numberOfPoints = 0; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected image Point2D point2D = *it; planarFigurePlaneGeometry->WorldToIndex(point2D, point2D); point2D[0] -= 0.5/m_UpsamplingFactor; point2D[1] -= 0.5/m_UpsamplingFactor; planarFigurePlaneGeometry->IndexToWorld(point2D, point2D); planarFigurePlaneGeometry->Map( point2D, point3D ); // Polygons (partially) outside of the image bounds can not be processed further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { float bounds[2] = {0,0}; bounds[0] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i0); bounds[1] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i1); imageGeometry3D->WorldToIndex( point3D, point3D ); if (point3D[i0]<0) point3D[i0] = 0.0; else if (point3D[i0]>bounds[0]) point3D[i0] = bounds[0]-0.001; if (point3D[i1]<0) point3D[i1] = 0.0; else if (point3D[i1]>bounds[1]) point3D[i1] = bounds[1]-0.001; points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } else { imageGeometry3D->WorldToIndex( point3D, point3D ); // Add point to polyline array points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } } polyline->SetPoints( points ); points->Delete(); vtkIdType *ptIds = new vtkIdType[numberOfPoints]; for ( vtkIdType i = 0; i < numberOfPoints; ++i ) ptIds[i] = i; polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); // Extrude the generated contour polygon vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); extrudeFilter->SetInputData( polyline ); extrudeFilter->SetScaleFactor( 1 ); extrudeFilter->SetExtrusionTypeToNormalExtrusion(); extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); // Make a stencil from the extruded polygon vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); polyDataToImageStencil->SetInputConnection( extrudeFilter->GetOutputPort() ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< ItkUCharImageType > ImageImportType; typedef itk::VTKImageExport< ItkUCharImageType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( newMaskImage ); vtkImageImport *vtkImporter = vtkImageImport::New(); this->ConnectPipelines( itkExporter, vtkImporter ); vtkImporter->Update(); // Apply the generated image stencil to the input image vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(polyDataToImageStencil->GetOutputPort() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkImageExport *vtkExporter = vtkImageExport::New(); vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // calculate cropping bounding box m_InternalImageMask3D = itkImporter->GetOutput(); m_InternalImageMask3D->SetDirection(image->GetDirection()); itk::ImageRegionConstIterator itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); itk::ImageRegionIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); typename ImageType::SizeType lowersize = {{itk::NumericTraits::max(),itk::NumericTraits::max(),itk::NumericTraits::max()}}; typename ImageType::SizeType uppersize = {{0,0,0}}; while( !itmask.IsAtEnd() ) { if(itmask.Get() == 0) itimage.Set(0); else { typename ImageType::IndexType index = itimage.GetIndex(); typename ImageType::SizeType signedindex; signedindex[0] = index[0]; signedindex[1] = index[1]; signedindex[2] = index[2]; lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; } ++itmask; ++itimage; } typename ImageType::IndexType index; index[0] = lowersize[0]; index[1] = lowersize[1]; index[2] = lowersize[2]; typename ImageType::SizeType size; size[0] = uppersize[0] - lowersize[0] + 1; size[1] = uppersize[1] - lowersize[1] + 1; size[2] = uppersize[2] - lowersize[2] + 1; itk::ImageRegion<3> cropRegion = itk::ImageRegion<3>(index, size); // crop internal mask typedef itk::RegionOfInterestImageFilter< ItkUCharImageType, ItkUCharImageType > ROIMaskFilterType; typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); roi2->SetRegionOfInterest(cropRegion); roi2->SetInput(m_InternalImageMask3D); roi2->Update(); m_InternalImageMask3D = roi2->GetOutput(); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer()); tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer()); Image::Pointer tmpImage2 = Image::New(); tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer()); const BaseGeometry *pfImageGeometry3D = tmpImage2->GetGeometry( 0 ); const BaseGeometry *intImageGeometry3D = tmpImage->GetGeometry( 0 ); typedef itk::ImageRegionIteratorWithIndex IteratorType; IteratorType imageIterator (m_InternalImageMask3D, m_InternalImageMask3D->GetRequestedRegion()); imageIterator.GoToBegin(); while ( !imageIterator.IsAtEnd() ) { unsigned char val = imageIterator.Value(); if (val>0) { itk::Index<3> index = imageIterator.GetIndex(); Point3D point; point[0] = index[0]; point[1] = index[1]; point[2] = index[2]; intImageGeometry3D->IndexToWorld(point, point); pfImageGeometry3D->WorldToIndex(point, point); point[i0] += 0.5; point[i1] += 0.5; index[0] = point[0]; index[1] = point[1]; index[2] = point[2]; if (pfImageGeometry3D->IsIndexInside(index)) m_PlanarFigureImage->SetPixel(index, 1); } ++imageIterator; } // Clean up VTK objects polyline->Delete(); extrudeFilter->Delete(); polyDataToImageStencil->Delete(); vtkImporter->Delete(); imageStencilFilter->Delete(); //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? delete[] ptIds; } void QmitkFiberProcessingView::UpdateGui() { m_Controls->m_FibLabel->setText("mandatory"); m_Controls->m_PfLabel->setText("needed for extraction"); m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->m_RemoveButton->setEnabled(false); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); m_Controls->PFCompoANDButton->setEnabled(false); m_Controls->PFCompoORButton->setEnabled(false); m_Controls->PFCompoNOTButton->setEnabled(false); m_Controls->m_GenerateRoiImage->setEnabled(false); m_Controls->m_ExtractFibersButton->setEnabled(false); m_Controls->m_ModifyButton->setEnabled(false); m_Controls->m_CopyBundle->setEnabled(false); m_Controls->m_JoinBundles->setEnabled(false); m_Controls->m_SubstractBundles->setEnabled(false); // disable alle frames m_Controls->m_BundleWeightFrame->setVisible(false); m_Controls->m_ExtactionFramePF->setVisible(false); m_Controls->m_RemoveDirectionFrame->setVisible(false); m_Controls->m_RemoveLengthFrame->setVisible(false); m_Controls->m_RemoveCurvatureFrame->setVisible(false); m_Controls->m_RemoveByWeightFrame->setVisible(false); + m_Controls->m_RemoveByDensityFrame->setVisible(false); m_Controls->m_SmoothFibersFrame->setVisible(false); m_Controls->m_CompressFibersFrame->setVisible(false); m_Controls->m_ColorFibersFrame->setVisible(false); m_Controls->m_MirrorFibersFrame->setVisible(false); m_Controls->m_MaskExtractionFrame->setVisible(false); m_Controls->m_ColorMapBox->setVisible(false); bool pfSelected = !m_SelectedPF.empty(); bool fibSelected = !m_SelectedFB.empty(); bool multipleFibsSelected = (m_SelectedFB.size()>1); bool maskSelected = m_RoiImageNode.IsNotNull(); bool imageSelected = m_SelectedImage.IsNotNull(); // toggle visibility of elements according to selected method switch ( m_Controls->m_ExtractionMethodBox->currentIndex() ) { case 0: m_Controls->m_ExtactionFramePF->setVisible(true); break; case 1: m_Controls->m_MaskExtractionFrame->setVisible(true); break; } switch ( m_Controls->m_RemovalMethodBox->currentIndex() ) { case 0: m_Controls->m_RemoveDirectionFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 1: m_Controls->m_RemoveLengthFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 2: m_Controls->m_RemoveCurvatureFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 3: break; case 4: break; case 5: m_Controls->m_RemoveByWeightFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; + case 6: + m_Controls->m_RemoveByDensityFrame->setVisible(true); + if ( fibSelected ) + m_Controls->m_RemoveButton->setEnabled(true); + break; } switch ( m_Controls->m_ModificationMethodBox->currentIndex() ) { case 0: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 1: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 2: m_Controls->m_CompressFibersFrame->setVisible(true); break; case 3: m_Controls->m_ColorFibersFrame->setVisible(true); m_Controls->m_ColorMapBox->setVisible(true); break; case 4: m_Controls->m_MirrorFibersFrame->setVisible(true); if (m_SelectedSurfaces.size()>0) m_Controls->m_ModifyButton->setEnabled(true); break; case 5: m_Controls->m_BundleWeightFrame->setVisible(true); break; case 6: m_Controls->m_ColorFibersFrame->setVisible(true); break; case 7: m_Controls->m_ColorFibersFrame->setVisible(true); break; case 8: m_Controls->m_ColorFibersFrame->setVisible(true); break; } // are fiber bundles selected? if ( fibSelected ) { m_Controls->m_CopyBundle->setEnabled(true); m_Controls->m_ModifyButton->setEnabled(true); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); m_Controls->m_FibLabel->setText(QString(m_SelectedFB.at(0)->GetName().c_str())); // one bundle and one planar figure needed to extract fibers if (pfSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==0) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_SelectedPF.at(0)->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } // more than two bundles needed to join/subtract if (multipleFibsSelected) { m_Controls->m_FibLabel->setText("multiple bundles selected"); m_Controls->m_JoinBundles->setEnabled(true); m_Controls->m_SubstractBundles->setEnabled(true); } if (maskSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==1) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_RoiImageNode->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } if (maskSelected && (m_Controls->m_RemovalMethodBox->currentIndex()==3 || m_Controls->m_RemovalMethodBox->currentIndex()==4) ) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_RoiImageNode->GetName().c_str())); m_Controls->m_RemoveButton->setEnabled(true); } } // are planar figures selected? if (pfSelected) { if ( fibSelected || m_SelectedImage.IsNotNull()) m_Controls->m_GenerateRoiImage->setEnabled(true); if (m_SelectedPF.size() > 1) { m_Controls->PFCompoANDButton->setEnabled(true); m_Controls->PFCompoORButton->setEnabled(true); } else m_Controls->PFCompoNOTButton->setEnabled(true); } // is image selected if (imageSelected || maskSelected) { m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); } } void QmitkFiberProcessingView::NodeRemoved(const mitk::DataNode* node ) { for (auto fnode: m_SelectedFB) if (node == fnode) { m_SelectedFB.clear(); break; } berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } void QmitkFiberProcessingView::NodeAdded(const mitk::DataNode* ) { if (!m_Controls->m_InteractiveBox->isChecked()) { berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } } void QmitkFiberProcessingView::OnEndInteraction() { if (m_Controls->m_InteractiveBox->isChecked()) ExtractWithPlanarFigure(true); } void QmitkFiberProcessingView::AddObservers() { typedef itk::SimpleMemberCommand< QmitkFiberProcessingView > SimpleCommandType; for (auto node : m_SelectedPF) { mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); if (figure!=nullptr) { figure->RemoveAllObservers(); // add observer for event when interaction with figure starts SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); endInteractionCommand->SetCallbackFunction( this, &QmitkFiberProcessingView::OnEndInteraction); m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); } } } void QmitkFiberProcessingView::RemoveObservers() { for (auto node : m_SelectedPF) { mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); if (figure!=nullptr) figure->RemoveAllObservers(); } } void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { RemoveObservers(); //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection std::vector lastSelectedFB = m_SelectedFB; m_SelectedFB.clear(); m_SelectedPF.clear(); m_SelectedSurfaces.clear(); m_SelectedImage = nullptr; m_RoiImageNode = nullptr; for (auto node: nodes) { if ( dynamic_cast(node->GetData()) ) m_SelectedFB.push_back(node); else if (dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData())) m_SelectedPF.push_back(node); else if (dynamic_cast(node->GetData())) { m_SelectedImage = dynamic_cast(node->GetData()); if (m_SelectedImage->GetDimension()==3) m_RoiImageNode = node; } else if (dynamic_cast(node->GetData())) m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); } // if we perform interactive fiber extraction, we want to avoid auto-selection of the extracted bundle if (m_SelectedFB.empty() && m_Controls->m_InteractiveBox->isChecked()) m_SelectedFB = lastSelectedFB; // if no fibers or surfaces are selected, select topmost if (m_SelectedFB.empty() && m_SelectedSurfaces.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedFB.clear(); m_SelectedFB.push_back(nodes->at(i)); } } } // if no plar figure is selected, select topmost if (m_SelectedPF.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedPF.clear(); m_SelectedPF.push_back(nodes->at(i)); } } } AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::OnDrawPolygon() { mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); figure->ClosedOn(); this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); } void QmitkFiberProcessingView::OnDrawCircle() { mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter)); } void QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *, mitk::BaseProperty* ) { // initialize figure's geometry with empty geometry mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New(); figure->SetPlaneGeometry( emptygeometry ); //set desired data to DataNode where Planarfigure is stored mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetName(name.toStdString()); newNode->SetData(figure); newNode->SetBoolProperty("planarfigure.3drendering", true); newNode->SetBoolProperty("planarfigure.3drendering.fill", true); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(newNode->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode(newNode); } // figure drawn on the topmost layer / image GetDataStorage()->Add(newNode ); RemoveObservers(); for(unsigned int i = 0; i < m_SelectedPF.size(); i++) m_SelectedPF[i]->SetSelected(false); newNode->SetSelected(true); m_SelectedPF.clear(); m_SelectedPF.push_back(newNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::ExtractWithPlanarFigure(bool interactive) { if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){ QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); return; } try { std::vector fiberBundles = m_SelectedFB; mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0); for (unsigned int i=0; i(fiberBundles.at(i)->GetData()); mitk::FiberBundle::Pointer extFB = fib->ExtractFiberSubset(planarFigure, GetDataStorage()); if (interactive && m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } float op = 5.0/sqrt(fib->GetNumFibers()); float currentOp = 0; fiberBundles.at(i)->GetFloatProperty("opacity", currentOp); if (currentOp!=op) { fib->SetFiberColors(255, 255, 255); fiberBundles.at(i)->SetFloatProperty("opacity", op); fiberBundles.at(i)->SetBoolProperty("Fiber2DfadeEFX", false); } m_InteractiveNode->SetData(extFB); } else { if (extFB->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } mitk::DataNode::Pointer node; node = mitk::DataNode::New(); node->SetData(extFB); QString name(fiberBundles.at(i)->GetName().c_str()); name += "*"; node->SetName(name.toStdString()); fiberBundles.at(i)->SetVisibility(false); GetDataStorage()->Add(node); } } } catch(const std::out_of_range& ) { QMessageBox::warning( nullptr, "Fiber extraction failed", "Did you only create the planar figure, using the circle or polygon button, but forgot to actually place it in the image afterwards? \nAfter creating a planar figure, simply left-click at the desired position in the image or on the tractogram to place it."); } } void QmitkFiberProcessingView::GenerateAndComposite() { mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); PFCAnd->setOperationType(mitk::PlanarFigureComposite::AND); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("AND"); newPFCNode->SetData(PFCAnd); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); RemoveObservers(); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::GenerateOrComposite() { mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); PFCOr->setOperationType(mitk::PlanarFigureComposite::OR); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("OR"); newPFCNode->SetData(PFCOr); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); UpdateGui(); } void QmitkFiberProcessingView::GenerateNotComposite() { mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); PFCNot->setOperationType(mitk::PlanarFigureComposite::NOT); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("NOT"); newPFCNode->SetData(PFCNot); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode ) { pfc->SetSelected(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(pfc, parentNode); else GetDataStorage()->Add(pfc); for (auto child : children) { if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); newChild->SetBoolProperty("planarfigure.3drendering", true); newChild->SetBoolProperty("planarfigure.3drendering.fill", true); GetDataStorage()->Add(newChild, pfc); GetDataStorage()->Remove(child); } else if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); std::vector< mitk::DataNode::Pointer > grandChildVector; mitk::DataStorage::SetOfObjects::ConstPointer grandchildren = GetDataStorage()->GetDerivations(child); for( mitk::DataStorage::SetOfObjects::const_iterator it = grandchildren->begin(); it != grandchildren->end(); ++it ) grandChildVector.push_back(*it); AddCompositeToDatastorage(newChild, grandChildVector, pfc); GetDataStorage()->Remove(child); } } UpdateGui(); } void QmitkFiberProcessingView::CopyBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( nullptr, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->GetDeepCopy(); node->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newFib); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); } UpdateGui(); } void QmitkFiberProcessingView::JoinBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } m_SelectedFB.at(0)->SetVisibility(false); mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); std::vector< mitk::FiberBundle::Pointer > tractograms; for (unsigned int i=1; iSetVisibility(false); tractograms.push_back(dynamic_cast(m_SelectedFB.at(i)->GetData())); } newBundle = newBundle->AddBundles(tractograms); mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName("Joined_Tractograms"); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } void QmitkFiberProcessingView::SubstractBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); m_SelectedFB.at(0)->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); for (unsigned int i=1; iSubtractBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); if (newBundle.IsNull()) break; name += "-"+QString(m_SelectedFB.at(i)->GetName().c_str()); m_SelectedFB.at(i)->SetVisibility(false); } if (newBundle.IsNull()) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers. Did you select the fiber bundles in the correct order? X-Y is not equal to Y-X!"); return; } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } void QmitkFiberProcessingView::ResampleSelectedBundlesSpline() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleSpline(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ResampleSelectedBundlesLinear() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleLinear(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::CompressSelectedBundles() { double factor = this->m_Controls->m_ErrorThresholdBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->Compress(factor); fib->ColorFibersByOrientation(); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::DoImageColorCoding() { if (m_Controls->m_ColorMapBox->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Bundle coloring aborted:", "No image providing the scalar values for coloring the selected bundle available."); return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByScalarMap(dynamic_cast(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoCurvatureColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByCurvature(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoLengthColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByLength(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoWeightColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByFiberWeights(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::MirrorFibers() { unsigned int axis = this->m_Controls->m_MirrorFibersBox->currentIndex(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); if (m_SelectedImage.IsNotNull()) fib->SetReferenceGeometry(m_SelectedImage->GetGeometry()); fib->MirrorFibers(axis); } for (auto surf : m_SelectedSurfaces) { vtkSmartPointer poly = surf->GetVtkPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); for (int i=0; iGetNumberOfPoints(); i++) { double* point = poly->GetPoint(i); point[axis] *= -1; vtkNewPoints->InsertNextPoint(point); } poly->SetPoints(vtkNewPoints); surf->CalculateBoundingBox(); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h index 07e9b8660b..50e7649deb 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.h @@ -1,202 +1,203 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef QmitkFiberProcessingView_h #define QmitkFiberProcessingView_h #include #include "ui_QmitkFiberProcessingViewControls.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /*! \brief View to process fiber bundles. Supplies methods to extract fibers from the bundle, fiber resampling, mirroring, join and subtract bundles and much more. */ class QmitkFiberProcessingView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: typedef itk::Image< unsigned char, 3 > ItkUCharImageType; typedef itk::Image< float, 3 > ItkFloatImageType; static const std::string VIEW_ID; QmitkFiberProcessingView(); virtual ~QmitkFiberProcessingView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; protected slots: void OnDrawCircle(); ///< add circle interactors etc. void OnDrawPolygon(); ///< add circle interactors etc. void GenerateAndComposite(); void GenerateOrComposite(); void GenerateNotComposite(); void CopyBundles(); ///< add copies of selected bundles to data storage void JoinBundles(); ///< merge selected fiber bundles void SubstractBundles(); ///< subtract bundle A from bundle B. Not commutative! Defined by order of selection. void GenerateRoiImage(); ///< generate binary image of selected planar figures. void Remove(); void Extract(); void Modify(); void UpdateGui(); ///< update button activity etc. dpending on current datamanager selection void OnMaskExtractionChanged(); virtual void AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *propertyKey = nullptr, mitk::BaseProperty *property = nullptr ); protected: void MirrorFibers(); ///< mirror bundle on the specified plane void ResampleSelectedBundlesSpline(); ///< void ResampleSelectedBundlesLinear(); ///< void DoImageColorCoding(); ///< color fibers by selected scalar image void DoWeightColorCoding(); ///< color fibers by their respective weights void DoLengthColorCoding(); ///< color fibers by length void DoCurvatureColorCoding(); ///< color fibers by curvature void CompressSelectedBundles(); ///< remove points below certain error threshold void WeightFibers(); void ApplyWeightThreshold(); + void ApplyDensityThreshold(); void RemoveWithMask(bool removeInside); void RemoveDir(); void ApplyCurvatureThreshold(); ///< remove/split fibers with a too high curvature threshold void PruneBundle(); ///< remove too short/too long fibers - void ExtractWithMask(bool onlyEnds, bool invert); + void ExtractWithMask(bool onlyEnds, bool invert, bool labelmap); void ExtractWithPlanarFigure(bool interactive=false); void OnEndInteraction(); /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkFiberProcessingViewControls* m_Controls; /** Connection from VTK to ITK */ template void ConnectPipelines(VTK_Exporter* exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } template void ConnectPipelines(ITK_Exporter exporter, VTK_Importer* importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string nodeName ); template < typename TPixel, unsigned int VImageDimension > void InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex ); int m_CircleCounter; ///< used for data node naming int m_PolygonCounter; ///< used for data node naming std::vector m_SelectedFB; ///< selected fiber bundle nodes std::vector m_SelectedPF; ///< selected planar figure nodes std::vector m_SelectedSurfaces; mitk::Image::Pointer m_SelectedImage; mitk::Image::Pointer m_InternalImage; mitk::PlanarFigure::Pointer m_PlanarFigure; ItkUCharImageType::Pointer m_InternalImageMask3D; ItkUCharImageType::Pointer m_PlanarFigureImage; float m_UpsamplingFactor; ///< upsampling factor for all image generations mitk::DataNode::Pointer m_RoiImageNode; unsigned int m_StartInteractionObserverTag; unsigned int m_EndInteractionObserverTag; mitk::DataNode::Pointer m_InteractiveNode; void AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode=nullptr); void debugPFComposition(mitk::PlanarFigureComposite::Pointer , int ); void WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image); mitk::DataNode::Pointer GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute); mitk::DataNode::Pointer GenerateColorHeatmap(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib); void NodeAdded( const mitk::DataNode* node ) override; void NodeRemoved(const mitk::DataNode* node) override; void RemoveObservers(); void AddObservers(); }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui index 8a38651dfc..e56f21b27a 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingViewControls.ui @@ -1,1603 +1,1732 @@ QmitkFiberProcessingViewControls 0 0 385 - 684 + 711 Form QCommandLinkButton:disabled { border: none; } QGroupBox { background-color: transparent; } 9 9 9 9 75 true 0 5 0 - 0 + -18 353 - 441 + 503 Fiber Extraction Extract a fiber subset from the selected fiber bundle using manually placed planar figures as waypoints or binary regions of interest. - + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + 0 + + + 0 + + + + + Min. overlap: + + + + + + + Threshold: + + + + + + + Min. num. fibers: + + + + + + + Threshold on ROI image for positions to be considered as positive. + + + 3 + + + 9999.000000000000000 + + + 0.100000000000000 + + + 0.500000000000000 + + + + + + + + 0 + 0 + + + + + Ending in ROI + + + + + Not ending in ROI + + + + + Passing ROI + + + + + Not passing ROI + + + + + Labelmap + + + + + + + + Both ends + + + true + + + + + + + Interpolate ROI + + + true + + + + + + + 1 + + + 999999999 + + + + + + + Extract fibers: + + + + + + + Minimum overlap of streamlines and ROI in terms of streamline length. Zero means that one streamline point inside the ROI is enough to be considered as "overlapping". + + + 3 + + + 1.000000000000000 + + + 0.100000000000000 + + + + + + + Labels: + + + + + + + Label values seperated by whitespace. + + + ALL + + + + + + + false 0 0 200 16777215 11 Extract fibers passing through selected ROI or composite ROI. Select ROI and fiber bundle to execute. Extract - - - - Qt::Vertical - - - - 20 - 40 - + + + + + 0 + 0 + - + + + Extract using planar figures + + + + + Extract using ROI image + + + QFrame::NoFrame QFrame::Raised 0 0 0 0 6 Interactive Extraction 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 0 0 0 30 30 Draw circular ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/circle.png:/QmitkDiffusionImaging/circle.png 32 32 false true Qt::Horizontal 40 20 30 30 Draw polygonal ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/polygon.png:/QmitkDiffusionImaging/polygon.png 32 32 true true 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 0 0 0 false 60 16777215 Create NOT composition from selected ROI. NOT false 60 16777215 Create OR composition with selected ROIs. OR Qt::Horizontal 40 20 false 60 16777215 Create AND composition with selected ROIs. AND false 0 0 16777215 16777215 11 Generate a binary image containing all selected ROIs. Select at least one ROI (planar figure) and a reference fiber bundle or image. Generate ROI Image - - - - - 0 - 0 - - - - - Extract using planar figures - - - - - Extract using ROI image - - - - - - + + + + + + 0 + -34 + 353 + 456 + + + + Fiber Removal + + + Remove fibers that satisfy certain criteria from the selected bundle. + + + + QFrame::NoFrame QFrame::Raised - + 0 0 0 0 - - + + 0 + + + - Min. overlap: + X: - - + + - Extract fibers: + Y: - - + + + + + + + + - Both ends + Z: - - true + + + + + + + + + Angle: - - + + - Minimum overlap of streamlines and ROI in terms of streamline length. Zero means that one streamline point inside the ROI is enough to be considered as "overlapping". + Angular deviation threshold in degree - 3 + 1 - 1.000000000000000 + 90.000000000000000 - 0.100000000000000 - - - - - - - - 0 - 0 - - - - - Ending in ROI - - - - - Not ending in ROI - - - - - Passing ROI - - - - - Not passing ROI - - - - - - - - Interpolate ROI + 1.000000000000000 - - true + + 25.000000000000000 - - + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + 0 + + + 0 + + + 0 + + + - Threshold: + Weight threshold: - - + + - Threshold on ROI image for positions to be considered as positive. + Only fibers with weight larger than this threshold are kept. - 3 + 5 - 9999.000000000000000 + 99999.000000000000000 0.100000000000000 - - 0.500000000000000 - - - - - - - 0 - 0 - 367 - 408 - - - - Fiber Removal - - - Remove fibers that satisfy certain criteria from the selected bundle. - - + + + + false + + + + 0 + 0 + + + + + 200 + 16777215 + + + + + 11 + + + + + + + Remove + + + + + + + + 0 + 0 + + + + + Remove fibers in direction + + + + + Remove fibers by length + + + + + Remove fibers by curvature + + + + + Remove fiber parts outside mask + + + + + Remove fiber parts inside mask + + + + + Remove fibers by weight + + + + + Remove fibers by tract density + + + + QFrame::NoFrame QFrame::Raised 0 0 0 0 If unchecked, the fiber exceeding the threshold will be split in two instead of removed. Remove Fiber false QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Max. Angular Deviation: Qt::Horizontal 40 20 Maximum angular deviation in degree 180.000000000000000 0.100000000000000 30.000000000000000 Distance: Distance in mm 1 999.000000000000000 1.000000000000000 10.000000000000000 - - - - QFrame::NoFrame + + + + Qt::Vertical - - QFrame::Raised + + + 20 + 40 + - - - 0 - - - 0 - - - 0 - - - 0 - - - 0 - - - - - X: - - - - - - - Y: - - - - - - - - - - - - - Z: - - - - - - - - - - Angle: - - - - - - - Angular deviation threshold in degree - - - 1 - - - 90.000000000000000 - - - 1.000000000000000 - - - 25.000000000000000 - - - - - + QFrame::NoFrame QFrame::Raised 0 0 0 0 Qt::Horizontal 40 20 Minimum fiber length in mm 0 999999999 20 Max. Length: Min. Length: Maximum fiber length in mm 0 999999999 300 - - - false - - - - 0 - 0 - - - - - 200 - 16777215 - - - - - 11 - - - - - - - Remove - - - - - - - Qt::Vertical - - - - 20 - 40 - - - - - - - - - 0 - 0 - - - - - Remove fibers in direction - - - - - Remove fibers by length - - - - - Remove fibers by curvature - - - - - Remove fiber parts outside mask - - - - - Remove fiber parts inside mask - - - - - Remove fibers by weight - - - - - - + QFrame::NoFrame QFrame::Raised - + 0 0 0 0 - + 0 + + 6 + + + + + + + + 5 + + + 0.000000000000000 + + + 10.000000000000000 + + + 0.100000000000000 + + + 0.010000000000000 + + + - + - Weight threshold: + Density threshold: - - + + + + Min. overlap: + + + + + - Only fibers with weight larger than this threshold are kept. + Tracts have to spend at least this fraction of their length inside the specified density region. - 5 + 3 - 99999.000000000000000 + 1.000000000000000 0.100000000000000 + + 0.500000000000000 + 0 0 367 408 Bundle Modification Modify the selected bundle with operations such as fiber resampling, FA coloring, etc. QFrame::NoFrame QFrame::Raised 0 0 0 0 0 6 Error threshold in mm: 999999999.000000000000000 0.100000000000000 0.100000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 6 Sagittal Coronal Axial Select direction: QFrame::NoFrame QFrame::Raised 0 0 0 0 0 6 If checked, the image values are not only used to color the fibers but are also used as opaxity values. Values as opacity false Scalar map: The values used to color the fibers are min-max normalized. If not checked, the values should be between 0 and 1. Normalize values true 0 0 Resample fibers (spline) Resample fibers (linear) Compress fibers Color fibers by scalar map (e.g. FA) Mirror fibers Weight bundle Color fibers by curvature Color fibers by fiber weights Color fibers by length QFrame::NoFrame QFrame::Raised 0 0 0 0 0 6 0.010000000000000 999999999.000000000000000 0.100000000000000 1.000000000000000 Point distance in mm: Qt::Vertical 20 40 false 0 0 200 16777215 11 Execute QFrame::NoFrame QFrame::Raised 0 0 0 0 0 Weight: 7 999999999.000000000000000 0.100000000000000 1.000000000000000 0 0 367 - 168 + 182 Bundle Operations Join, subtract or copy bundles. false 0 0 200 16777215 11 Returns all fibers contained in bundle X that are not contained in bundle Y (not commutative!). Select at least two fiber bundles to execute. Substract Qt::Vertical 20 40 false 0 0 200 16777215 11 Merge selected fiber bundles. Select at least two fiber bundles to execute. Join false 0 0 200 16777215 11 Merge selected fiber bundles. Select at least two fiber bundles to execute. Copy Please Select Input Data 6 6 6 6 <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true <html><head/><body><p><span style=" color:#969696;">needed for extraction</span></p></body></html> true Input DTI Fiber Bundle: Binary seed ROI. If not specified, the whole image area is seeded. ROI: Qt::Vertical 20 40 QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h