diff --git a/Modules/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp b/Modules/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp index 98235b3..da7a692 100644 --- a/Modules/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp +++ b/Modules/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp @@ -1,321 +1,321 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include #include #include typedef itk::Point PointType4; typedef itk::Image< float, 4 > PeakImgType; /*! \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[]) { mitkDiffusionCommandLineParser 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", mitkDiffusionCommandLineParser::StringList, "Input tractograms:", "input tractograms (files or folder)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "i2", mitkDiffusionCommandLineParser::String, "Input image:", "input image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output root", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("max_iter", "", mitkDiffusionCommandLineParser::Int, "Max. iterations:", "maximum number of optimizer iterations", 20); parser.addArgument("bundle_based", "", mitkDiffusionCommandLineParser::Bool, "Bundle based fit:", "fit one weight per input tractogram/bundle, not for each fiber", false); parser.addArgument("min_g", "", mitkDiffusionCommandLineParser::Float, "Min. g:", "lower termination threshold for gradient magnitude", 1e-5); parser.addArgument("lambda", "", mitkDiffusionCommandLineParser::Float, "Lambda:", "modifier for regularization", 1.0); parser.addArgument("save_res", "", mitkDiffusionCommandLineParser::Bool, "Save Residuals:", "save residual images", false); parser.addArgument("save_weights", "", mitkDiffusionCommandLineParser::Bool, "Save Weights:", "save fiber weights in a separate text file", false); parser.addArgument("filter_zero", "", mitkDiffusionCommandLineParser::Bool, "Filter Zero Weights:", "filter fibers with zero weight", false); parser.addArgument("filter_outliers", "", mitkDiffusionCommandLineParser::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", "", mitkDiffusionCommandLineParser::Bool, "Join output tracts:", "outout tracts are merged into a single tractogram", false); parser.addArgument("regu", "", mitkDiffusionCommandLineParser::String, "Regularization:", "MSM; Variance; VoxelVariance; Lasso; GroupLasso; GroupVariance; NONE", std::string("VoxelVariance")); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::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 filter_zero = false; if (parsedArgs.count("filter_zero")) filter_zero = us::any_cast(parsedArgs["filter_zero"]); 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 = 1.0; 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"; - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, std::vector()); std::vector< std::string > fib_names; 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(); mitk::LocaleSwitch localeSwitch("C"); 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; bundleFilterByWeights(0.0); if (fib->GetNumFibers()>0) { fib->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(fib, outRoot + name + "_fitted.fib"); if (save_weights) { ofstream logfile; logfile.open (outRoot + name + "_weights.txt"); for (unsigned int f=0; fGetNumFibers(); ++f) logfile << output_tracts.at(bundle)->GetFiberWeight(f) << "\n"; logfile.close(); } } else MITK_INFO << "Output contains no fibers!"; } } else { mitk::FiberBundle::Pointer out = mitk::FiberBundle::New(); out = out->AddBundles(output_tracts); if (filter_zero) out = out->FilterByWeights(0.0); if (out->GetNumFibers()>0) { out->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(out, outRoot + "_fitted.fib"); if (save_weights) { ofstream logfile; logfile.open (outRoot + "_weights.txt"); for (unsigned int f=0; fGetNumFibers(); ++f) logfile << out->GetFiberWeight(f) << "\n"; logfile.close(); } } else MITK_INFO << "Output contains no fibers!"; } } 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/Fiberfox/Fiberfox.cpp b/Modules/DiffusionCmdApps/Fiberfox/Fiberfox.cpp index cd9beb9..b9fc4d5 100644 --- a/Modules/DiffusionCmdApps/Fiberfox/Fiberfox.cpp +++ b/Modules/DiffusionCmdApps/Fiberfox/Fiberfox.cpp @@ -1,289 +1,289 @@ /*=================================================================== 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 #include #include #include #include #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include using namespace mitk; /*! * \brief Command line interface to Fiberfox. * Simulate a diffusion-weighted image from a tractogram using the specified parameter file. */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Fiberfox"); parser.setCategory("Diffusion Simulation Tools"); parser.setContributor("MIC"); parser.setDescription("Command line interface to Fiberfox." " Simulate a diffusion-weighted image from a tractogram using the specified parameter file."); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output root:", "output folder and file prefix", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input tractogram or diffusion-weighted image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("parameters", "p", mitkDiffusionCommandLineParser::String, "Parameter file:", "fiberfox parameter file (.ffp)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("template", "t", mitkDiffusionCommandLineParser::String, "Template image:", "use parameters of the template image", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("verbose", "v", mitkDiffusionCommandLineParser::Bool, "Output additional images:", "output volume fraction images etc.", us::Any()); parser.addArgument("dont_apply_direction_matrix", "", mitkDiffusionCommandLineParser::Bool, "Don't apply direction matrix:", "don't rotate gradients by image direction matrix", us::Any()); parser.addArgument("fix_seed", "", mitkDiffusionCommandLineParser::Bool, "Use fix random seed:", "always use same sequence of random numbers", us::Any()); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) { return EXIT_FAILURE; } std::string outName = us::any_cast(parsedArgs["o"]); std::string paramName = us::any_cast(parsedArgs["parameters"]); std::string input=""; if (parsedArgs.count("i")) input = us::any_cast(parsedArgs["i"]); bool fix_seed = false; if (parsedArgs.count("fix_seed")) fix_seed = us::any_cast(parsedArgs["fix_seed"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = us::any_cast(parsedArgs["verbose"]); bool apply_direction_matrix = true; if (parsedArgs.count("dont_apply_direction_matrix")) apply_direction_matrix = false; FiberfoxParameters parameters; parameters.LoadParameters(paramName, fix_seed); // Test if /path/dir is an existing directory: std::string file_extension = ""; if( itksys::SystemTools::FileIsDirectory( outName ) ) { while( *(--(outName.cend())) == '/') { outName.pop_back(); } outName = outName + '/'; parameters.m_Misc.m_OutputPath = outName; outName = outName + parameters.m_Misc.m_OutputPrefix; // using default m_OutputPrefix as initialized. } else { // outName is NOT an existing directory, so we need to remove all trailing slashes: while( *(--(outName.cend())) == '/') { outName.pop_back(); } // now split up the given outName into directory and (prefix of) filename: if( ! itksys::SystemTools::GetFilenamePath( outName ).empty() && itksys::SystemTools::FileIsDirectory(itksys::SystemTools::GetFilenamePath( outName ) ) ) { parameters.m_Misc.m_OutputPath = itksys::SystemTools::GetFilenamePath( outName ) + '/'; } else { parameters.m_Misc.m_OutputPath = itksys::SystemTools::GetCurrentWorkingDirectory() + '/'; } file_extension = itksys::SystemTools::GetFilenameExtension(outName); if( ! itksys::SystemTools::GetFilenameWithoutExtension( outName ).empty() ) { parameters.m_Misc.m_OutputPrefix = itksys::SystemTools::GetFilenameWithoutExtension( outName ); } else { parameters.m_Misc.m_OutputPrefix = "fiberfox"; } outName = parameters.m_Misc.m_OutputPath + parameters.m_Misc.m_OutputPrefix; } - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "Fiberbundles"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "Fiberbundles"}, std::vector()); mitk::BaseData::Pointer inputData = mitk::IOUtil::Load(input, &functor)[0]; itk::TractsToDWIImageFilter< short >::Pointer tractsToDwiFilter = itk::TractsToDWIImageFilter< short >::New(); if ( dynamic_cast(inputData.GetPointer()) ) // simulate dataset from fibers { tractsToDwiFilter->SetFiberBundle(dynamic_cast(inputData.GetPointer())); if (parsedArgs.count("template")) { MITK_INFO << "Loading template image"; typedef itk::VectorImage< short, 3 > ItkDwiType; typedef itk::Image< short, 3 > ItkImageType; mitk::BaseData::Pointer templateData = mitk::IOUtil::Load(us::any_cast(parsedArgs["template"]), &functor)[0]; mitk::Image::Pointer template_image = dynamic_cast(templateData.GetPointer()); if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(template_image)) { ItkDwiType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(template_image); parameters.m_SignalGen.m_ImageRegion = itkVectorImagePointer->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkVectorImagePointer->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkVectorImagePointer->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkVectorImagePointer->GetDirection(); parameters.SetBvalue(mitk::DiffusionPropertyHelper::GetReferenceBValue(template_image)); parameters.SetGradienDirections(mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(template_image)); } else { ItkImageType::Pointer itkImagePointer = ItkImageType::New(); mitk::CastToItkImage(template_image, itkImagePointer); parameters.m_SignalGen.m_ImageRegion = itkImagePointer->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkImagePointer->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkImagePointer->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkImagePointer->GetDirection(); } } } else if ( dynamic_cast(inputData.GetPointer()) ) // add artifacts to existing image { typedef itk::VectorImage< short, 3 > ItkDwiType; mitk::Image::Pointer diffImg = dynamic_cast(inputData.GetPointer()); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(diffImg, itkVectorImagePointer); parameters.m_SignalGen.m_SignalScale = 1; parameters.m_SignalGen.m_ImageRegion = itkVectorImagePointer->GetLargestPossibleRegion(); parameters.m_SignalGen.m_ImageSpacing = itkVectorImagePointer->GetSpacing(); parameters.m_SignalGen.m_ImageOrigin = itkVectorImagePointer->GetOrigin(); parameters.m_SignalGen.m_ImageDirection = itkVectorImagePointer->GetDirection(); parameters.SetBvalue(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImg)); parameters.SetGradienDirections(mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(diffImg)); tractsToDwiFilter->SetInputImage(itkVectorImagePointer); } if (verbose) { MITK_DEBUG << outName << ".ffp"; parameters.m_Misc.m_OutputAdditionalImages = true; parameters.SaveParameters(outName+".ffp"); } else parameters.m_Misc.m_OutputAdditionalImages = false; if (apply_direction_matrix) { MITK_INFO << "Applying direction matrix to gradient directions."; parameters.ApplyDirectionMatrix(); } tractsToDwiFilter->SetParameters(parameters); tractsToDwiFilter->SetUseConstantRandSeed(fix_seed); tractsToDwiFilter->Update(); mitk::Image::Pointer image = mitk::GrabItkImageMemory(tractsToDwiFilter->GetOutput()); if (parameters.m_SignalGen.GetNumWeightedVolumes()>0) { if (apply_direction_matrix) mitk::DiffusionPropertyHelper::SetGradientContainer(image, parameters.m_SignalGen.GetItkGradientContainer()); else mitk::DiffusionPropertyHelper::SetOriginalGradientContainer(image, parameters.m_SignalGen.GetItkGradientContainer()); mitk::DiffusionPropertyHelper::SetReferenceBValue(image, parameters.m_SignalGen.GetBvalue()); mitk::DiffusionPropertyHelper::InitializeImage(image); if (file_extension=="") mitk::IOUtil::Save(image, "DWI_NIFTI", outName+".nii.gz"); else if (file_extension==".nii" || file_extension==".nii.gz") mitk::IOUtil::Save(image, "DWI_NIFTI", outName+file_extension); else mitk::IOUtil::Save(image, outName+file_extension); } else mitk::IOUtil::Save(image, outName+".nii.gz"); if (verbose) { if (tractsToDwiFilter->GetTickImage().IsNotNull()) { mitk::Image::Pointer mitkImage = mitk::Image::New(); itk::TractsToDWIImageFilter< short >::Float2DImageType::Pointer itkImage = tractsToDwiFilter->GetTickImage(); mitkImage = mitk::GrabItkImageMemory( itkImage.GetPointer() ); mitk::IOUtil::Save(mitkImage, outName+"_Ticks.nii.gz"); } if (tractsToDwiFilter->GetRfImage().IsNotNull()) { mitk::Image::Pointer mitkImage = mitk::Image::New(); itk::TractsToDWIImageFilter< short >::Float2DImageType::Pointer itkImage = tractsToDwiFilter->GetRfImage(); mitkImage = mitk::GrabItkImageMemory( itkImage.GetPointer() ); mitk::IOUtil::Save(mitkImage, outName+"_TimeFromRf.nii.gz"); } std::vector< itk::TractsToDWIImageFilter< short >::ItkDoubleImgType::Pointer > volumeFractions = tractsToDwiFilter->GetVolumeFractions(); for (unsigned int k=0; kInitializeByItk(volumeFractions.at(k).GetPointer()); image->SetVolume(volumeFractions.at(k)->GetBufferPointer()); mitk::IOUtil::Save(image, outName+"_Compartment"+boost::lexical_cast(k+1)+".nii.gz"); } if (tractsToDwiFilter->GetPhaseImage().IsNotNull()) { mitk::Image::Pointer image = mitk::Image::New(); itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer itkPhase = tractsToDwiFilter->GetPhaseImage(); image = mitk::GrabItkImageMemory( itkPhase.GetPointer() ); mitk::IOUtil::Save(image, outName+"_Phase.nii.gz"); } if (tractsToDwiFilter->GetKspaceImage().IsNotNull()) { mitk::Image::Pointer image = mitk::Image::New(); itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer itkImage = tractsToDwiFilter->GetKspaceImage(); image = mitk::GrabItkImageMemory( itkImage.GetPointer() ); mitk::IOUtil::Save(image, outName+"_kSpace.nii.gz"); } int c = 1; std::vector< itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer > output_real = tractsToDwiFilter->GetOutputImagesReal(); for (auto real : output_real) { mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(real.GetPointer()); image->SetVolume(real->GetBufferPointer()); mitk::IOUtil::Save(image, outName+"_Coil-"+boost::lexical_cast(c)+"-real.nii.gz"); ++c; } c = 1; std::vector< itk::TractsToDWIImageFilter< short >::DoubleDwiType::Pointer > output_imag = tractsToDwiFilter->GetOutputImagesImag(); for (auto imag : output_imag) { mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk(imag.GetPointer()); image->SetVolume(imag->GetBufferPointer()); mitk::IOUtil::Save(image, outName+"_Coil-"+boost::lexical_cast(c)+"-imag.nii.gz"); ++c; } } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/Misc/DImp.cpp b/Modules/DiffusionCmdApps/Misc/DImp.cpp index 4bb5d25..74b9011 100644 --- a/Modules/DiffusionCmdApps/Misc/DImp.cpp +++ b/Modules/DiffusionCmdApps/Misc/DImp.cpp @@ -1,72 +1,72 @@ /*=================================================================== 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 #include #include "mitkDiffusionCommandLineParser.h" #include #include int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("DIMP"); parser.setCategory("Preprocessing Tools"); parser.setDescription("TEMPORARY: Converts DICOM to other image types"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image", 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; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); try { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); mitk::Image::Pointer source = mitk::IOUtil::Load(imageName, &functor); std::string ext = itksys::SystemTools::GetFilenameExtension(outImage); if (ext==".nii" || ext==".nii.gz") mitk::IOUtil::Save(source, "DWI_NIFTI", outImage); else mitk::IOUtil::Save(source, outImage); } 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/Misc/DReg.cpp b/Modules/DiffusionCmdApps/Misc/DReg.cpp index 83dd4c2..d1b8961 100644 --- a/Modules/DiffusionCmdApps/Misc/DReg.cpp +++ b/Modules/DiffusionCmdApps/Misc/DReg.cpp @@ -1,222 +1,222 @@ /*=================================================================== 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 #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include typedef mitk::DiffusionPropertyHelper DPH; mitk::Image::Pointer apply_transform(mitk::Image::Pointer moving, mitk::Image::Pointer fixed_single, mitk::MAPRegistrationWrapper::Pointer reg, bool resample) { mitk::Image::Pointer registered_image; if (!resample) { registered_image = mitk::ImageMappingHelper::refineGeometry(moving, reg, true); } else { if (!mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(moving)) { registered_image = mitk::ImageMappingHelper::map(moving, reg, false, 0, fixed_single->GetGeometry(), false, 0, mitk::ImageMappingInterpolator::BSpline_3); } else { typedef itk::Image ITKDiffusionVolumeType; typedef itk::ComposeImageFilter < ITKDiffusionVolumeType > ComposeFilterType; auto composer = ComposeFilterType::New(); auto itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(moving); for (unsigned int i=0; iGetVectorLength(); ++i) { itk::ExtractDwiChannelFilter< short >::Pointer filter = itk::ExtractDwiChannelFilter< short >::New(); filter->SetInput( itkVectorImagePointer); filter->SetChannelIndex(i); filter->Update(); mitk::Image::Pointer gradientVolume = mitk::Image::New(); gradientVolume->InitializeByItk( filter->GetOutput() ); gradientVolume->SetImportChannel( filter->GetOutput()->GetBufferPointer() ); mitk::Image::Pointer registered_mitk_image = mitk::ImageMappingHelper::map(gradientVolume, reg, false, 0, fixed_single->GetGeometry(), false, 0, mitk::ImageMappingInterpolator::BSpline_3); auto registered_itk_image = ITKDiffusionVolumeType::New(); mitk::CastToItkImage(registered_mitk_image, registered_itk_image); composer->SetInput(i, registered_itk_image); } composer->Update(); registered_image = mitk::GrabItkImageMemory( composer->GetOutput() ); mitk::DiffusionPropertyHelper::CopyProperties(moving, registered_image, true); typedef mitk::DiffusionImageCorrectionFilter CorrectionFilterType; CorrectionFilterType::Pointer corrector = CorrectionFilterType::New(); corrector->SetImage( registered_image ); corrector->CorrectDirections( mitk::MITKRegistrationHelper::getAffineMatrix(reg, false)->GetMatrix().GetVnlMatrix() ); } } return registered_image; } int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("DREG"); parser.setCategory("Preprocessing Tools"); parser.setDescription("TEMPORARY: Rigid registration of two images"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "f", mitkDiffusionCommandLineParser::String, "Fixed:", "fixed image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "m", mitkDiffusionCommandLineParser::String, "Moving:", "moving image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("resample", "", mitkDiffusionCommandLineParser::Bool, "Resample:", "resample moving image", false); parser.addArgument("coreg", "", mitkDiffusionCommandLineParser::StringList, "", "additionally apply transform to these images", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string f = us::any_cast(parsedArgs["f"]); std::string m = us::any_cast(parsedArgs["m"]); std::string o = us::any_cast(parsedArgs["o"]); bool resample = false; if (parsedArgs.count("resample")) resample = true; mitkDiffusionCommandLineParser::StringContainerType coreg; if (parsedArgs.count("coreg")) coreg = us::any_cast(parsedArgs["coreg"]); try { typedef itk::Image< float, 3 > ItkFloatImageType; - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); mitk::Image::Pointer fixed = mitk::IOUtil::Load(f, &functor); mitk::Image::Pointer moving = mitk::IOUtil::Load(m, &functor); mitk::Image::Pointer fixed_single = fixed; mitk::Image::Pointer moving_single = moving; mitk::MultiModalRigidDefaultRegistrationAlgorithm< ItkFloatImageType >::Pointer algo = mitk::MultiModalRigidDefaultRegistrationAlgorithm< ItkFloatImageType >::New(); mitk::MAPAlgorithmHelper helper(algo); if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(fixed)) { DPH::ImageType::Pointer itkVectorImagePointer = DPH::ImageType::New(); mitk::CastToItkImage(fixed, itkVectorImagePointer); itk::ExtractDwiChannelFilter< short >::Pointer filter = itk::ExtractDwiChannelFilter< short >::New(); filter->SetInput( itkVectorImagePointer); filter->SetChannelIndex(0); filter->Update(); fixed_single = mitk::Image::New(); fixed_single->InitializeByItk( filter->GetOutput() ); fixed_single->SetImportChannel( filter->GetOutput()->GetBufferPointer() ); } if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(moving)) { DPH::ImageType::Pointer itkVectorImagePointer = DPH::ImageType::New(); mitk::CastToItkImage(moving, itkVectorImagePointer); itk::ExtractDwiChannelFilter< short >::Pointer filter = itk::ExtractDwiChannelFilter< short >::New(); filter->SetInput( itkVectorImagePointer); filter->SetChannelIndex(0); filter->Update(); moving_single = mitk::Image::New(); moving_single->InitializeByItk( filter->GetOutput() ); moving_single->SetImportChannel( filter->GetOutput()->GetBufferPointer() ); } helper.SetData(moving_single, fixed_single); mitk::MAPRegistrationWrapper::Pointer reg = helper.GetMITKRegistrationWrapper(); mitk::Image::Pointer registered_image = apply_transform(moving, fixed_single, reg, resample); if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(registered_image)) { mitk::DiffusionPropertyHelper::InitializeImage( registered_image ); std::string file_extension = itksys::SystemTools::GetFilenameExtension(o); if (file_extension==".nii" || file_extension==".nii.gz") mitk::IOUtil::Save(registered_image, "DWI_NIFTI", o); else mitk::IOUtil::Save(registered_image, o); } else mitk::IOUtil::Save(registered_image, o); std::string path = ist::GetFilenamePath(o) + "/"; std::vector< std::string > file_names; auto coreg_images = mitk::DiffusionDataIOHelper::load_mitk_images(coreg, &file_names); for (unsigned int i=0; i #include #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("DmriDenoising"); parser.setCategory("Preprocessing Tools"); parser.setDescription("dMRI denoising tool"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("type", "", mitkDiffusionCommandLineParser::Int, "Type:", "0 (TotalVariation), 1 (Gauss), 2 (NLM)", 0); parser.endGroup(); parser.beginGroup("2. Total variation parameters:"); parser.addArgument("tv_iterations", "", mitkDiffusionCommandLineParser::Int, "Iterations:", "", 1); parser.addArgument("lambda", "", mitkDiffusionCommandLineParser::Float, "Lambda:", "", 0.1); parser.endGroup(); parser.beginGroup("3. Gauss parameters:"); parser.addArgument("variance", "", mitkDiffusionCommandLineParser::Float, "Variance:", "", 1.0); parser.endGroup(); parser.beginGroup("4. NLM parameters:"); parser.addArgument("nlm_iterations", "", mitkDiffusionCommandLineParser::Int, "Iterations:", "", 4); parser.addArgument("sampling_radius", "", mitkDiffusionCommandLineParser::Int, "Sampling radius:", "", 4); parser.addArgument("patch_radius", "", mitkDiffusionCommandLineParser::Int, "Patch radius:", "", 1); parser.addArgument("num_patches", "", mitkDiffusionCommandLineParser::Int, "Num. patches:", "", 10); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); int type = 0; if (parsedArgs.count("type")) type = us::any_cast(parsedArgs["type"]); int tv_iterations = 1; if (parsedArgs.count("tv_iterations")) tv_iterations = us::any_cast(parsedArgs["tv_iterations"]); float lambda = 0.1; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); float variance = 1.0; if (parsedArgs.count("variance")) variance = us::any_cast(parsedArgs["variance"]); int nlm_iterations = 4; if (parsedArgs.count("nlm_iterations")) nlm_iterations = us::any_cast(parsedArgs["nlm_iterations"]); int sampling_radius = 4; if (parsedArgs.count("sampling_radius")) sampling_radius = us::any_cast(parsedArgs["sampling_radius"]); int patch_radius = 1; if (parsedArgs.count("patch_radius")) patch_radius = us::any_cast(parsedArgs["patch_radius"]); int num_patches = 10; if (parsedArgs.count("num_patches")) num_patches = us::any_cast(parsedArgs["num_patches"]); try { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); mitk::Image::Pointer input_image = mitk::IOUtil::Load(imageName, &functor); typedef short DiffusionPixelType; typedef itk::VectorImage DwiImageType; typedef itk::Image DwiVolumeType; typedef itk::DiscreteGaussianImageFilter < DwiVolumeType, DwiVolumeType > GaussianFilterType; typedef itk::PatchBasedDenoisingImageFilter < DwiVolumeType, DwiVolumeType > NlmFilterType; typedef itk::VectorImageToImageFilter < DiffusionPixelType > ExtractFilterType; typedef itk::ComposeImageFilter < itk::Image > ComposeFilterType; if (!mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(input_image)) mitkThrow() << "Input is not a diffusion-weighted image!"; DwiImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(input_image); mitk::Image::Pointer denoised_image = nullptr; switch (type) { case 0: { ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i=0; iGetVectorLength(); ++i) { ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); extractor->SetIndex( i ); extractor->Update(); DwiVolumeType::Pointer gradient_volume = extractor->GetOutput(); itk::TotalVariationDenoisingImageFilter< DwiVolumeType, DwiVolumeType >::Pointer filter = itk::TotalVariationDenoisingImageFilter< DwiVolumeType, DwiVolumeType >::New(); filter->SetInput(gradient_volume); filter->SetLambda(lambda); filter->SetNumberIterations(tv_iterations); filter->Update(); composer->SetInput(i, filter->GetOutput()); } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); break; } case 1: { ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i = 0; i < itkVectorImagePointer->GetVectorLength(); ++i) { extractor->SetIndex(i); extractor->Update(); GaussianFilterType::Pointer filter = GaussianFilterType::New(); filter->SetVariance(variance); filter->SetInput(extractor->GetOutput()); filter->Update(); composer->SetInput(i, filter->GetOutput()); } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); break; } case 2: { typedef itk::Statistics::GaussianRandomSpatialNeighborSubsampler< NlmFilterType::PatchSampleType, DwiVolumeType::RegionType > SamplerType; // sampling the image to find similar patches SamplerType::Pointer sampler = SamplerType::New(); sampler->SetRadius( sampling_radius ); sampler->SetVariance( sampling_radius*sampling_radius ); sampler->SetNumberOfResultsRequested( num_patches ); MITK_INFO << "Starting NLM denoising"; ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i = 0; i < itkVectorImagePointer->GetVectorLength(); ++i) { extractor->SetIndex(i); extractor->Update(); NlmFilterType::Pointer filter = NlmFilterType::New(); filter->SetInput(extractor->GetOutput()); filter->SetPatchRadius(patch_radius); filter->SetNoiseModel(NlmFilterType::RICIAN); filter->UseSmoothDiscPatchWeightsOn(); filter->UseFastTensorComputationsOn(); filter->SetNumberOfIterations(nlm_iterations); filter->SetSmoothingWeight( 1 ); filter->SetKernelBandwidthEstimation(true); filter->SetSampler( sampler ); filter->Update(); composer->SetInput(i, filter->GetOutput()); MITK_INFO << "Gradient " << i << " finished"; } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); break; } } mitk::DiffusionPropertyHelper::SetGradientContainer(denoised_image, mitk::DiffusionPropertyHelper::GetGradientContainer(input_image)); mitk::DiffusionPropertyHelper::SetReferenceBValue(denoised_image, mitk::DiffusionPropertyHelper::GetReferenceBValue(input_image)); mitk::DiffusionPropertyHelper::InitializeImage( denoised_image ); std::string ext = itksys::SystemTools::GetFilenameExtension(outImage); if (ext==".nii" || ext==".nii.gz") mitk::IOUtil::Save(denoised_image, "DWI_NIFTI", outImage); else mitk::IOUtil::Save(denoised_image, outImage); } 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/Misc/FileFormatConverter.cpp b/Modules/DiffusionCmdApps/Misc/FileFormatConverter.cpp index fd2b1a0..55d2cee 100644 --- a/Modules/DiffusionCmdApps/Misc/FileFormatConverter.cpp +++ b/Modules/DiffusionCmdApps/Misc/FileFormatConverter.cpp @@ -1,80 +1,80 @@ /*=================================================================== 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 #include #include "mitkDiffusionCommandLineParser.h" #include #include /*! \brief Load image and save as specified file type. */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Format Converter"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Load image and save as specified file type."); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input file", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output file", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string inName = us::any_cast(parsedArgs["i"]); std::string outName = us::any_cast(parsedArgs["o"]); try { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); std::vector baseData = mitk::IOUtil::Load(inName, &functor); if ( baseData.size()>0 && dynamic_cast(baseData[0].GetPointer()) ) { mitk::IOUtil::Save(dynamic_cast(baseData[0].GetPointer()), outName.c_str()); } else if ( baseData.size()>0 && dynamic_cast(baseData[0].GetPointer()) ) { mitk::IOUtil::Save(dynamic_cast(baseData[0].GetPointer()) ,outName.c_str()); } else std::cout << "File type currently not supported!"; } 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/Misc/FlipPeaks.cpp b/Modules/DiffusionCmdApps/Misc/FlipPeaks.cpp index 6d60810..4936c96 100644 --- a/Modules/DiffusionCmdApps/Misc/FlipPeaks.cpp +++ b/Modules/DiffusionCmdApps/Misc/FlipPeaks.cpp @@ -1,104 +1,104 @@ /*=================================================================== 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 #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include /*! \brief Copies transformation matrix of one image to another */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Flip Peaks"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Flips the peaks of the input peak image along the given dimensions."); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input", "input image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output", "output image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("", "x", mitkDiffusionCommandLineParser::Bool, "Flip x", "flip along x-axis"); parser.addArgument("", "y", mitkDiffusionCommandLineParser::Bool, "Flip y", "flip along y-axis"); parser.addArgument("", "z", mitkDiffusionCommandLineParser::Bool, "Flip z", "flip along z-axis"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); bool x = false; if (parsedArgs.count("x")) x = us::any_cast(parsedArgs["x"]); bool y = false; if (parsedArgs.count("y")) y = us::any_cast(parsedArgs["y"]); bool z = false; if (parsedArgs.count("z")) z = us::any_cast(parsedArgs["z"]); try { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, std::vector()); mitk::PeakImage::Pointer image = mitk::IOUtil::Load(imageName, &functor); typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); mitk::PeakImage::ItkPeakImageType::Pointer itkImg = caster->GetOutput(); itk::FlipPeaksFilter< float >::Pointer flipper = itk::FlipPeaksFilter< float >::New(); flipper->SetInput(itkImg); flipper->SetFlipX(x); flipper->SetFlipY(y); flipper->SetFlipZ(z); flipper->Update(); mitk::Image::Pointer resultImage = dynamic_cast(mitk::PeakImage::New().GetPointer()); mitk::CastToMitkImage(flipper->GetOutput(), resultImage); resultImage->SetVolume(flipper->GetOutput()->GetBufferPointer()); mitk::IOUtil::Save(resultImage, outImage); } 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/Misc/HeadMotionCorrection.cpp b/Modules/DiffusionCmdApps/Misc/HeadMotionCorrection.cpp index 46c5d8d..970d171 100644 --- a/Modules/DiffusionCmdApps/Misc/HeadMotionCorrection.cpp +++ b/Modules/DiffusionCmdApps/Misc/HeadMotionCorrection.cpp @@ -1,78 +1,78 @@ /*=================================================================== 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 #include #include "mitkDiffusionCommandLineParser.h" #include #include #include int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("HeadMotionCorrection"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Simple affine head-motion correction tool"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image", 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; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); try { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); mitk::Image::Pointer in_image = mitk::IOUtil::Load(imageName, &functor); mitk::DWIHeadMotionCorrectionFilter::Pointer registerer = mitk::DWIHeadMotionCorrectionFilter::New(); registerer->SetInput(in_image); registerer->Update(); mitk::Image::Pointer out_image = registerer->GetCorrectedImage(); std::string ext = itksys::SystemTools::GetFilenameExtension(outImage); if (ext==".nii" || ext==".nii.gz") mitk::IOUtil::Save(out_image, "DWI_NIFTI", outImage); else mitk::IOUtil::Save(out_image, outImage); } 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/Misc/ResampleGradients.cpp b/Modules/DiffusionCmdApps/Misc/ResampleGradients.cpp index d7771d8..d465856 100644 --- a/Modules/DiffusionCmdApps/Misc/ResampleGradients.cpp +++ b/Modules/DiffusionCmdApps/Misc/ResampleGradients.cpp @@ -1,230 +1,230 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include "itkDWIVoxelFunctor.h" #include typedef short DiffusionPixelType; typedef itk::VectorImage< short, 3 > ItkDwiType; // itk includes #include "itkTimeProbe.h" #include "itkB0ImageExtractionImageFilter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" #include "itkBrainMaskExtractionImageFilter.h" #include "itkCastImageFilter.h" #include "itkVectorContainer.h" #include #include #include #include #include #include // Multishell includes #include // Multishell Functors #include #include #include #include // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "mitkProperties.h" #include "mitkVtkResliceInterpolationProperty.h" #include "mitkLookupTable.h" #include "mitkLookupTableProperty.h" #include "mitkTransferFunction.h" #include "mitkTransferFunctionProperty.h" //#include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include #include #include #include #include //#include //#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "mitkPreferenceListReaderOptionsFunctor.h" mitk::Image::Pointer DoReduceGradientDirections(mitk::Image::Pointer image, double BValue, unsigned int numOfGradientsToKeep, bool use_first_n) { bool isDiffusionImage( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(image) ); if ( !isDiffusionImage ) { std::cout << "Image is not a Diffusion Weighted Image" << std::endl; //return; } typedef itk::ElectrostaticRepulsionDiffusionGradientReductionFilter FilterType; typedef mitk::BValueMapProperty::BValueMap BValueMap; BValueMap shellSlectionMap; BValueMap originalShellMap = mitk::DiffusionPropertyHelper::GetBValueMap(image); std::vector newNumGradientDirections; //Keeps 1 b0 gradient double B0Value = 0; shellSlectionMap[B0Value] = originalShellMap[B0Value]; unsigned int num = 1; newNumGradientDirections.push_back(num); //BValue = 1000; shellSlectionMap[BValue] = originalShellMap[BValue]; //numOfGradientsToKeep = 32; newNumGradientDirections.push_back(numOfGradientsToKeep); if (newNumGradientDirections.empty()) { std::cout << "newNumGradientDirections is empty" << std::endl; //return; } auto gradientContainer = mitk::DiffusionPropertyHelper::GetGradientContainer(image); ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(image, itkVectorImagePointer); std::cout << "1" << std::endl; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkVectorImagePointer ); filter->SetOriginalGradientDirections(gradientContainer); filter->SetNumGradientDirections(newNumGradientDirections); filter->SetOriginalBValueMap(originalShellMap); filter->SetShellSelectionBValueMap(shellSlectionMap); filter->SetUseFirstN(use_first_n); filter->Update(); std::cout << "2" << std::endl; if( filter->GetOutput() == nullptr) { std::cout << "filter get output is nullptr" << std::endl; } mitk::Image::Pointer newImage = mitk::GrabItkImageMemory( filter->GetOutput() ); mitk::DiffusionPropertyHelper::CopyProperties(image, newImage, true); mitk::DiffusionPropertyHelper::SetGradientContainer(newImage, filter->GetGradientDirections()); mitk::DiffusionPropertyHelper::InitializeImage( newImage ); return newImage; } /*! \brief Resample gradients of input DWI image. */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Resample Gradients"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Resample gradients of input DWI image. You can select one b-value shell and the number of gradients within this shell you want to have. It will also keep one b0 image."); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image", us::Any(), false); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output image", us::Any(), false); parser.addArgument("b_value", "", mitkDiffusionCommandLineParser::Float, "b-value:", "float", 1000, false); parser.addArgument("num_gradients", "", mitkDiffusionCommandLineParser::Int, "Nr of gradients:", "integer", 32, false); parser.addArgument("use_first_n", "", mitkDiffusionCommandLineParser::Bool, "Use first N:", "no optimization, simply use first n gradients", 0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string outFileName = us::any_cast(parsedArgs["o"]); double bValue = us::any_cast(parsedArgs["b_value"]); unsigned int nrOfGradients = us::any_cast(parsedArgs["num_gradients"]); bool use_first_n = false; if (parsedArgs.count("use_first_n")) use_first_n = true; try { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({ "Diffusion Weighted Images" }, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({ "Diffusion Weighted Images" }, std::vector()); mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(inFileName, &functor); mitk::Image::Pointer newImage = DoReduceGradientDirections(mitkImage, bValue, nrOfGradients, use_first_n); //mitk::IOUtil::Save(newImage, outFileName); //save as dwi image mitk::IOUtil::Save(newImage, "DWI_NIFTI", outFileName); //save as nifti image } 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/Misc/RoundBvalues.cpp b/Modules/DiffusionCmdApps/Misc/RoundBvalues.cpp index 18ec92e..7685c64 100644 --- a/Modules/DiffusionCmdApps/Misc/RoundBvalues.cpp +++ b/Modules/DiffusionCmdApps/Misc/RoundBvalues.cpp @@ -1,106 +1,106 @@ /*=================================================================== 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 #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("RoundBvalues"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Round b-values"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("to_nearest", "", mitkDiffusionCommandLineParser::Int, "To nearest:", "integer", 1000); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); int to_nearest = 1000; if (parsedArgs.count("to_nearest")) to_nearest = us::any_cast(parsedArgs["to_nearest"]); try { typedef mitk::DiffusionPropertyHelper PropHelper; - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector({})); mitk::Image::Pointer in_image = mitk::IOUtil::Load(imageName, &functor); if (!PropHelper::IsDiffusionWeightedImage(in_image)) { mitkThrow() << "Input is not a diffusion weighted image: " << imageName; } typedef itk::DwiGradientLengthCorrectionFilter FilterType; auto itkVectorImagePointer = PropHelper::GetItkVectorImage(in_image); FilterType::Pointer filter = FilterType::New(); filter->SetRoundingValue(to_nearest); filter->SetReferenceBValue(PropHelper::GetReferenceBValue(in_image)); filter->SetReferenceGradientDirectionContainer(PropHelper::GetGradientContainer(in_image)); filter->Update(); mitk::Image::Pointer newImage = mitk::Image::New(); newImage->InitializeByItk( itkVectorImagePointer.GetPointer() ); newImage->SetImportVolume( itkVectorImagePointer->GetBufferPointer(), 0, 0, mitk::Image::CopyMemory); itkVectorImagePointer->GetPixelContainer()->ContainerManageMemoryOff(); PropHelper::CopyProperties(in_image, newImage, true); PropHelper::SetReferenceBValue(newImage, filter->GetNewBValue()); PropHelper::SetGradientContainer(newImage, filter->GetOutputGradientDirectionContainer()); PropHelper::InitializeImage(newImage); std::string ext = itksys::SystemTools::GetFilenameExtension(outImage); if (ext==".nii" || ext==".nii.gz") mitk::IOUtil::Save(newImage, "DWI_NIFTI", outImage); else mitk::IOUtil::Save(newImage, outImage); } 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/Misc/ShToOdfImage.cpp b/Modules/DiffusionCmdApps/Misc/ShToOdfImage.cpp index f3e6487..4d8d32c 100644 --- a/Modules/DiffusionCmdApps/Misc/ShToOdfImage.cpp +++ b/Modules/DiffusionCmdApps/Misc/ShToOdfImage.cpp @@ -1,72 +1,72 @@ /*=================================================================== 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 #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("ShToOdfImage"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Calculate discrete ODF image from SH coefficient image"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image", 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; // mandatory arguments std::string imageName = us::any_cast(parsedArgs["i"]); std::string outImage = us::any_cast(parsedArgs["o"]); try { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, std::vector()); mitk::ShImage::Pointer source = mitk::IOUtil::Load(imageName, &functor); mitk::Image::Pointer mitkImage = dynamic_cast(source.GetPointer()); mitk::OdfImage::Pointer out_image = mitk::convert::GetOdfFromShImage(mitkImage); if (out_image.IsNotNull()) mitk::IOUtil::Save(out_image, outImage); } 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/Python/BrainExtraction.cpp b/Modules/DiffusionCmdApps/Python/BrainExtraction.cpp index acd0d5d..bbb5934 100644 --- a/Modules/DiffusionCmdApps/Python/BrainExtraction.cpp +++ b/Modules/DiffusionCmdApps/Python/BrainExtraction.cpp @@ -1,179 +1,179 @@ /*=================================================================== 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 #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include #include typedef mitk::DiffusionPropertyHelper DPH; typedef itksys::SystemTools ist; std::string GetPythonFile(std::string filename, std::string exec_dir) { std::string out = ""; for (auto dir : mitk::bet::relative_search_dirs) { if ( ist::FileExists( exec_dir + dir + filename) ) { out = exec_dir + dir + filename; return out; } if ( ist::FileExists( ist::GetCurrentWorkingDirectory() + dir + filename) ) { out = ist::GetCurrentWorkingDirectory() + dir + filename; return out; } } for (auto dir : mitk::bet::absolute_search_dirs) { if ( ist::FileExists( dir + filename) ) { out = dir + filename; return out; } } return out; } int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("BrainExtraction"); parser.setCategory("Preprocessing Tools"); parser.setDescription("Performs brain extraction using a deep learning model"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output root", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string i = us::any_cast(parsedArgs["i"]); std::string o = us::any_cast(parsedArgs["o"]); std::string exec_dir = ist::GetFilenamePath(argv[0]); - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); mitk::Image::Pointer mitk_image = mitk::IOUtil::Load(i, &functor); bool missing_file = false; std::string missing_file_string = ""; if ( GetPythonFile("run_mitk.py", exec_dir).empty() ) { missing_file_string += "Brain extraction script file missing: run_mitk.py\n\n"; missing_file = true; } if ( GetPythonFile("model_final.model", exec_dir).empty() ) { missing_file_string += "Brain extraction model file missing: model_final.model\n\n"; missing_file = true; } if ( GetPythonFile("basic_config_just_like_braintumor.py", exec_dir).empty() ) { missing_file_string += "Config file missing: basic_config_just_like_braintumor.py\n\n"; missing_file = true; } if (missing_file) { mitkThrow() << missing_file_string; } 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->AddAbsoluteSearchDirs(mitk::bet::absolute_search_dirs); m_PythonService->AddRelativeSearchDirs(mitk::bet::relative_search_dirs); m_PythonService->Execute("paths=[]"); // set input files (model and config) m_PythonService->Execute("model_file=\""+GetPythonFile("model_final.model", exec_dir)+"\""); m_PythonService->Execute("config_file=\""+GetPythonFile("basic_config_just_like_braintumor.py", exec_dir)+"\""); // copy input image to python m_PythonService->CopyToPythonAsSimpleItkImage( mitk_image, "in_image"); // run segmentation script m_PythonService->ExecuteScript( GetPythonFile("run_mitk.py", exec_dir) ); // clean up after running script (better way than deleting individual variables?) if(m_PythonService->DoesVariableExist("in_image")) m_PythonService->Execute("del in_image"); // check for errors if(!m_PythonService->GetVariable("error_string").empty()) mitkThrow() << m_PythonService->GetVariable("error_string"); // get output images and add to datastorage std::string output_variables = m_PythonService->GetVariable("output_variables"); std::vector outputs; boost::split(outputs, output_variables, boost::is_any_of(",")); std::string output_types = m_PythonService->GetVariable("output_types"); std::vector types; boost::split(types, output_types, boost::is_any_of(",")); for (unsigned int i=0; iDoesVariableExist(outputs.at(i))) { mitk::Image::Pointer image = m_PythonService->CopySimpleItkImageFromPython(outputs.at(i)); if(types.at(i)=="input" && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(mitk_image)) { mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, image, true); mitk::DiffusionPropertyHelper::InitializeImage(image); } mitk::DataNode::Pointer corrected_node = mitk::DataNode::New(); corrected_node->SetData( image ); std::string name = o + "_"; name += outputs.at(i); if(types.at(i)=="input" && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(mitk_image)) mitk::IOUtil::Save(image, "DWI_NIFTI", name+".nii.gz"); else mitk::IOUtil::Save(image, name+".nii.gz"); } } MITK_INFO << "Finished brain extraction"; return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/Quantification/DiffusionIndices.cpp b/Modules/DiffusionCmdApps/Quantification/DiffusionIndices.cpp index 42f915c..4ab390a 100644 --- a/Modules/DiffusionCmdApps/Quantification/DiffusionIndices.cpp +++ b/Modules/DiffusionCmdApps/Quantification/DiffusionIndices.cpp @@ -1,196 +1,196 @@ /*=================================================================== 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 #include #include #include #include #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include #include /** * */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Diffusion Indices"); parser.setCategory("Diffusion Related Measures"); parser.setDescription("Computes requested diffusion related measures"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image (tensor, ODF or SH-coefficient image)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("index", "idx", mitkDiffusionCommandLineParser::String, "Index:", "index (fa, gfa, ra, ad, rd, ca, l2, l3, md, adc)", us::Any(), false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string index = us::any_cast(parsedArgs["index"]); std::string outFileName = us::any_cast(parsedArgs["o"]); std::string ext = itksys::SystemTools::GetFilenameLastExtension(outFileName); if (ext.empty()) outFileName += ".nii.gz"; try { // load input image - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "SH Image", "ODF Image", "Tensor Image"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "SH Image", "ODF Image", "Tensor Image"}, std::vector()); auto input = mitk::IOUtil::Load(inFileName, &functor); bool is_odf = (dynamic_cast(input.GetPointer()) || dynamic_cast(input.GetPointer())); bool is_dt = dynamic_cast(input.GetPointer()); bool is_dw = mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(input); if (is_odf) MITK_INFO << "Input is ODF image"; else if (is_dt) MITK_INFO << "Input is tensor image"; else if (is_dw) MITK_INFO << "Input is dMRI"; else { MITK_WARN << "Input is no ODF, SH, tensor or raw dMRI."; return EXIT_FAILURE; } mitk::LocaleSwitch localeSwitch("C"); if( is_odf && index=="gfa" ) { typedef itk::Vector OdfVectorType; typedef itk::Image OdfVectorImgType; OdfVectorImgType::Pointer itkvol; if (dynamic_cast(input.GetPointer())) { MITK_INFO << "Assuming MITK/MRtrix style SH convention!"; itkvol = mitk::convert::GetItkOdfFromShImage(input); } else itkvol = mitk::convert::GetItkOdfFromOdfImage(input); typedef itk::DiffusionOdfGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(itkvol); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(gfaFilter->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else if( is_dt ) { typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(input.GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); typedef itk::TensorDerivedMeasurementsFilter MeasurementsType; MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itk_dti.GetPointer() ); if(index=="fa") measurementsCalculator->SetMeasure(MeasurementsType::FA); else if(index=="ra") measurementsCalculator->SetMeasure(MeasurementsType::RA); else if(index=="ad") measurementsCalculator->SetMeasure(MeasurementsType::AD); else if(index=="rd") measurementsCalculator->SetMeasure(MeasurementsType::RD); else if(index=="ca") measurementsCalculator->SetMeasure(MeasurementsType::CA); else if(index=="l2") measurementsCalculator->SetMeasure(MeasurementsType::L2); else if(index=="l3") measurementsCalculator->SetMeasure(MeasurementsType::L3); else if(index=="md") measurementsCalculator->SetMeasure(MeasurementsType::MD); else { MITK_WARN << "No valid diffusion index for input image (tensor image) defined"; return EXIT_FAILURE; } measurementsCalculator->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(measurementsCalculator->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else if(is_dw && (index=="adc" || index=="md")) { typedef itk::AdcImageFilter< short, double > FilterType; auto itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(input); FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkVectorImagePointer ); filter->SetGradientDirections( mitk::DiffusionPropertyHelper::GetGradientContainer(input) ); filter->SetB_value( static_cast(mitk::DiffusionPropertyHelper::GetReferenceBValue(input)) ); if (index=="adc") filter->SetFitSignal(true); else filter->SetFitSignal(false); filter->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(filter->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else std::cout << "Diffusion index " << index << " not supported for supplied file type."; } 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/Quantification/DiffusionIvimFit.cpp b/Modules/DiffusionCmdApps/Quantification/DiffusionIvimFit.cpp index 256ce90..e5460a5 100644 --- a/Modules/DiffusionCmdApps/Quantification/DiffusionIvimFit.cpp +++ b/Modules/DiffusionCmdApps/Quantification/DiffusionIvimFit.cpp @@ -1,189 +1,189 @@ /*=================================================================== 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 "mitkDiffusionCommandLineParser.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 #include void IvimMapComputation( mitk::Image::Pointer input, std::string output_prefix , std::string output_type, double b_thresh, int type) { MITK_INFO << "Starting fit"; DPH::ImageType::Pointer vectorImage = DPH::ImageType::New(); mitk::CastToItkImage( input, vectorImage ); typedef itk::DiffusionIntravoxelIncoherentMotionReconstructionImageFilter IVIMFilterType; IVIMFilterType::Pointer ivim_filter = IVIMFilterType::New(); ivim_filter->SetInput( vectorImage ); ivim_filter->SetBValue( DPH::GetReferenceBValue( input.GetPointer() ) ); ivim_filter->SetGradientDirections( DPH::GetGradientContainer( input.GetPointer() ) ); switch (type) { case 0: ivim_filter->SetMethod(IVIMFilterType::IVIM_FIT_ALL); break; case 1: ivim_filter->SetMethod(IVIMFilterType::IVIM_DSTAR_FIX); break; case 2: ivim_filter->SetMethod(IVIMFilterType::IVIM_D_THEN_DSTAR); break; case 3: ivim_filter->SetMethod(IVIMFilterType::IVIM_LINEAR_D_THEN_F); break; default: ivim_filter->SetMethod(IVIMFilterType::IVIM_D_THEN_DSTAR); } ivim_filter->SetBThres(b_thresh); ivim_filter->SetS0Thres(0); ivim_filter->SetFitDStar(true); ivim_filter->SetNumberOfThreads(1); try { ivim_filter->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "IVIM fit failed with an ITK Exception: " << e.what(); } mitk::Image::Pointer f_image = mitk::Image::New(); f_image->InitializeByItk( ivim_filter->GetOutput() ); f_image->SetVolume( ivim_filter->GetOutput()->GetBufferPointer() ); mitk::Image::Pointer d_image = mitk::Image::New(); d_image->InitializeByItk( ivim_filter->GetOutput(1) ); d_image->SetVolume( ivim_filter->GetOutput(1)->GetBufferPointer() ); mitk::Image::Pointer dstar_image = mitk::Image::New(); dstar_image->InitializeByItk( ivim_filter->GetOutput(1) ); dstar_image->SetVolume( ivim_filter->GetOutput(2)->GetBufferPointer() ); std::string outputf_FileName = output_prefix + "_f_map." + output_type; std::string outputD_FileName = output_prefix + "_D_map." + output_type; std::string outputDstar_FileName = output_prefix + "_Dstar_map." + output_type; try { mitk::IOUtil::Save( dstar_image, outputDstar_FileName ); mitk::IOUtil::Save( d_image, outputD_FileName ); mitk::IOUtil::Save( f_image, outputf_FileName ); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Failed to save the KurtosisFit Results due to exception: " << e.what(); } } int main( int argc, char* argv[] ) { mitkDiffusionCommandLineParser parser; parser.setTitle("Diffusion IVIM Fit"); parser.setCategory("Diffusion Related Measures"); parser.setContributor("MIC"); parser.setDescription("Fitting IVIM"); parser.setArgumentPrefix("--","-"); // mandatory arguments parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input: ", "input image (DWI)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output Preifx: ", "Prefix for the output images, will append _ADC, _AKC accordingly ", us::Any(), false); parser.addArgument("output_type", "", mitkDiffusionCommandLineParser::String, "Output Type: ", "choose data type of output image, e.g. '.nii' or '.nrrd' ", std::string(".nrrd")); parser.addArgument("b_threshold", "", mitkDiffusionCommandLineParser::Float, "b-threshold:", "Omit smaller b-values for first fit^", 170.0); parser.addArgument("fit_type", "", mitkDiffusionCommandLineParser::Int, "Fit:", "Jointly fit D, f and D* (0); Fit D&f with fixed D* (1); Fit D&f (high b), then fit D* (2); Linearly fit D&f (high b), then fit D* (3)", 2); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string inFileName = us::any_cast(parsedArgs["i"]); std::string out_prefix = us::any_cast(parsedArgs["o"]); - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); mitk::Image::Pointer inputImage = mitk::IOUtil::Load(inFileName, &functor); double b_thresh = 170; int fit_type = 2; std::string out_type = "nrrd"; if (parsedArgs.count("output_type")) out_type = us::any_cast(parsedArgs["output_type"]); if (parsedArgs.count("b_threshold")) b_thresh = us::any_cast(parsedArgs["b_threshold"]); if (parsedArgs.count("fit_type")) fit_type = us::any_cast(parsedArgs["fit_type"]); if( !DPH::IsDiffusionWeightedImage( inputImage ) ) { MITK_ERROR("DiffusionIVIMFit.Input") << "No valid diffusion-weighted image provided, failed to load " << inFileName << " as DW Image. Aborting..."; return EXIT_FAILURE; } IvimMapComputation( inputImage, out_prefix , out_type, b_thresh, fit_type); } diff --git a/Modules/DiffusionCmdApps/Quantification/DiffusionKurtosisFit.cpp b/Modules/DiffusionCmdApps/Quantification/DiffusionKurtosisFit.cpp index 72b5cc4..f34637b 100644 --- a/Modules/DiffusionCmdApps/Quantification/DiffusionKurtosisFit.cpp +++ b/Modules/DiffusionCmdApps/Quantification/DiffusionKurtosisFit.cpp @@ -1,254 +1,254 @@ /*=================================================================== 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 "mitkDiffusionCommandLineParser.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 #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 , 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->SetOmitUnweightedValue(omitBZero); kurtosis_filter->SetBoundariesForKurtosis(-lower,upper); // kurtosis_filter->SetInitialSolution(const vnl_vector& x0 ); if(maskPath != "") { mitk::Image::Pointer segmentation; segmentation = mitk::IOUtil::Load(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." + 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 ); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Failed to save the KurtosisFit Results due to exception: " << e.what(); } } int main( int argc, char* argv[] ) { mitkDiffusionCommandLineParser parser; parser.setTitle("Diffusion Kurtosis Fit"); parser.setCategory("Diffusion Related Measures"); parser.setContributor("MIC"); parser.setDescription("Fitting Kurtosis"); parser.setArgumentPrefix("--","-"); // mandatory arguments parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input: ", "input image (DWI)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output Preifx: ", "Prefix for the output images, will append _ADC, _AKC accordingly ", us::Any(), false); parser.addArgument("output_type", "", mitkDiffusionCommandLineParser::String, "Output Type: ", "choose data type of output image, e.g. '.nii' or '.nrrd' "); // optional arguments parser.addArgument("mask", "m", mitkDiffusionCommandLineParser::String, "Masking Image: ", "ROI (segmentation)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("help", "h", mitkDiffusionCommandLineParser::Bool, "Help", "Show this help text"); parser.addArgument("omitbzero", "om", mitkDiffusionCommandLineParser::Bool, "Omit b0:", "Omit b0 value during fit (default = false)", us::Any()); parser.addArgument("lowerkbound", "kl", mitkDiffusionCommandLineParser::Float, "lower Kbound:", "Set (unsigned) lower boundary for Kurtosis parameter (default = -1000)", us::Any()); parser.addArgument("upperkbound", "ku", mitkDiffusionCommandLineParser::Float, "upper Kbound:", "Set upper boundary for Kurtosis parameter (default = 1000)", us::Any()); std::map parsedArgs = parser.parseArguments(argc, argv); 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["i"]); std::string out_prefix = us::any_cast(parsedArgs["o"]); std::string maskPath = ""; - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); mitk::Image::Pointer inputImage = mitk::IOUtil::Load(inFileName, &functor); bool omitBZero = false; double lower = -1000; double upper = 1000; std::string out_type = "nrrd"; if (parsedArgs.count("mask") || parsedArgs.count("m")) { maskPath = us::any_cast(parsedArgs["mask"]); } if (parsedArgs.count("output_type") || parsedArgs.count("ot")) { out_type = us::any_cast(parsedArgs["output_type"]); } 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"]); } if (parsedArgs.count("upperkbound") || parsedArgs.count("ku")) { upper = us::any_cast(parsedArgs["upperkbound"]); } if( !DPH::IsDiffusionWeightedImage( inputImage ) ) { MITK_ERROR("DiffusionKurtosisFit.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); } diff --git a/Modules/DiffusionCmdApps/Quantification/MultishellMethods.cpp b/Modules/DiffusionCmdApps/Quantification/MultishellMethods.cpp index 118890c..a8b2c33 100644 --- a/Modules/DiffusionCmdApps/Quantification/MultishellMethods.cpp +++ b/Modules/DiffusionCmdApps/Quantification/MultishellMethods.cpp @@ -1,215 +1,215 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Multishell Methods"); parser.setCategory("Preprocessing Tools"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input file", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output file", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("adc", "D", mitkDiffusionCommandLineParser::Bool, "ADC:", "ADC Average", us::Any(), false); parser.addArgument("akc", "K", mitkDiffusionCommandLineParser::Bool, "Kurtosis fit:", "Kurtosis Fit", us::Any(), false); parser.addArgument("biexp", "B", mitkDiffusionCommandLineParser::Bool, "BiExp fit:", "BiExp fit", us::Any(), false); parser.addArgument("targetbvalue", "b", mitkDiffusionCommandLineParser::String, "b Value:", "target bValue (mean, min, max)", us::Any(), false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments std::string inName = us::any_cast(parsedArgs["i"]); std::string outName = us::any_cast(parsedArgs["o"]); bool applyADC = us::any_cast(parsedArgs["adc"]); bool applyAKC = us::any_cast(parsedArgs["akc"]); bool applyBiExp = us::any_cast(parsedArgs["biexp"]); std::string targetType = us::any_cast(parsedArgs["targetbvalue"]); try { std::cout << "Loading " << inName; - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); mitk::Image::Pointer dwi = mitk::IOUtil::Load(inName, &functor); if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dwi ) ) { typedef itk::RadialMultishellToSingleshellImageFilter FilterType; typedef itk::DwiGradientLengthCorrectionFilter CorrectionFilterType; CorrectionFilterType::Pointer roundfilter = CorrectionFilterType::New(); roundfilter->SetRoundingValue( 1000 ); roundfilter->SetReferenceBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi )); roundfilter->SetReferenceGradientDirectionContainer(mitk::DiffusionPropertyHelper::GetGradientContainer(dwi)); roundfilter->Update(); mitk::DiffusionPropertyHelper::SetReferenceBValue(dwi, roundfilter->GetNewBValue()); mitk::DiffusionPropertyHelper::SetGradientContainer(dwi, roundfilter->GetOutputGradientDirectionContainer()); // filter input parameter const mitk::DiffusionPropertyHelper::BValueMapType &originalShellMap = mitk::DiffusionPropertyHelper::GetBValueMap(dwi); mitk::DiffusionPropertyHelper::ImageType::Pointer vectorImage = mitk::DiffusionPropertyHelper::ImageType::New(); mitk::CastToItkImage(dwi, vectorImage); const mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientContainer = mitk::DiffusionPropertyHelper::GetGradientContainer(dwi); const unsigned int &bValue = mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi ); // filter call vnl_vector bValueList(originalShellMap.size()-1); double targetBValue = bValueList.mean(); mitk::DiffusionPropertyHelper::BValueMapType::const_iterator it = originalShellMap.begin(); ++it; int i = 0 ; for(; it != originalShellMap.end(); ++it) bValueList.put(i++,it->first); if( targetType == "mean" ) targetBValue = bValueList.mean(); else if( targetType == "min" ) targetBValue = bValueList.min_value(); else if( targetType == "max" ) targetBValue = bValueList.max_value(); if(applyADC) { FilterType::Pointer filter = FilterType::New(); filter->SetInput(vectorImage); filter->SetOriginalGradientDirections(gradientContainer); filter->SetOriginalBValueMap(originalShellMap); filter->SetOriginalBValue(bValue); itk::ADCAverageFunctor::Pointer functor = itk::ADCAverageFunctor::New(); functor->setListOfBValues(bValueList); functor->setTargetBValue(targetBValue); filter->SetFunctor(functor); filter->Update(); // create new DWI image mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() ); mitk::DiffusionPropertyHelper::SetReferenceBValue(outImage, targetBValue); mitk::DiffusionPropertyHelper::SetGradientContainer(outImage, filter->GetTargetGradientDirections()); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, (outName + "_ADC.dwi").c_str()); } if(applyAKC) { FilterType::Pointer filter = FilterType::New(); filter->SetInput(vectorImage); filter->SetOriginalGradientDirections(gradientContainer); filter->SetOriginalBValueMap(originalShellMap); filter->SetOriginalBValue(bValue); itk::KurtosisFitFunctor::Pointer functor = itk::KurtosisFitFunctor::New(); functor->setListOfBValues(bValueList); functor->setTargetBValue(targetBValue); filter->SetFunctor(functor); filter->Update(); // create new DWI image mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() ); mitk::DiffusionPropertyHelper::SetReferenceBValue(outImage, targetBValue); mitk::DiffusionPropertyHelper::SetGradientContainer(outImage, filter->GetTargetGradientDirections()); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, (std::string(outName) + "_AKC.dwi").c_str()); } if(applyBiExp) { FilterType::Pointer filter = FilterType::New(); filter->SetInput(vectorImage); filter->SetOriginalGradientDirections(gradientContainer); filter->SetOriginalBValueMap(originalShellMap); filter->SetOriginalBValue(bValue); itk::BiExpFitFunctor::Pointer functor = itk::BiExpFitFunctor::New(); functor->setListOfBValues(bValueList); functor->setTargetBValue(targetBValue); filter->SetFunctor(functor); filter->Update(); // create new DWI image mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() ); mitk::DiffusionPropertyHelper::SetReferenceBValue(outImage, targetBValue); mitk::DiffusionPropertyHelper::SetGradientContainer(outImage, filter->GetTargetGradientDirections()); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, (std::string(outName) + "_BiExp.dwi").c_str()); } } } 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/Quantification/QballReconstruction.cpp b/Modules/DiffusionCmdApps/Quantification/QballReconstruction.cpp index 18a1b9a..d8e2240 100644 --- a/Modules/DiffusionCmdApps/Quantification/QballReconstruction.cpp +++ b/Modules/DiffusionCmdApps/Quantification/QballReconstruction.cpp @@ -1,266 +1,266 @@ /*=================================================================== 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 "mitkImage.h" #include "itkAnalyticalDiffusionQballReconstructionImageFilter.h" #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include template void TemplatedMultishellQBallReconstruction(float lambda, mitk::Image::Pointer dwi, bool output_sampled, int threshold, std::string outfilename) { typedef itk::DiffusionMultiShellQballReconstructionImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); auto bMap = mitk::DiffusionPropertyHelper::GetBValueMap(dwi); auto it1 = bMap.rbegin(); auto it2 = bMap.rbegin(); ++it2; // Get average distance int avdistance = 0; for(; it2 != bMap.rend(); ++it1, ++it2) avdistance += static_cast(it1->first - it2->first); avdistance /= bMap.size()-1; // Check if all shells are using the same averae distance it1 = bMap.rbegin(); it2 = bMap.rbegin(); ++it2; for(; it2 != bMap.rend(); ++it1,++it2) { if(avdistance != static_cast(it1->first - it2->first)) { mitkThrow() << "Shells are not equidistant."; } } auto itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(dwi); filter->SetBValueMap(bMap); filter->SetGradientImage(mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer, mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi)); filter->SetThreshold(static_cast(threshold)); filter->SetLambda(static_cast(lambda)); filter->Update(); mitk::OdfImage::Pointer image = mitk::OdfImage::New(); mitk::Image::Pointer coeffsImage = dynamic_cast(mitk::ShImage::New().GetPointer()); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); std::string coeffout = outfilename; coeffout += ".nii.gz"; mitk::IOUtil::Save(coeffsImage, "SH_IMAGE", coeffout); outfilename += ".odf"; if (output_sampled) mitk::IOUtil::Save(image, outfilename); } template void TemplatedCsaQBallReconstruction(float lambda, mitk::Image::Pointer dwi, bool output_sampled, int threshold, std::string outfilename) { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; auto itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(dwi); FilterType::Pointer filter = FilterType::New(); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi)); filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetGradientContainer(dwi), itkVectorImagePointer ); filter->SetThreshold(static_cast(threshold)); filter->SetLambda(static_cast(lambda)); // filter->SetUseMrtrixBasis(mrTrix); filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); filter->Update(); mitk::OdfImage::Pointer image = mitk::OdfImage::New(); mitk::Image::Pointer coeffsImage = dynamic_cast(mitk::ShImage::New().GetPointer()); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); std::string coeffout = outfilename; coeffout += ".nii.gz"; mitk::IOUtil::Save(coeffsImage, "SH_IMAGE", coeffout); outfilename += ".odf"; if (output_sampled) mitk::IOUtil::Save(image, outfilename); } /** * Perform Q-ball reconstruction using a spherical harmonics basis */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input image", "input raw dwi (.dwi or .nii/.nii.gz)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output image", "output image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("sh_order", "", mitkDiffusionCommandLineParser::Int, "Spherical harmonics order", "spherical harmonics order", 4); parser.addArgument("b0_threshold", "", mitkDiffusionCommandLineParser::Int, "b0 threshold", "baseline image intensity threshold", 0); parser.addArgument("round_bvalues", "", mitkDiffusionCommandLineParser::Int, "Round b-values", "round to specified integer", 0); parser.addArgument("lambda", "", mitkDiffusionCommandLineParser::Float, "Lambda", "ragularization factor lambda", 0.006); parser.addArgument("output_sampled", "", mitkDiffusionCommandLineParser::Bool, "Output sampled ODFs", "output file containing the sampled ODFs"); parser.setCategory("Signal Modelling"); parser.setTitle("Qball Reconstruction"); parser.setDescription(""); parser.setContributor("MIC"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string outfilename = us::any_cast(parsedArgs["o"]); if (itksys::SystemTools::GetFilenamePath(outfilename).size()>0) outfilename = itksys::SystemTools::GetFilenamePath(outfilename)+"/"+itksys::SystemTools::GetFilenameWithoutExtension(outfilename); else outfilename = itksys::SystemTools::GetFilenameWithoutExtension(outfilename); int threshold = 0; if (parsedArgs.count("b0_threshold")) threshold = us::any_cast(parsedArgs["b0_threshold"]); int round_bvalues = 0; if (parsedArgs.count("round_bvalues")) round_bvalues = us::any_cast(parsedArgs["round_bvalues"]); int shOrder = 4; if (parsedArgs.count("sh_order")) shOrder = us::any_cast(parsedArgs["sh_order"]); float lambda = 0.006f; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); bool outCoeffs = false; if (parsedArgs.count("output_coeffs")) outCoeffs = us::any_cast(parsedArgs["output_coeffs"]); // bool mrTrix = false; // if (parsedArgs.count("mrtrix")) // mrTrix = us::any_cast(parsedArgs["mrtrix"]); try { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); std::vector< mitk::BaseData::Pointer > infile = mitk::IOUtil::Load(inFileName, &functor); mitk::Image::Pointer dwi = dynamic_cast(infile.at(0).GetPointer()); if (round_bvalues>0) { MITK_INFO << "Rounding b-values"; typedef itk::DwiGradientLengthCorrectionFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetRoundingValue(round_bvalues); filter->SetReferenceBValue(static_cast(mitk::DiffusionPropertyHelper::GetReferenceBValue(dwi))); filter->SetReferenceGradientDirectionContainer(mitk::DiffusionPropertyHelper::GetGradientContainer(dwi)); filter->Update(); mitk::DiffusionPropertyHelper::SetReferenceBValue(dwi, static_cast(filter->GetNewBValue())); mitk::DiffusionPropertyHelper::CopyProperties(dwi, dwi, true); mitk::DiffusionPropertyHelper::SetGradientContainer(dwi, filter->GetOutputGradientDirectionContainer()); mitk::DiffusionPropertyHelper::InitializeImage(dwi); } auto bMap = mitk::DiffusionPropertyHelper::GetBValueMap(dwi); if(bMap.size()!=4 && bMap.size()!=2) mitkThrow() << "Only three equidistant shells or a single shell are supported. Found " << bMap.size(); MITK_INFO << "Averaging redundant gradients"; mitk::DiffusionPropertyHelper::AverageRedundantGradients(dwi, 0.001); MITK_INFO << "SH order: " << shOrder; MITK_INFO << "lambda: " << lambda; MITK_INFO << "B0 threshold: " << threshold; MITK_INFO << "Round bvalues: " << round_bvalues; switch ( shOrder ) { case 4: { if(bMap.size()==2) TemplatedCsaQBallReconstruction<4>(lambda, dwi, outCoeffs, threshold, outfilename); else if(bMap.size()==4) TemplatedMultishellQBallReconstruction<4>(lambda, dwi, outCoeffs, threshold, outfilename); break; } case 6: { if(bMap.size()==2) TemplatedCsaQBallReconstruction<6>(lambda, dwi, outCoeffs, threshold, outfilename); else if(bMap.size()==4) TemplatedMultishellQBallReconstruction<6>(lambda, dwi, outCoeffs, threshold, outfilename); break; } case 8: { if(bMap.size()==2) TemplatedCsaQBallReconstruction<8>(lambda, dwi, outCoeffs, threshold, outfilename); else if(bMap.size()==4) TemplatedMultishellQBallReconstruction<8>(lambda, dwi, outCoeffs, threshold, outfilename); break; } case 10: { if(bMap.size()==2) TemplatedCsaQBallReconstruction<10>(lambda, dwi, outCoeffs, threshold, outfilename); else if(bMap.size()==4) TemplatedMultishellQBallReconstruction<10>(lambda, dwi, outCoeffs, threshold, outfilename); break; } case 12: { if(bMap.size()==2) TemplatedCsaQBallReconstruction<12>(lambda, dwi, outCoeffs, threshold, outfilename); else if(bMap.size()==4) TemplatedMultishellQBallReconstruction<12>(lambda, dwi, outCoeffs, threshold, outfilename); break; } default: { mitkThrow() << "SH order not supported"; } } } catch ( itk::ExceptionObject &err) { std::cout << "Exception: " << err; } catch ( std::exception& err) { std::cout << "Exception: " << err.what(); } catch ( ... ) { std::cout << "Exception!"; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp b/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp index fd20d77..40cd281 100644 --- a/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp +++ b/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp @@ -1,134 +1,134 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include /*! \brief Perform global fiber tractography (Gibbs tractography) */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Gibbs Tracking"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform global fiber tractography (Gibbs tractography)"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image (tensor, ODF or SH-coefficient image)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output tractogram", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("parameters", "", mitkDiffusionCommandLineParser::String, "Parameters:", "parameter file (.gtp)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("mask", "", mitkDiffusionCommandLineParser::String, "Mask:", "binary mask image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string paramFileName = us::any_cast(parsedArgs["parameters"]); std::string outFileName = us::any_cast(parsedArgs["o"]); try { // instantiate gibbs tracker typedef itk::Vector OdfVectorType; typedef itk::Image ItkOdfImageType; typedef itk::GibbsTrackingFilter GibbsTrackingFilterType; GibbsTrackingFilterType::Pointer gibbsTracker = GibbsTrackingFilterType::New(); // load input image - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, std::vector()); mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(inFileName, &functor); // try to cast to Odf image if( dynamic_cast(mitkImage.GetPointer()) ) { mitk::OdfImage::Pointer mitkOdfImage = dynamic_cast(mitkImage.GetPointer()); ItkOdfImageType::Pointer itk_odf = ItkOdfImageType::New(); mitk::CastToItkImage(mitkOdfImage, itk_odf); gibbsTracker->SetOdfImage(itk_odf.GetPointer()); } else if( dynamic_cast(mitkImage.GetPointer()) ) { typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(mitkImage.GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); gibbsTracker->SetTensorImage(itk_dti); } else if ( dynamic_cast(mitkImage.GetPointer()) ) { MITK_INFO << "Assuming MITK/MRtrix style SH convention!"; mitk::Image::Pointer shImage = dynamic_cast(mitkImage.GetPointer()); gibbsTracker->SetOdfImage(mitk::convert::GetItkOdfFromShImage(shImage)); } else return EXIT_FAILURE; // global tracking if (parsedArgs.count("mask")) { typedef itk::Image MaskImgType; mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::Load(us::any_cast(parsedArgs["mask"])); MaskImgType::Pointer itk_mask = MaskImgType::New(); mitk::CastToItkImage(mitkMaskImage, itk_mask); gibbsTracker->SetMaskImage(itk_mask); } gibbsTracker->SetDuplicateImage(false); gibbsTracker->SetLoadParameterFile( paramFileName ); // gibbsTracker->SetLutPath( "" ); gibbsTracker->Update(); mitk::FiberBundle::Pointer mitkFiberBundle = mitk::FiberBundle::New(gibbsTracker->GetFiberBundle()); mitkFiberBundle->SetTrackVisHeader(mitkImage->GetGeometry()); mitk::IOUtil::Save(mitkFiberBundle, outFileName ); } 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/Tractography/RfTraining.cpp b/Modules/DiffusionCmdApps/Tractography/RfTraining.cpp index 67d200e..bf16143 100644 --- a/Modules/DiffusionCmdApps/Tractography/RfTraining.cpp +++ b/Modules/DiffusionCmdApps/Tractography/RfTraining.cpp @@ -1,238 +1,239 @@ /*=================================================================== 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 #include #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include /*! \brief Train random forest classifier for machine learning based streamline tractography */ int main(int argc, char* argv[]) { MITK_INFO << "RfTraining"; mitkDiffusionCommandLineParser parser; parser.setTitle("Trains Random Forests for Machine Learning Based Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Train random forest classifier for machine learning based streamline tractography"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::StringList, "DWIs:", "input diffusion-weighted images", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "t", mitkDiffusionCommandLineParser::StringList, "Tractograms:", "input training tractograms", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Forest:", "output random forest (HDF5)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.endGroup(); parser.beginGroup("2. Additional input images:"); parser.addArgument("masks", "", mitkDiffusionCommandLineParser::StringList, "Masks:", "restrict training using a binary mask image", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("wm_masks", "", mitkDiffusionCommandLineParser::StringList, "WM-Masks:", "if no binary white matter mask is specified, the envelope of the input tractogram is used", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("volume_modification_images", "", mitkDiffusionCommandLineParser::StringList, "Volume modification images:", "specify a list of float images that modify the fiber density", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("additional_feature_images", "", mitkDiffusionCommandLineParser::StringList, "Additional feature images:", "specify a list of float images that hold additional features (float)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.endGroup(); parser.beginGroup("3. Forest parameters:"); parser.addArgument("num_trees", "", mitkDiffusionCommandLineParser::Int, "Number of trees:", "number of trees", 30); parser.addArgument("max_tree_depth", "", mitkDiffusionCommandLineParser::Int, "Max. tree depth:", "maximum tree depth", 25); parser.addArgument("sample_fraction", "", mitkDiffusionCommandLineParser::Float, "Sample fraction:", "fraction of samples used per tree", 0.7); parser.endGroup(); parser.beginGroup("4. Feature parameters:"); parser.addArgument("use_sh_features", "", mitkDiffusionCommandLineParser::Bool, "Use SH features:", "use SH features", false); parser.addArgument("sampling_distance", "", mitkDiffusionCommandLineParser::Float, "Sampling distance:", "resampling parameter for the input tractogram in mm (determines number of white-matter samples)", us::Any()); parser.addArgument("max_wm_samples", "", mitkDiffusionCommandLineParser::Int, "Max. num. WM samples:", "upper limit for the number of WM samples"); parser.addArgument("num_gm_samples", "", mitkDiffusionCommandLineParser::Int, "Number of gray matter samples per voxel:", "Number of gray matter samples per voxel", us::Any()); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; bool shfeatures = false; if (parsedArgs.count("use_sh_features")) shfeatures = us::any_cast(parsedArgs["use_sh_features"]); mitkDiffusionCommandLineParser::StringContainerType imageFiles = us::any_cast(parsedArgs["i"]); mitkDiffusionCommandLineParser::StringContainerType wmMaskFiles; if (parsedArgs.count("wm_masks")) wmMaskFiles = us::any_cast(parsedArgs["wm_masks"]); mitkDiffusionCommandLineParser::StringContainerType volModFiles; if (parsedArgs.count("volume_modification_images")) volModFiles = us::any_cast(parsedArgs["volume_modification_images"]); mitkDiffusionCommandLineParser::StringContainerType addFeatFiles; if (parsedArgs.count("additional_feature_images")) addFeatFiles = us::any_cast(parsedArgs["additional_feature_images"]); mitkDiffusionCommandLineParser::StringContainerType maskFiles; if (parsedArgs.count("masks")) maskFiles = us::any_cast(parsedArgs["masks"]); std::string forestFile = us::any_cast(parsedArgs["o"]); mitkDiffusionCommandLineParser::StringContainerType tractogramFiles; if (parsedArgs.count("t")) tractogramFiles = us::any_cast(parsedArgs["t"]); int num_trees = 30; if (parsedArgs.count("num_trees")) num_trees = us::any_cast(parsedArgs["num_trees"]); int gm_samples = -1; if (parsedArgs.count("num_gm_samples")) gm_samples = us::any_cast(parsedArgs["num_gm_samples"]); float sampling_distance = -1; if (parsedArgs.count("sampling_distance")) sampling_distance = us::any_cast(parsedArgs["sampling_distance"]); int max_tree_depth = 25; if (parsedArgs.count("max_tree_depth")) max_tree_depth = us::any_cast(parsedArgs["max_tree_depth"]); double sample_fraction = 0.7; if (parsedArgs.count("sample_fraction")) sample_fraction = us::any_cast(parsedArgs["sample_fraction"]); int maxWmSamples = -1; if (parsedArgs.count("max_wm_samples")) maxWmSamples = us::any_cast(parsedArgs["max_wm_samples"]); MITK_INFO << "loading diffusion-weighted images"; std::vector< mitk::Image::Pointer > rawData; - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, std::vector()); for (auto imgFile : imageFiles) { auto dwi = mitk::IOUtil::Load(imgFile, &functor); rawData.push_back(dwi); } typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUcharImgType; MITK_INFO << "loading mask images"; std::vector< ItkUcharImgType::Pointer > maskImageVector; for (auto maskFile : maskFiles) { mitk::Image::Pointer img = mitk::IOUtil::Load(maskFile); ItkUcharImgType::Pointer mask = ItkUcharImgType::New(); mitk::CastToItkImage(img, mask); maskImageVector.push_back(mask); } MITK_INFO << "loading white matter mask images"; std::vector< ItkUcharImgType::Pointer > wmMaskImageVector; for (auto wmFile : wmMaskFiles) { mitk::Image::Pointer img = mitk::IOUtil::Load(wmFile); ItkUcharImgType::Pointer wmmask = ItkUcharImgType::New(); mitk::CastToItkImage(img, wmmask); wmMaskImageVector.push_back(wmmask); } MITK_INFO << "loading tractograms"; std::vector< mitk::FiberBundle::Pointer > tractograms; for (auto tractFile : tractogramFiles) { mitk::FiberBundle::Pointer fib = mitk::IOUtil::Load(tractFile); tractograms.push_back(fib); } MITK_INFO << "loading white volume modification images"; std::vector< ItkFloatImgType::Pointer > volumeModImages; for (auto file : volModFiles) { mitk::Image::Pointer img = mitk::IOUtil::Load(file); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); volumeModImages.push_back(itkimg); } MITK_INFO << "loading additional feature images"; std::vector< std::vector< ItkFloatImgType::Pointer > > addFeatImages; for (std::size_t i=0; i()); int c = 0; for (auto file : addFeatFiles) { mitk::Image::Pointer img = mitk::IOUtil::Load(file); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); addFeatImages.at(c%addFeatImages.size()).push_back(itkimg); c++; } mitk::TractographyForest::Pointer forest = nullptr; if (shfeatures) { mitk::TrackingHandlerRandomForest<6,28> forestHandler; forestHandler.SetDwis(rawData); forestHandler.SetMaskImages(maskImageVector); forestHandler.SetWhiteMatterImages(wmMaskImageVector); forestHandler.SetFiberVolumeModImages(volumeModImages); forestHandler.SetAdditionalFeatureImages(addFeatImages); forestHandler.SetTractograms(tractograms); forestHandler.SetNumTrees(num_trees); forestHandler.SetMaxTreeDepth(max_tree_depth); forestHandler.SetGrayMatterSamplesPerVoxel(gm_samples); forestHandler.SetSampleFraction(sample_fraction); forestHandler.SetFiberSamplingStep(sampling_distance); forestHandler.SetMaxNumWmSamples(maxWmSamples); forestHandler.StartTraining(); forest = forestHandler.GetForest(); } else { mitk::TrackingHandlerRandomForest<6,100> forestHandler; forestHandler.SetDwis(rawData); forestHandler.SetMaskImages(maskImageVector); forestHandler.SetWhiteMatterImages(wmMaskImageVector); forestHandler.SetFiberVolumeModImages(volumeModImages); forestHandler.SetAdditionalFeatureImages(addFeatImages); forestHandler.SetTractograms(tractograms); forestHandler.SetNumTrees(num_trees); forestHandler.SetMaxTreeDepth(max_tree_depth); forestHandler.SetGrayMatterSamplesPerVoxel(gm_samples); forestHandler.SetSampleFraction(sample_fraction); forestHandler.SetFiberSamplingStep(sampling_distance); forestHandler.SetMaxNumWmSamples(maxWmSamples); forestHandler.StartTraining(); forest = forestHandler.GetForest(); } mitk::IOUtil::Save(forest, forestFile); return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp b/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp index e03c52b..1d6b327 100644 --- a/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp @@ -1,577 +1,580 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include -#include #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include const int numOdfSamples = 200; typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType; /*! \brief */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform streamline tractography"); parser.setContributor("MIC"); // parameters fo all methods parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::StringList, "Input:", "input image (multiple possible for 'DetTensor' algorithm)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output fiberbundle/probability map", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("type", "", mitkDiffusionCommandLineParser::String, "Type:", "which tracker to use (Peaks; Tensor; ODF; ODF-DIPY/FSL; RF)", us::Any(), false); parser.addArgument("probabilistic", "", mitkDiffusionCommandLineParser::Bool, "Probabilistic:", "Probabilistic tractography", us::Any(false)); parser.endGroup(); parser.beginGroup("2. Seeding:"); parser.addArgument("seeds", "", mitkDiffusionCommandLineParser::Int, "Seeds per voxel:", "number of seed points per voxel", 1); parser.addArgument("seed_image", "", mitkDiffusionCommandLineParser::String, "Seed image:", "mask image defining seed voxels", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("trials_per_seed", "", mitkDiffusionCommandLineParser::Int, "Max. trials per seed:", "try each seed N times until a valid streamline is obtained (only for probabilistic tractography)", 10); parser.addArgument("max_tracts", "", mitkDiffusionCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded", -1); parser.endGroup(); parser.beginGroup("3. Tractography constraints:"); parser.addArgument("tracking_mask", "", mitkDiffusionCommandLineParser::String, "Mask image:", "streamlines leaving the mask will stop immediately", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("stop_image", "", mitkDiffusionCommandLineParser::String, "Stop ROI image:", "streamlines entering the mask will stop immediately", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("exclusion_image", "", mitkDiffusionCommandLineParser::String, "Exclusion ROI image:", "streamlines entering the mask will be discarded", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("ep_constraint", "", mitkDiffusionCommandLineParser::String, "Endpoint constraint:", "determines which fibers are accepted based on their endpoint location - options are NONE, EPS_IN_TARGET, EPS_IN_TARGET_LABELDIFF, EPS_IN_SEED_AND_TARGET, MIN_ONE_EP_IN_TARGET, ONE_EP_IN_TARGET and NO_EP_IN_TARGET", us::Any()); parser.addArgument("target_image", "", mitkDiffusionCommandLineParser::String, "Target ROI image:", "effact depends on the chosen endpoint constraint (option ep_constraint)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.endGroup(); parser.beginGroup("4. Streamline integration parameters:"); parser.addArgument("sharpen_odfs", "", mitkDiffusionCommandLineParser::Bool, "SHarpen ODFs:", "if you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). this is not necessary for CSD fODFs, since they are narurally much sharper."); parser.addArgument("cutoff", "", mitkDiffusionCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 0.1); parser.addArgument("odf_cutoff", "", mitkDiffusionCommandLineParser::Float, "ODF Cutoff:", "threshold on the ODF magnitude. this is useful in case of CSD fODF tractography.", 0.0); parser.addArgument("step_size", "", mitkDiffusionCommandLineParser::Float, "Step size:", "step size (in voxels)", 0.5); parser.addArgument("min_tract_length", "", mitkDiffusionCommandLineParser::Float, "Min. tract length:", "minimum fiber length (in mm)", 20); parser.addArgument("angular_threshold", "", mitkDiffusionCommandLineParser::Float, "Angular threshold:", "angular threshold between two successive steps, (default: 90° * step_size, minimum 15°)"); parser.addArgument("loop_check", "", mitkDiffusionCommandLineParser::Float, "Check for loops:", "threshold on angular stdev over the last 4 voxel lengths"); parser.addArgument("peak_jitter", "", mitkDiffusionCommandLineParser::Float, "Peak jitter:", "important for probabilistic peak tractography and peak prior. actual jitter is drawn from a normal distribution with peak_jitter*fabs(direction_value) as standard deviation.", 0.01); parser.endGroup(); parser.beginGroup("5. Tractography prior:"); parser.addArgument("prior_image", "", mitkDiffusionCommandLineParser::String, "Peak prior:", "tractography prior in thr for of a peak image", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("prior_weight", "", mitkDiffusionCommandLineParser::Float, "Prior weight", "weighting factor between prior and data.", 0.5); parser.addArgument("dont_restrict_to_prior", "", mitkDiffusionCommandLineParser::Bool, "Don't restrict to prior:", "don't restrict tractography to regions where the prior is valid.", us::Any(false)); parser.addArgument("no_new_directions_from_prior", "", mitkDiffusionCommandLineParser::Bool, "No new directios from prior:", "the prior cannot create directions where there are none in the data.", us::Any(false)); parser.addArgument("prior_flip_x", "", mitkDiffusionCommandLineParser::Bool, "Prior Flip X:", "multiply x-coordinate of prior direction by -1"); parser.addArgument("prior_flip_y", "", mitkDiffusionCommandLineParser::Bool, "Prior Flip Y:", "multiply y-coordinate of prior direction by -1"); parser.addArgument("prior_flip_z", "", mitkDiffusionCommandLineParser::Bool, "Prior Flip Z:", "multiply z-coordinate of prior direction by -1"); parser.endGroup(); parser.beginGroup("6. Neighborhood sampling:"); parser.addArgument("num_samples", "", mitkDiffusionCommandLineParser::Int, "Num. neighborhood samples:", "number of neighborhood samples that are use to determine the next progression direction", 0); parser.addArgument("sampling_distance", "", mitkDiffusionCommandLineParser::Float, "Sampling distance:", "distance of neighborhood sampling points (in voxels)", 0.25); parser.addArgument("use_stop_votes", "", mitkDiffusionCommandLineParser::Bool, "Use stop votes:", "use stop votes"); parser.addArgument("use_only_forward_samples", "", mitkDiffusionCommandLineParser::Bool, "Use only forward samples:", "use only forward samples"); parser.endGroup(); parser.beginGroup("7. Tensor tractography specific:"); parser.addArgument("tend_f", "", mitkDiffusionCommandLineParser::Float, "Weight f", "weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0); parser.addArgument("tend_g", "", mitkDiffusionCommandLineParser::Float, "Weight g", "weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0); parser.endGroup(); parser.beginGroup("8. Random forest tractography specific:"); parser.addArgument("forest", "", mitkDiffusionCommandLineParser::String, "Forest:", "input random forest (HDF5 file)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("use_sh_features", "", mitkDiffusionCommandLineParser::Bool, "Use SH features:", "use SH features"); parser.endGroup(); parser.beginGroup("9. Additional input:"); parser.addArgument("additional_images", "", mitkDiffusionCommandLineParser::StringList, "Additional images:", "specify a list of float images that hold additional information (FA, GFA, additional features for RF tractography)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.endGroup(); parser.beginGroup("10. Misc:"); parser.addArgument("flip_x", "", mitkDiffusionCommandLineParser::Bool, "Flip X:", "multiply x-coordinate of direction proposal by -1"); parser.addArgument("flip_y", "", mitkDiffusionCommandLineParser::Bool, "Flip Y:", "multiply y-coordinate of direction proposal by -1"); parser.addArgument("flip_z", "", mitkDiffusionCommandLineParser::Bool, "Flip Z:", "multiply z-coordinate of direction proposal by -1"); parser.addArgument("no_data_interpolation", "", mitkDiffusionCommandLineParser::Bool, "Don't interpolate input data:", "don't interpolate input image values"); parser.addArgument("no_mask_interpolation", "", mitkDiffusionCommandLineParser::Bool, "Don't interpolate masks:", "don't interpolate mask image values"); parser.addArgument("compress", "", mitkDiffusionCommandLineParser::Float, "Compress:", "compress output fibers using the given error threshold (in mm)"); parser.addArgument("fix_seed", "", mitkDiffusionCommandLineParser::Bool, "Fix Random Seed:", "always use the same random numbers"); parser.addArgument("parameter_file", "", mitkDiffusionCommandLineParser::String, "Parameter File:", "load parameters from json file (svae using MITK Diffusion GUI). the parameters loaded form this file are overwritten by the manually set parameters.", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["i"]); std::string outFile = us::any_cast(parsedArgs["o"]); std::string type = us::any_cast(parsedArgs["type"]); std::shared_ptr< mitk::StreamlineTractographyParameters > params = std::make_shared(); if (parsedArgs.count("parameter_file")) { auto parameter_file = us::any_cast(parsedArgs["parameter_file"]); params->LoadParameters(parameter_file); } if (parsedArgs.count("probabilistic")) params->m_Mode = mitk::StreamlineTractographyParameters::MODE::PROBABILISTIC; else { params->m_Mode = mitk::StreamlineTractographyParameters::MODE::DETERMINISTIC; } std::string prior_image = ""; if (parsedArgs.count("prior_image")) prior_image = us::any_cast(parsedArgs["prior_image"]); if (parsedArgs.count("prior_weight")) params->m_Weight = us::any_cast(parsedArgs["prior_weight"]); if (parsedArgs.count("fix_seed")) params->m_FixRandomSeed = us::any_cast(parsedArgs["fix_seed"]); params->m_RestrictToPrior = true; if (parsedArgs.count("dont_restrict_to_prior")) params->m_RestrictToPrior = !us::any_cast(parsedArgs["dont_restrict_to_prior"]); params->m_NewDirectionsFromPrior = true; if (parsedArgs.count("no_new_directions_from_prior")) params->m_NewDirectionsFromPrior = !us::any_cast(parsedArgs["no_new_directions_from_prior"]); params->m_SharpenOdfs = false; if (parsedArgs.count("sharpen_odfs")) params->m_SharpenOdfs = us::any_cast(parsedArgs["sharpen_odfs"]); params->m_InterpolateTractographyData = true; if (parsedArgs.count("no_data_interpolation")) params->m_InterpolateTractographyData = !us::any_cast(parsedArgs["no_data_interpolation"]); params->m_InterpolateRoiImages = true; if (parsedArgs.count("no_mask_interpolation")) params->m_InterpolateRoiImages = !us::any_cast(parsedArgs["no_mask_interpolation"]); bool use_sh_features = false; if (parsedArgs.count("use_sh_features")) use_sh_features = us::any_cast(parsedArgs["use_sh_features"]); params->m_StopVotes = false; if (parsedArgs.count("use_stop_votes")) params->m_StopVotes = us::any_cast(parsedArgs["use_stop_votes"]); params->m_OnlyForwardSamples = false; if (parsedArgs.count("use_only_forward_samples")) params->m_OnlyForwardSamples = us::any_cast(parsedArgs["use_only_forward_samples"]); params->m_FlipX = false; if (parsedArgs.count("flip_x")) params->m_FlipX = us::any_cast(parsedArgs["flip_x"]); params->m_FlipY = false; if (parsedArgs.count("flip_y")) params->m_FlipY = us::any_cast(parsedArgs["flip_y"]); params->m_FlipZ = false; if (parsedArgs.count("flip_z")) params->m_FlipZ = us::any_cast(parsedArgs["flip_z"]); bool prior_flip_x = false; if (parsedArgs.count("prior_flip_x")) prior_flip_x = us::any_cast(parsedArgs["prior_flip_x"]); bool prior_flip_y = false; if (parsedArgs.count("prior_flip_y")) prior_flip_y = us::any_cast(parsedArgs["prior_flip_y"]); bool prior_flip_z = false; if (parsedArgs.count("prior_flip_z")) prior_flip_z = us::any_cast(parsedArgs["prior_flip_z"]); params->m_ApplyDirectionMatrix = false; if (parsedArgs.count("apply_image_rotation")) params->m_ApplyDirectionMatrix = us::any_cast(parsedArgs["apply_image_rotation"]); float compress = -1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); params->m_MinTractLengthMm = 20; if (parsedArgs.count("min_tract_length")) params->m_MinTractLengthMm = us::any_cast(parsedArgs["min_tract_length"]); params->SetLoopCheckDeg(-1); if (parsedArgs.count("loop_check")) params->SetLoopCheckDeg(us::any_cast(parsedArgs["loop_check"])); std::string forestFile; if (parsedArgs.count("forest")) forestFile = us::any_cast(parsedArgs["forest"]); std::string maskFile = ""; if (parsedArgs.count("tracking_mask")) maskFile = us::any_cast(parsedArgs["tracking_mask"]); std::string seedFile = ""; if (parsedArgs.count("seed_image")) seedFile = us::any_cast(parsedArgs["seed_image"]); std::string targetFile = ""; if (parsedArgs.count("target_image")) targetFile = us::any_cast(parsedArgs["target_image"]); std::string exclusionFile = ""; if (parsedArgs.count("exclusion_image")) exclusionFile = us::any_cast(parsedArgs["exclusion_image"]); std::string stopFile = ""; if (parsedArgs.count("stop_image")) stopFile = us::any_cast(parsedArgs["stop_image"]); std::string ep_constraint = "NONE"; if (parsedArgs.count("ep_constraint")) ep_constraint = us::any_cast(parsedArgs["ep_constraint"]); params->m_Cutoff = 0.1f; if (parsedArgs.count("cutoff")) params->m_Cutoff = us::any_cast(parsedArgs["cutoff"]); params->m_OdfCutoff = 0.0; if (parsedArgs.count("odf_cutoff")) params->m_OdfCutoff = us::any_cast(parsedArgs["odf_cutoff"]); params->m_PeakJitter = 0.01; if (parsedArgs.count("peak_jitter")) params->m_PeakJitter = us::any_cast(parsedArgs["peak_jitter"]); params->SetStepSizeVox(-1); if (parsedArgs.count("step_size")) params->SetStepSizeVox(us::any_cast(parsedArgs["step_size"])); params->SetSamplingDistanceVox(-1); if (parsedArgs.count("sampling_distance")) params->SetSamplingDistanceVox(us::any_cast(parsedArgs["sampling_distance"])); params->m_NumSamples = 0; if (parsedArgs.count("num_samples")) params->m_NumSamples = static_cast(us::any_cast(parsedArgs["num_samples"])); params->m_SeedsPerVoxel = 1; if (parsedArgs.count("seeds")) params->m_SeedsPerVoxel = us::any_cast(parsedArgs["seeds"]); params->m_TrialsPerSeed = 10; if (parsedArgs.count("trials_per_seed")) params->m_TrialsPerSeed = static_cast(us::any_cast(parsedArgs["trials_per_seed"])); params->m_F = 1; if (parsedArgs.count("tend_f")) params->m_F = us::any_cast(parsedArgs["tend_f"]); params->m_G = 0; if (parsedArgs.count("tend_g")) params->m_G = us::any_cast(parsedArgs["tend_g"]); params->SetAngularThresholdDeg(-1); if (parsedArgs.count("angular_threshold")) params->SetAngularThresholdDeg(us::any_cast(parsedArgs["angular_threshold"])); params->m_MaxNumFibers = -1; if (parsedArgs.count("max_tracts")) params->m_MaxNumFibers = us::any_cast(parsedArgs["max_tracts"]); std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); if (ext != ".fib" && ext != ".trk") { MITK_INFO << "Output file format not supported. Use one of .fib, .trk, .nii, .nii.gz, .nrrd"; return EXIT_FAILURE; } // LOAD DATASETS mitkDiffusionCommandLineParser::StringContainerType addFiles; if (parsedArgs.count("additional_images")) addFiles = us::any_cast(parsedArgs["additional_images"]); typedef itk::Image ItkFloatImgType; ItkFloatImgType::Pointer mask = nullptr; if (!maskFile.empty()) { MITK_INFO << "loading mask image"; mitk::Image::Pointer img = mitk::IOUtil::Load(maskFile); mask = ItkFloatImgType::New(); mitk::CastToItkImage(img, mask); } ItkFloatImgType::Pointer seed = nullptr; if (!seedFile.empty()) { MITK_INFO << "loading seed ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(seedFile); seed = ItkFloatImgType::New(); mitk::CastToItkImage(img, seed); } ItkFloatImgType::Pointer stop = nullptr; if (!stopFile.empty()) { MITK_INFO << "loading stop ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(stopFile); stop = ItkFloatImgType::New(); mitk::CastToItkImage(img, stop); } ItkFloatImgType::Pointer target = nullptr; if (!targetFile.empty()) { MITK_INFO << "loading target ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(targetFile); target = ItkFloatImgType::New(); mitk::CastToItkImage(img, target); } ItkFloatImgType::Pointer exclusion = nullptr; if (!exclusionFile.empty()) { MITK_INFO << "loading exclusion ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(exclusionFile); exclusion = ItkFloatImgType::New(); mitk::CastToItkImage(img, exclusion); } MITK_INFO << "loading additional images"; std::vector< std::vector< ItkFloatImgType::Pointer > > addImages; addImages.push_back(std::vector< ItkFloatImgType::Pointer >()); for (auto file : addFiles) { mitk::Image::Pointer img = mitk::IOUtil::Load(file); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); addImages.at(0).push_back(itkimg); } // ////////////////////////////////////////////////////////////////// // omp_set_num_threads(1); typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); if (!prior_image.empty()) { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, std::vector()); mitk::PeakImage::Pointer priorImage = mitk::IOUtil::Load(prior_image, &functor); if (priorImage.IsNull()) { MITK_INFO << "Only peak images are supported as prior at the moment!"; return EXIT_FAILURE; } mitk::TrackingDataHandler* priorhandler = new mitk::TrackingHandlerPeaks(); typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(priorImage); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); std::shared_ptr< mitk::StreamlineTractographyParameters > prior_params = std::make_shared< mitk::StreamlineTractographyParameters >(*params); prior_params->m_FlipX = prior_flip_x; prior_params->m_FlipY = prior_flip_y; prior_params->m_FlipZ = prior_flip_z; prior_params->m_Cutoff = 0.0; dynamic_cast(priorhandler)->SetPeakImage(itkImg); priorhandler->SetParameters(prior_params); tracker->SetTrackingPriorHandler(priorhandler); } mitk::TrackingDataHandler* handler; mitk::Image::Pointer reference_image; if (type == "RF") { mitk::TractographyForest::Pointer forest = mitk::IOUtil::Load(forestFile); if (forest.IsNull()) mitkThrow() << "Forest file " << forestFile << " could not be read."; - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); + std::vector include = {"Diffusion Weighted Images"}; + std::vector exclude = {}; + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor(include, exclude); auto input = mitk::IOUtil::Load(input_files.at(0), &functor); reference_image = input; if (use_sh_features) { handler = new mitk::TrackingHandlerRandomForest<6,28>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } else { handler = new mitk::TrackingHandlerRandomForest<6,100>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } } else if (type == "Peaks") { handler = new mitk::TrackingHandlerPeaks(); MITK_INFO << "loading input peak image"; mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(input_files.at(0)); reference_image = mitkImage; mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = mitk::convert::GetItkPeakFromPeakImage(mitkImage); dynamic_cast(handler)->SetPeakImage(itkImg); } else if (type == "Tensor" && params->m_Mode == mitk::StreamlineTractographyParameters::MODE::DETERMINISTIC) { handler = new mitk::TrackingHandlerTensor(); MITK_INFO << "loading input tensor images"; std::vector< mitk::Image::Pointer > input_images; for (unsigned int i=0; i(input_files.at(i)); reference_image = mitkImage; mitk::TensorImage::ItkTensorImageType::Pointer itkImg = mitk::convert::GetItkTensorFromTensorImage(mitkImage); dynamic_cast(handler)->AddTensorImage(itkImg.GetPointer()); } if (addImages.at(0).size()>0) dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); } else if (type == "ODF" || type == "ODF-DIPY/FSL" || (type == "Tensor" && params->m_Mode == mitk::StreamlineTractographyParameters::MODE::PROBABILISTIC)) { handler = new mitk::TrackingHandlerOdf(); mitk::OdfImage::ItkOdfImageType::Pointer itkImg = nullptr; if (type == "Tensor") { MITK_INFO << "Converting Tensor to ODF image"; auto input = mitk::IOUtil::Load(input_files.at(0)); reference_image = input; itkImg = mitk::convert::GetItkOdfFromTensorImage(input); dynamic_cast(handler)->SetIsOdfFromTensor(true); } else { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image", "ODF Image"}, {}); + std::vector include = {"SH Image", "ODF Image"}; + std::vector exclude = {}; + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor(include, exclude); auto input = mitk::IOUtil::Load(input_files.at(0), &functor)[0]; reference_image = dynamic_cast(input.GetPointer()); if (dynamic_cast(input.GetPointer())) { MITK_INFO << "Converting SH to ODF image"; mitk::ShImage::Pointer mitkShImage = dynamic_cast(input.GetPointer()); if (type == "ODF-DIPY/FSL") mitkShImage->SetShConvention(mitk::ShImage::SH_CONVENTION::FSL); mitk::Image::Pointer mitkImg = dynamic_cast(mitkShImage.GetPointer()); itkImg = mitk::convert::GetItkOdfFromShImage(mitkImg); } else if (dynamic_cast(input.GetPointer())) { mitk::Image::Pointer mitkImg = dynamic_cast(input.GetPointer()); itkImg = mitk::convert::GetItkOdfFromOdfImage(mitkImg); } else mitkThrow() << ""; } dynamic_cast(handler)->SetOdfImage(itkImg); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); } else { MITK_INFO << "Unknown tractography algorithm (" + type+"). Known types are Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF."; return EXIT_FAILURE; } if (ep_constraint=="NONE") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::NONE; else if (ep_constraint=="EPS_IN_TARGET") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET; else if (ep_constraint=="EPS_IN_TARGET_LABELDIFF") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF; else if (ep_constraint=="EPS_IN_SEED_AND_TARGET") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET; else if (ep_constraint=="MIN_ONE_EP_IN_TARGET") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET; else if (ep_constraint=="ONE_EP_IN_TARGET") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET; else if (ep_constraint=="NO_EP_IN_TARGET") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET; MITK_INFO << "Tractography algorithm: " << type; tracker->SetMaskImage(mask); tracker->SetSeedImage(seed); tracker->SetStoppingRegions(stop); tracker->SetTargetRegions(target); tracker->SetExclusionRegions(exclusion); tracker->SetTrackingHandler(handler); if (ext != ".fib" && ext != ".trk") params->m_OutputProbMap = true; tracker->SetParameters(params); tracker->Update(); if (ext == ".fib" || ext == ".trk") { vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); if (compress > 0) outFib->Compress(compress); outFib->SetTrackVisHeader(reference_image->GetGeometry()); mitk::IOUtil::Save(outFib, outFile); } else { TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (ext != ".nii" && ext != ".nii.gz" && ext != ".nrrd") outFile += ".nii.gz"; mitk::IOUtil::Save(img, outFile); } delete handler; return EXIT_SUCCESS; } diff --git a/Modules/FiberTracking/Testing/mitkPeakShImageReaderTest.cpp b/Modules/FiberTracking/Testing/mitkPeakShImageReaderTest.cpp index 7d5b9af..24fa53e 100644 --- a/Modules/FiberTracking/Testing/mitkPeakShImageReaderTest.cpp +++ b/Modules/FiberTracking/Testing/mitkPeakShImageReaderTest.cpp @@ -1,127 +1,127 @@ /*=================================================================== 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 #include #include #include #include "mitkTestingMacros.h" #include #include #include #include #include #include #include #include #include class mitkPeakShImageReaderTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPeakShImageReaderTestSuite); MITK_TEST(PeakReader); MITK_TEST(ShReader); CPPUNIT_TEST_SUITE_END(); typedef itk::Image ItkFloatImgType; private: /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ public: void setUp() override { } void tearDown() override { } void PeakReader() { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, std::vector()); auto inputData = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/peak_image_test.nii.gz"), &functor); std::string class_name = "PeakImage"; bool ok = true; if (class_name.compare(inputData->GetNameOfClass())!=0) ok = false; MITK_TEST_CONDITION_REQUIRED(ok, "Check class name = PeakImage"); mitk::Vector3D sp; sp[0] = 1.23; sp[1] = 2.001; sp[2] = 1.409; MITK_TEST_CONDITION_REQUIRED(mitk::Equal(sp, inputData->GetGeometry()->GetSpacing(), 0.0001, true), "Check spacing ok"); mitk::Point3D o; o[0] = 0.615; o[1] = 1.0005; o[2] = 0.7045; MITK_TEST_CONDITION_REQUIRED(mitk::Equal(o, inputData->GetGeometry()->GetOrigin(), 0.0001, true), "Check origin ok"); mitk::ImagePixelReadAccessor readAccess(inputData); itk::Index<4> idx; idx[0] = 5; idx[1] = 4; idx[2] = 1; idx[3] = 0; MITK_TEST_CONDITION_REQUIRED(fabs(-0.997147f-readAccess.GetPixelByIndex(idx))<0.0001f, "Check pixel value ok"); } void ShReader() { - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, std::vector()); auto inputData = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/sh_image_test.nii.gz"), &functor); std::string class_name = "ShImage"; bool ok = true; if (class_name.compare(inputData->GetNameOfClass())!=0) ok = false; MITK_TEST_CONDITION_REQUIRED(ok, "Check class name = ShImage"); mitk::Vector3D sp; sp[0] = 1.23; sp[1] = 2.001; sp[2] = 1.409; MITK_TEST_CONDITION_REQUIRED(mitk::Equal(sp, inputData->GetGeometry()->GetSpacing(), 0.0001, true), "Check spacing ok"); mitk::Point3D o; o[0] = 0.615; o[1] = 1.0005; o[2] = 0.7045; MITK_TEST_CONDITION_REQUIRED(mitk::Equal(o, inputData->GetGeometry()->GetOrigin(), 0.0001, true), "Check origin ok"); mitk::ImagePixelReadAccessor, 3> readAccess(inputData); itk::Index<3> idx; idx[0] = 5; idx[1] = 4; idx[2] = 1; itk::Vector pixel = readAccess.GetPixelByIndex(idx); MITK_TEST_CONDITION_REQUIRED(fabs(0.282095f-pixel[0])<0.000001f, "Check SH coeff value 0 ok"); MITK_TEST_CONDITION_REQUIRED(fabs(-0.00900571f-pixel[2])<0.000001f, "Check SH coeff value 2 ok"); MITK_TEST_CONDITION_REQUIRED(fabs(0.0513185f-pixel[6])<0.000001f, "Check SH coeff value 6 ok"); } }; MITK_TEST_SUITE_REGISTRATION(mitkPeakShImageReader) diff --git a/Modules/MriSimulation/Testing/mitkFiberFitTest.cpp b/Modules/MriSimulation/Testing/mitkFiberFitTest.cpp index b49846b..d2c6014 100644 --- a/Modules/MriSimulation/Testing/mitkFiberFitTest.cpp +++ b/Modules/MriSimulation/Testing/mitkFiberFitTest.cpp @@ -1,277 +1,277 @@ /*=================================================================== 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 "mitkTestingMacros.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include class mitkFiberFitTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkFiberFitTestSuite); MITK_TEST(Fit1); MITK_TEST(Fit2); MITK_TEST(Fit3); MITK_TEST(Fit4); MITK_TEST(Fit5); MITK_TEST(Fit6); CPPUNIT_TEST_SUITE_END(); typedef itk::Image ItkFloatImgType; private: /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ typedef itk::FitFibersToImageFilter FitterType; FitterType::Pointer fitter; public: mitk::FiberBundle::Pointer LoadFib(std::string fib_name) { std::vector fibInfile = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberFit/" + fib_name)); mitk::BaseData::Pointer baseData = fibInfile.at(0); mitk::FiberBundle::Pointer fib = dynamic_cast(baseData.GetPointer()); return fib; } void setUp() override { std::vector tracts; tracts.push_back(LoadFib("Cluster_0.fib")); tracts.push_back(LoadFib("Cluster_1.fib")); tracts.push_back(LoadFib("Cluster_2.fib")); tracts.push_back(LoadFib("Cluster_3.fib")); tracts.push_back(LoadFib("Cluster_4.fib")); - mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, {}); + mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, std::vector()); mitk::PeakImage::Pointer peaks = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberFit/csd_peak_image.nii.gz"), &functor); typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(peaks); caster->Update(); mitk::PeakImage::ItkPeakImageType::Pointer peak_image = caster->GetOutput(); fitter = FitterType::New(); fitter->SetPeakImage(peak_image); fitter->SetTractograms(tracts); } void tearDown() override { } void CompareFibs(mitk::FiberBundle::Pointer test, mitk::FiberBundle::Pointer ref, std::string out_name) { vtkSmartPointer weights = test->GetFiberWeights(); vtkSmartPointer ref_weights = ref->GetFiberWeights(); CPPUNIT_ASSERT_MESSAGE("Number of weights should be equal", weights->GetSize()==ref_weights->GetSize()); for (int i=0; iGetSize(); ++i) { if (ref_weights->GetValue(i)>0) { if (fabs( weights->GetValue(i)/ref_weights->GetValue(i)-1 )>0.01) { mitk::IOUtil::Save(test, mitk::IOUtil::GetTempPath()+out_name); CPPUNIT_ASSERT_MESSAGE("Weights should be equal", false); } } else if (weights->GetValue(i)>0) { mitk::IOUtil::Save(test, mitk::IOUtil::GetTempPath()+out_name); CPPUNIT_ASSERT_MESSAGE("Weights should be equal", false); } } } void CompareImages(mitk::PeakImage::ItkPeakImageType::Pointer testImage, std::string name) { mitk::LocaleSwitch localeSwitch("C"); typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberFit/out/" + name))); caster->Update(); mitk::PeakImage::ItkPeakImageType::Pointer refImage = caster->GetOutput(); itk::ImageRegionConstIterator< mitk::PeakImage::ItkPeakImageType > it1(testImage, testImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator< mitk::PeakImage::ItkPeakImageType > it2(refImage, refImage->GetLargestPossibleRegion()); while(!it1.IsAtEnd()) { if (it2.Get()>0.0001) { if (fabs( it1.Get()/it2.Get()-1 )>0.01) { itk::ImageFileWriter< mitk::PeakImage::ItkPeakImageType >::Pointer writer = itk::ImageFileWriter< mitk::PeakImage::ItkPeakImageType >::New(); writer->SetInput(testImage); writer->SetFileName(mitk::IOUtil::GetTempPath()+name); writer->Update(); MITK_INFO << it1.Get() << " - " << it2.Get(); CPPUNIT_ASSERT_MESSAGE("Peak images should be equal 1", false); } } else if (it1.Get()>0.0001) { itk::ImageFileWriter< mitk::PeakImage::ItkPeakImageType >::Pointer writer = itk::ImageFileWriter< mitk::PeakImage::ItkPeakImageType >::New(); writer->SetInput(testImage); writer->SetFileName(mitk::IOUtil::GetTempPath()+name); writer->Update(); CPPUNIT_ASSERT_MESSAGE("Peak images should be equal 2", false); } ++it1; ++it2; } } void Fit1() { omp_set_num_threads(1); // fitter->SetLambda(0.1); fitter->SetFilterOutliers(false); fitter->SetRegularization(VnlCostFunction::NONE); fitter->Update(); std::vector< mitk::FiberBundle::Pointer > output_tracts = fitter->GetTractograms(); mitk::FiberBundle::Pointer test = mitk::FiberBundle::New(); test = test->AddBundles(output_tracts); mitk::FiberBundle::Pointer ref = LoadFib("out/NONE_fitted.fib"); CompareFibs(test, ref, "NONE_fitted.fib"); CompareImages(fitter->GetFittedImage(), "NONE_fitted_image.nrrd"); CompareImages(fitter->GetResidualImage(), "NONE_residual_image.nrrd"); CompareImages(fitter->GetOverexplainedImage(), "NONE_underexplained_image.nrrd"); CompareImages(fitter->GetUnderexplainedImage(), "NONE_overexplained_image.nrrd"); } void Fit2() { omp_set_num_threads(1); // fitter->SetLambda(0.1); fitter->SetFilterOutliers(false); fitter->SetRegularization(VnlCostFunction::MSM); fitter->Update(); std::vector< mitk::FiberBundle::Pointer > output_tracts = fitter->GetTractograms(); mitk::FiberBundle::Pointer test = mitk::FiberBundle::New(); test = test->AddBundles(output_tracts); mitk::FiberBundle::Pointer ref = LoadFib("out/MSM_fitted.fib"); CompareFibs(test, ref, "MSM_fitted.fib"); CompareImages(fitter->GetFittedImage(), "MSM_fitted_image.nrrd"); CompareImages(fitter->GetResidualImage(), "MSM_residual_image.nrrd"); CompareImages(fitter->GetOverexplainedImage(), "MSM_underexplained_image.nrrd"); CompareImages(fitter->GetUnderexplainedImage(), "MSM_overexplained_image.nrrd"); } void Fit3() { omp_set_num_threads(1); // fitter->SetLambda(0.1); fitter->SetFilterOutliers(false); fitter->SetRegularization(VnlCostFunction::VARIANCE); fitter->Update(); std::vector< mitk::FiberBundle::Pointer > output_tracts = fitter->GetTractograms(); mitk::FiberBundle::Pointer test = mitk::FiberBundle::New(); test = test->AddBundles(output_tracts); mitk::FiberBundle::Pointer ref = LoadFib("out/MSE_fitted.fib"); CompareFibs(test, ref, "MSE_fitted.fib"); CompareImages(fitter->GetFittedImage(), "MSE_fitted_image.nrrd"); CompareImages(fitter->GetResidualImage(), "MSE_residual_image.nrrd"); CompareImages(fitter->GetOverexplainedImage(), "MSE_underexplained_image.nrrd"); CompareImages(fitter->GetUnderexplainedImage(), "MSE_overexplained_image.nrrd"); } void Fit4() { omp_set_num_threads(1); fitter->SetLambda(100); fitter->SetFilterOutliers(false); fitter->SetRegularization(VnlCostFunction::VOXEL_VARIANCE); fitter->Update(); std::vector< mitk::FiberBundle::Pointer > output_tracts = fitter->GetTractograms(); mitk::FiberBundle::Pointer test = mitk::FiberBundle::New(); test = test->AddBundles(output_tracts); mitk::FiberBundle::Pointer ref = LoadFib("out/LocalMSE_fitted.fib"); CompareFibs(test, ref, "LocalMSE_fitted.fib"); CompareImages(fitter->GetFittedImage(), "LocalMSE_fitted_image.nrrd"); CompareImages(fitter->GetResidualImage(), "LocalMSE_residual_image.nrrd"); CompareImages(fitter->GetOverexplainedImage(), "LocalMSE_underexplained_image.nrrd"); CompareImages(fitter->GetUnderexplainedImage(), "LocalMSE_overexplained_image.nrrd"); } void Fit5() { omp_set_num_threads(1); fitter->SetFilterOutliers(false); fitter->SetRegularization(VnlCostFunction::GROUP_VARIANCE); fitter->Update(); std::vector< mitk::FiberBundle::Pointer > output_tracts = fitter->GetTractograms(); mitk::FiberBundle::Pointer test = mitk::FiberBundle::New(); test = test->AddBundles(output_tracts); mitk::FiberBundle::Pointer ref = LoadFib("out/GroupMSE_fitted.fib"); CompareFibs(test, ref, "GroupMSE_fitted.fib"); CompareImages(fitter->GetFittedImage(), "GroupMSE_fitted_image.nrrd"); CompareImages(fitter->GetResidualImage(), "GroupMSE_residual_image.nrrd"); CompareImages(fitter->GetOverexplainedImage(), "GroupMSE_underexplained_image.nrrd"); CompareImages(fitter->GetUnderexplainedImage(), "GroupMSE_overexplained_image.nrrd"); } void Fit6() { omp_set_num_threads(1); fitter->SetLambda(1000); fitter->SetFilterOutliers(false); fitter->SetRegularization(VnlCostFunction::GROUP_LASSO); fitter->Update(); std::vector< mitk::FiberBundle::Pointer > output_tracts = fitter->GetTractograms(); mitk::FiberBundle::Pointer test = mitk::FiberBundle::New(); test = test->AddBundles(output_tracts); mitk::FiberBundle::Pointer ref = LoadFib("out/GroupLasso_fitted.fib"); CompareFibs(test, ref, "GroupLasso_fitted.fib"); CompareImages(fitter->GetFittedImage(), "GroupLasso_fitted_image.nrrd"); CompareImages(fitter->GetResidualImage(), "GroupLasso_residual_image.nrrd"); CompareImages(fitter->GetOverexplainedImage(), "GroupLasso_underexplained_image.nrrd"); CompareImages(fitter->GetUnderexplainedImage(), "GroupLasso_overexplained_image.nrrd"); } }; MITK_TEST_SUITE_REGISTRATION(mitkFiberFit) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.partialvolume/CMakeLists.txt b/Plugins/org.mitk.gui.qt.diffusionimaging.partialvolume/CMakeLists.txt index 3aaaf92..155c715 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.partialvolume/CMakeLists.txt +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.partialvolume/CMakeLists.txt @@ -1,10 +1,11 @@ # The project name must correspond to the directory name of your plug-in # and must not contain periods. project(org_mitk_gui_qt_diffusionimaging_partialvolume) mitk_create_plugin( SUBPROJECTS MITK-Diffusion EXPORT_DIRECTIVE DIFFUSIONIMAGING_PARTIALVOLUME_EXPORT EXPORTED_INCLUDE_SUFFIXES src MODULE_DEPENDS MitkQtWidgetsExt MitkDiffusionCore MitkFiberTracking + WARNINGS_NO_ERRORS )