diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt index 11f0c2aa2b..7fbf48d40b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCore/cmdapps/CMakeLists.txt @@ -1,44 +1,44 @@ option(BUILD_DiffusionMiniApps "Build commandline tools for diffusion" OFF) if(BUILD_DiffusionMiniApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion miniapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionminiapps DwiDenoising^^ ImageResampler^^ ExportShImage^^ CopyGeometry^^ DiffusionIndices^^ QballReconstruction^^ Registration^^ TensorReconstruction^^ TensorDerivedMapsExtraction^^ DiffusionDICOMLoader^^ - DiffusionIVIMFit^^ + DiffusionKurtosisFit^^ ) foreach(diffusionminiapp ${diffusionminiapps}) # extract mini app name and dependencies string(REPLACE "^^" "\\;" miniapp_info ${diffusionminiapp}) set(miniapp_info_list ${miniapp_info}) list(GET miniapp_info_list 0 appname) list(GET miniapp_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 ) endforeach() endif() diff --git a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionIVIMFit.cpp b/Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionKurtosisFit.cpp similarity index 60% rename from Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionIVIMFit.cpp rename to Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionKurtosisFit.cpp index 6c49d7a7d0..d1886d0469 100644 --- a/Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionIVIMFit.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/cmdapps/DiffusionKurtosisFit.cpp @@ -1,207 +1,253 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include +#include "mitkImage.h" +#include +#include +#include +#include "mitkIOUtil.h" + +#include +#include + + //vnl_includes #include "vnl/vnl_math.h" #include "vnl/vnl_cost_function.h" #include "vnl/vnl_least_squares_function.h" #include "vnl/algo/vnl_lbfgsb.h" #include "vnl/algo/vnl_lbfgs.h" #include "vnl/algo/vnl_levenberg_marquardt.h" typedef mitk::DiffusionPropertyHelper DPH; #include #include #include #include #include #include + + DPH::ImageType::Pointer GetBlurredVectorImage( DPH::ImageType::Pointer vectorImage, double sigma) { typedef itk::DiscreteGaussianImageFilter< itk::Image, itk::Image > GaussianFilterType; typedef itk::VectorIndexSelectionCastImageFilter< DPH::ImageType, itk::Image > IndexSelectionType; IndexSelectionType::Pointer indexSelectionFilter = IndexSelectionType::New(); indexSelectionFilter->SetInput( vectorImage ); typedef itk::ComposeImageFilter< itk::Image, DPH::ImageType > ComposeFilterType; ComposeFilterType::Pointer vec_composer = ComposeFilterType::New(); for( unsigned int i=0; iGetVectorLength(); ++i) { GaussianFilterType::Pointer gaussian_filter = GaussianFilterType::New(); indexSelectionFilter->SetIndex( i ); gaussian_filter->SetInput( indexSelectionFilter->GetOutput() ); gaussian_filter->SetVariance( sigma ); vec_composer->SetInput(i, gaussian_filter->GetOutput() ); gaussian_filter->Update(); } try { vec_composer->Update(); } catch(const itk::ExceptionObject &e) { mitkThrow() << "[VectorImage.GaussianSmoothing] !! Failed with ITK Exception: " << e.what(); } DPH::ImageType::Pointer smoothed_vector = vec_composer->GetOutput(); -/* + + /* itk::ImageFileWriter< DPH::ImageType >::Pointer writer = itk::ImageFileWriter< DPH::ImageType >::New(); writer->SetInput( smoothed_vector ); writer->SetFileName( "/tmp/itk_smoothed_vector.nrrd"); writer->Update();*/ return smoothed_vector; } -void KurtosisMapComputation( mitk::Image::Pointer input, std::string output_prefix ) +void KurtosisMapComputation( mitk::Image::Pointer input, + std::string output_prefix , + std::string output_type, + std::string maskPath, + bool omitBZero, + double lower, + double upper ) { DPH::ImageType::Pointer vectorImage = DPH::ImageType::New(); mitk::CastToItkImage( input, vectorImage ); - - typedef itk::DiffusionKurtosisReconstructionImageFilter< short, double > KurtosisFilterType; KurtosisFilterType::Pointer kurtosis_filter = KurtosisFilterType::New(); - kurtosis_filter->SetInput( GetBlurredVectorImage( vectorImage, 1.5 ) ); kurtosis_filter->SetReferenceBValue( DPH::GetReferenceBValue( input.GetPointer() ) ); kurtosis_filter->SetGradientDirections( DPH::GetGradientContainer( input.GetPointer() ) ); - kurtosis_filter->SetNumberOfThreads(1); +// kurtosis_filter->SetNumberOfThreads(1); + kurtosis_filter->SetOmitUnweightedValue(omitBZero); + kurtosis_filter->SetBoundariesForKurtosis(-lower,upper); +// kurtosis_filter->SetInitialSolution(const vnl_vector& x0 ); - KurtosisFilterType::OutputImageRegionType o_region; - KurtosisFilterType::OutputImageRegionType::SizeType o_size; - KurtosisFilterType::OutputImageRegionType::IndexType o_index; - o_index.Fill(0); o_size.Fill(0); - o_index[0] = 48; o_index[1] = 18; o_index[2] = 12; - o_size[0] = 16; o_size[1] = 24; o_size[2] = 1; - - o_region.SetSize( o_size ); - o_region.SetIndex( o_index ); - kurtosis_filter->SetMapOutputRegion( o_region ); + if(maskPath != "") + { + mitk::Image::Pointer segmentation; + segmentation = mitk::IOUtil::LoadImage(maskPath); + typedef itk::Image< short , 3> MaskImageType; + MaskImageType::Pointer vectorSeg = MaskImageType::New() ; + mitk::CastToItkImage( segmentation, vectorSeg ); + kurtosis_filter->SetImageMask(vectorSeg) ; + } try { kurtosis_filter->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Kurtosis fit failed with an ITK Exception: " << e.what(); } mitk::Image::Pointer d_image = mitk::Image::New(); d_image->InitializeByItk( kurtosis_filter->GetOutput(0) ); d_image->SetVolume( kurtosis_filter->GetOutput(0)->GetBufferPointer() ); mitk::Image::Pointer k_image = mitk::Image::New(); k_image->InitializeByItk( kurtosis_filter->GetOutput(1) ); k_image->SetVolume( kurtosis_filter->GetOutput(1)->GetBufferPointer() ); - std::string outputD_FileName = output_prefix + "_ADC_map.nrrd"; - std::string outputK_FileName = output_prefix + "_AKC_map.nrrd"; + std::string outputD_FileName = output_prefix + "_ADC_map." + output_type; + std::string outputK_FileName = output_prefix + "_AKC_map." + output_type; try { - mitk::IOUtil::Save( d_image, outputD_FileName ); - mitk::IOUtil::Save( k_image, outputK_FileName ); + mitk::IOUtil::Save( d_image, outputD_FileName ); + mitk::IOUtil::Save( k_image, outputK_FileName ); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Failed to save the KurtosisFit Results due to exception: " << e.what(); } } int main( int argc, char* argv[] ) { - mitkCommandLineParser parser; - parser.setTitle("Diffusion IVIM (Kurtosis) Fit"); + parser.setTitle("Diffusion Kurtosis Fit"); parser.setCategory("Signal Reconstruction"); parser.setContributor("MIC"); - parser.setDescription("Fitting of IVIM / Kurtosis"); - + parser.setDescription("Fitting Kurtosis"); parser.setArgumentPrefix("--","-"); + // mandatory arguments parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input: ", "input image (DWI)", us::Any(), false); parser.addArgument("output", "o", mitkCommandLineParser::String, "Output Preifx: ", "Prefix for the output images, will append _f, _K, _D accordingly ", us::Any(), false); - parser.addArgument("fit", "f", mitkCommandLineParser::String, "Input: ", "input image (DWI)", us::Any(), false); + parser.addArgument("output_type", "ot", mitkCommandLineParser::String, "Output Type: ", "choose data type of output image, e.g. '.nii' or '.nrrd' ", us::Any(), false); // optional arguments - parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Masking Image: ", "ROI (segmentation)", us::Any(), true); + parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Masking Image: ", "ROI (segmentation)", us::Any()); + parser.addArgument("help", "h", mitkCommandLineParser::Bool, "Help", "Show this help text"); + parser.addArgument("omitbzero", "om", mitkCommandLineParser::Bool, "Omit b0:", "Omit b0 value during fit (default = false)", us::Any()); + parser.addArgument("lowerkbound", "kl", mitkCommandLineParser::Float, "lower Kbound:", "Set (unsigned) lower boundary for Kurtosis parameter (default = -1000)", us::Any()); + parser.addArgument("upperkbound", "ku", mitkCommandLineParser::Float, "upper Kbound:", "Set upper boundary for Kurtosis parameter (default = 1000)", us::Any()); std::map parsedArgs = parser.parseArguments(argc, argv); - if (parsedArgs.size()==0) - return EXIT_FAILURE; + + if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")){ + std::cout << parser.helpText(); + return EXIT_SUCCESS; + } // mandatory arguments std::string inFileName = us::any_cast(parsedArgs["input"]); std::string out_prefix = us::any_cast(parsedArgs["output"]); - std::string fit_name = us::any_cast(parsedArgs["fit"]); + std::string maskPath = ""; mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inFileName); - if( !DPH::IsDiffusionWeightedImage( inputImage ) ) + + bool omitBZero = false; + double lower = -1000; + double upper = 1000; + std::string out_type = "nrrd"; + + if (parsedArgs.count("mask") || parsedArgs.count("m")) { - MITK_ERROR("DiffusionIVIMFit.Input") << "No valid diffusion-weighted image provided, failed to load " << inFileName << " as DW Image. Aborting..."; - return EXIT_FAILURE; + maskPath = us::any_cast(parsedArgs["mask"]); } - if( fit_name == "Kurtosis" ) + if (parsedArgs.count("output_type") || parsedArgs.count("ot")) { - MITK_INFO("DiffusionIVIMFit.Main") << "-----[ Kurtosis Fit ]-----"; + out_type = us::any_cast(parsedArgs["output_type"]); + } - KurtosisMapComputation( inputImage, out_prefix ); + if (parsedArgs.count("omitbzero") || parsedArgs.count("om")) + { + omitBZero = us::any_cast(parsedArgs["omitbzero"]); + } + if (parsedArgs.count("lowerkbound") || parsedArgs.count("kl")) + { + lower = us::any_cast(parsedArgs["lowerkbound"]); } - else if (fit_name == "IVIM" ) + + if (parsedArgs.count("upperkbound") || parsedArgs.count("ku")) { - MITK_INFO("DiffusionIVIMFit.Main") << "IVIM Fit not fully implemented yet. Aborting..."; - return EXIT_FAILURE; + upper = us::any_cast(parsedArgs["upperkbound"]); } - else + + if( !DPH::IsDiffusionWeightedImage( inputImage ) ) { - MITK_ERROR("DiffusionIVIMFit.Main") << "Unrecognized option: " << fit_name << ". Valid values [\"IVIM\", \"Kurtosis\"] \n Aborting... \n"; + MITK_ERROR("DiffusionIVIMFit.Input") << "No valid diffusion-weighted image provided, failed to load " << inFileName << " as DW Image. Aborting..."; return EXIT_FAILURE; - } + + +KurtosisMapComputation( inputImage, + out_prefix , + out_type, + maskPath, + omitBZero, + lower, + upper); + }