diff --git a/Modules/DiffusionCmdApps/FiberProcessing/CMakeLists.txt b/Modules/DiffusionCmdApps/FiberProcessing/CMakeLists.txt index 7cbb616..4f5a308 100644 --- a/Modules/DiffusionCmdApps/FiberProcessing/CMakeLists.txt +++ b/Modules/DiffusionCmdApps/FiberProcessing/CMakeLists.txt @@ -1,49 +1,49 @@ 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^^ Sift2WeightCopy^^ FiberExtraction^^ FiberExtractionRoi^^ FiberProcessing^^ FitFibersToImage^^ FiberDirectionExtraction^^ FiberJoin^^ FiberClustering^^ TractDensityFilter^^ PrintFiberStatistics^^ - CalculateDistanceToSegmentation^^ + DistanceToSegmentation^^ ) 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 MitkDiffusionCmdApps MitkMriSimulation MitkFiberTracking ${dependencies_list} PACKAGE_DEPENDS ) endforeach() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionCmdApps/FiberProcessing/CalculateDistanceToSegmentation.cpp b/Modules/DiffusionCmdApps/FiberProcessing/CalculateDistanceToSegmentation.cpp deleted file mode 100644 index ded8d8b..0000000 --- a/Modules/DiffusionCmdApps/FiberProcessing/CalculateDistanceToSegmentation.cpp +++ /dev/null @@ -1,90 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center. - -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 "mitkDiffusionCommandLineParser.h" -#include -#include -#include -#include -#include - -typedef itk::Image ItkFloatImgType; - -/*! -\brief -*/ -int main(int argc, char* argv[]) -{ - mitkDiffusionCommandLineParser parser; - - parser.setTitle("Calcualte distance between segmentation mesh and fibers in form of a mask or TDI"); - parser.setCategory("Fiber Tracking and Processing Methods"); - parser.setContributor("MIC"); - - parser.setArgumentPrefix("--", "-"); - parser.addArgument("", "t", mitkDiffusionCommandLineParser::String, "TDI:", "input tract density image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); - parser.addArgument("", "s", mitkDiffusionCommandLineParser::String, "Segmentation:", "input segmentation mesh (.vtp)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); - parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); - - std::map parsedArgs = parser.parseArguments(argc, argv); - if (parsedArgs.size()==0) - return EXIT_FAILURE; - - std::string inTDI = us::any_cast(parsedArgs["t"]); - std::string inSeg = us::any_cast(parsedArgs["s"]); - std::string outImageFile = us::any_cast(parsedArgs["o"]); - - try - { - mitk::Image::Pointer inputTDI= mitk::IOUtil::Load(inTDI); - - ItkFloatImgType::Pointer inputItkTDI = ItkFloatImgType::New(); - mitk::CastToItkImage(inputTDI, inputItkTDI); - - mitk::Surface::Pointer inputSeg = mitk::IOUtil::Load(inSeg); - - typedef itk::DistanceFromSegmentationImageFilter< float > ImageGeneratorType; - ImageGeneratorType::Pointer filter = ImageGeneratorType::New(); - filter->SetInput(inputItkTDI); - filter->SetSegmentationSurface(inputSeg); - filter->Update(); - auto outImg = filter->GetOutput(); - - // get output image - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk(outImg); - img->SetVolume(outImg->GetBufferPointer()); - - mitk::IOUtil::Save(img, outImageFile ); - } - catch (const itk::ExceptionObject& e) - { - std::cout << e.what(); - 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/DiffusionCmdApps/FiberProcessing/DistanceToSegmentation.cpp b/Modules/DiffusionCmdApps/FiberProcessing/DistanceToSegmentation.cpp new file mode 100644 index 0000000..72aa59c --- /dev/null +++ b/Modules/DiffusionCmdApps/FiberProcessing/DistanceToSegmentation.cpp @@ -0,0 +1,148 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center. + +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 "mitkDiffusionCommandLineParser.h" +#include +#include +#include +#include +#include +#include +#include + +typedef itk::Image ItkFloatImgType; + +/*! +\brief +*/ +int main(int argc, char* argv[]) +{ + mitkDiffusionCommandLineParser parser; + + parser.setTitle("Calcualte distance between segmentation mesh and fibers in form of a mask or TDI"); + parser.setCategory("Fiber Tracking and Processing Methods"); + parser.setContributor("MIC"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("tracts", "", mitkDiffusionCommandLineParser::String, "Tracts:", "input fiber tracts (as tractogram or as tract density iumage (TDI))", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); + parser.addArgument("surface", "", mitkDiffusionCommandLineParser::String, "Segmentation:", "input segmentation mesh (.vtp)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); + parser.addArgument("out_file", "", mitkDiffusionCommandLineParser::String, "Output:", "output text file", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); + parser.addArgument("thresholds", "", mitkDiffusionCommandLineParser::StringList, "Thresholds:", "distances for which voxels are counted", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); + + std::map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + std::string in_tracts = us::any_cast(parsedArgs["tracts"]); + std::string in_seg = us::any_cast(parsedArgs["surface"]); + std::string outImageFile = us::any_cast(parsedArgs["out_file"]); + + std::vector< std::string > thresholds_str = {"0.0", "3.0", "5.0", "7.0"}; + std::vector< float > thresholds; + if (parsedArgs.count("thresholds")) + thresholds_str = us::any_cast(parsedArgs["thresholds"]); + for (auto s : thresholds_str) + thresholds.push_back(boost::lexical_cast(s)); + + try + { + ItkFloatImgType::Pointer inputItkTDI; + + auto input= mitk::IOUtil::Load(in_tracts); + + if (dynamic_cast(input.GetPointer())) + { + mitk::CastToItkImage(dynamic_cast(input.GetPointer()), inputItkTDI); + } + else + { + mitk::FiberBundle::Pointer intput_tracts= dynamic_cast(input.GetPointer()); + + itk::TractDensityImageFilter< ItkFloatImgType >::Pointer tdi_filter = itk::TractDensityImageFilter< ItkFloatImgType >::New(); + tdi_filter->SetFiberBundle(intput_tracts); + tdi_filter->SetBinaryOutput(false); + tdi_filter->SetOutputAbsoluteValues(false); + tdi_filter->Update(); + + inputItkTDI = tdi_filter->GetOutput(); + } + + mitk::Surface::Pointer inputSeg = mitk::IOUtil::Load(in_seg); + + typedef itk::DistanceFromSegmentationImageFilter< float > ImageGeneratorType; + ImageGeneratorType::Pointer filter = ImageGeneratorType::New(); + filter->SetInput(inputItkTDI); + filter->SetSegmentationSurface(inputSeg); + filter->SetThresholds(thresholds); + filter->Update(); +// auto outImg = filter->GetOutput(); + + auto volume = inputItkTDI->GetSpacing()[0]*inputItkTDI->GetSpacing()[1]*inputItkTDI->GetSpacing()[2]; + + + + std::ofstream statistics_file; + statistics_file.open(outImageFile, std::ios_base::out | std::ios_base::app); + + + statistics_file << "MITK Diffusion git commit hash: " << MITKDIFFUSION_REVISION << std::endl; + statistics_file << "MITK Diffusion branch name: " << MITKDIFFUSION_REVISION_NAME << std::endl; + statistics_file << "MITK git commit hash: " << MITK_REVISION << std::endl; + statistics_file << "MITK branch name: " << MITK_REVISION_NAME << std::endl; + + statistics_file << "Tractogram file: " << in_tracts << std::endl; + statistics_file << "Segmentation surface file: " << in_seg << std::endl; + + MITK_INFO << "Reference image spacing: " << inputItkTDI->GetSpacing(); + statistics_file << "Reference image spacing: " << inputItkTDI->GetSpacing() << std::endl; + + MITK_INFO << "Distance between surface and closest fiber-containing voxel-center: " << filter->GetMinDistance() << " mm"; + statistics_file << "Distance between surface and closest fiber-containing voxel-center: " << filter->GetMinDistance() << " mm" << std::endl; + + for (unsigned int i=0; iGetThresholds().size(); ++i) + { + MITK_INFO << "Voxels closer than " << filter->GetThresholds().at(i) << " mm: " << filter->GetCounts().at(i) << " (" << filter->GetCounts().at(i)*volume << " mm³)"; + statistics_file << "Voxels closer than " << filter->GetThresholds().at(i) << " mm: " << filter->GetCounts().at(i) << " (" << filter->GetCounts().at(i)*volume << " mm³)" << std::endl; + } + + statistics_file << "-------------------------------------------------\n\n" << std::endl; + statistics_file.close(); + +// // get output image +// mitk::Image::Pointer img = mitk::Image::New(); +// img->InitializeByItk(outImg); +// img->SetVolume(outImg->GetBufferPointer()); + +// mitk::IOUtil::Save(img, outImageFile ); + } + catch (const itk::ExceptionObject& e) + { + std::cout << e.what(); + 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/FiberTracking/Algorithms/itkDistanceFromSegmentationImageFilter.cpp b/Modules/FiberTracking/Algorithms/itkDistanceFromSegmentationImageFilter.cpp index 410f21f..56ffc3e 100644 --- a/Modules/FiberTracking/Algorithms/itkDistanceFromSegmentationImageFilter.cpp +++ b/Modules/FiberTracking/Algorithms/itkDistanceFromSegmentationImageFilter.cpp @@ -1,78 +1,107 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. 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 __itkDistanceFromSegmentationImageFilter_txx #define __itkDistanceFromSegmentationImageFilter_txx #include "itkDistanceFromSegmentationImageFilter.h" #include #include #include #include #include namespace itk { template< class TPixelType > DistanceFromSegmentationImageFilter< TPixelType >::DistanceFromSegmentationImageFilter() { - + m_Thresholds = {0.0, 3.0, 5.0, 7.0}; } template< class TPixelType > void DistanceFromSegmentationImageFilter< TPixelType >::GenerateData() { typename InputImageType::Pointer input_tdi = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); ImageRegionConstIterator< InputImageType > tdi_it(input_tdi, input_tdi->GetLargestPossibleRegion()); typename OutputImageType::Pointer outputImage = this->GetOutput(); outputImage->SetOrigin( input_tdi->GetOrigin() ); outputImage->SetRegions( input_tdi->GetLargestPossibleRegion() ); outputImage->SetSpacing( input_tdi->GetSpacing() ); outputImage->SetDirection( input_tdi->GetDirection() ); outputImage->Allocate(); outputImage->FillBuffer(0.0); ImageRegionIterator< InputImageType > out_it(outputImage, outputImage->GetLargestPossibleRegion()); vtkSmartPointer vtkFilter = vtkSmartPointer::New(); vtkFilter->SetInput(m_SegmentationSurface->GetVtkPolyData()); - float min_distance = 999999; + m_MinDistance = 999999; + m_Counts.resize(m_Thresholds.size(), 0); while( !tdi_it.IsAtEnd() ) { TPixelType tdi_val = tdi_it.Get(); if (tdi_val>0) { auto idx = tdi_it.GetIndex(); itk::Point point3D; input_tdi->TransformIndexToPhysicalPoint(idx, point3D); double dist = vtkFilter->EvaluateFunction(point3D[0], point3D[1], point3D[2]); out_it.Set(dist); - if (dist +float DistanceFromSegmentationImageFilter< TPixelType >::GetMinDistance() const +{ + return m_MinDistance; +} + +template< class TPixelType > +std::vector DistanceFromSegmentationImageFilter< TPixelType >::GetThresholds() const +{ + return m_Thresholds; +} + +template< class TPixelType > +std::vector DistanceFromSegmentationImageFilter< TPixelType >::GetCounts() const +{ + return m_Counts; +} + +template< class TPixelType > +void DistanceFromSegmentationImageFilter< TPixelType >::SetThresholds(const std::vector &Thresholds) +{ + m_Thresholds = Thresholds; } } #endif diff --git a/Modules/FiberTracking/Algorithms/itkDistanceFromSegmentationImageFilter.h b/Modules/FiberTracking/Algorithms/itkDistanceFromSegmentationImageFilter.h index 35b40d6..f489c1f 100644 --- a/Modules/FiberTracking/Algorithms/itkDistanceFromSegmentationImageFilter.h +++ b/Modules/FiberTracking/Algorithms/itkDistanceFromSegmentationImageFilter.h @@ -1,76 +1,86 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. 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. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkDistanceFromSegmentationImageFilter_h_ #define __itkDistanceFromSegmentationImageFilter_h_ #include #include #include namespace itk{ /** * \brief */ template< class TPixelType > class DistanceFromSegmentationImageFilter : public ImageToImageFilter< Image< TPixelType, 3 >, Image< TPixelType, 3 > > { public: typedef DistanceFromSegmentationImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< TPixelType, 3 >, Image< TPixelType, 3 > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(DistanceFromSegmentationImageFilter, ImageToImageFilter) itkSetMacro( SegmentationSurface, mitk::Surface::Pointer ) typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + void SetThresholds(const std::vector &Thresholds); + std::vector GetThresholds() const; + + std::vector GetCounts() const; + + float GetMinDistance() const; + protected: DistanceFromSegmentationImageFilter(); ~DistanceFromSegmentationImageFilter() override {} void GenerateData() override; private: mitk::Surface::Pointer m_SegmentationSurface; + std::vector m_Thresholds; + std::vector m_Counts; + float m_MinDistance; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDistanceFromSegmentationImageFilter.cpp" #endif #endif //__itkDistanceFromSegmentationImageFilter_h_