diff --git a/Modules/DiffusionCmdApps/FiberProcessing/DistanceToSegmentation.cpp b/Modules/DiffusionCmdApps/FiberProcessing/DistanceToSegmentation.cpp index 7b217d5..2fdd86f 100644 --- a/Modules/DiffusionCmdApps/FiberProcessing/DistanceToSegmentation.cpp +++ b/Modules/DiffusionCmdApps/FiberProcessing/DistanceToSegmentation.cpp @@ -1,169 +1,169 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include typedef itk::Image ItkFloatImgType; /*! \brief */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Calcualte distance between segmentation mesh and fibers in form of a mask or TDI"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("tracts", "", mitkDiffusionCommandLineParser::String, "Tracts:", "input fiber tracts (as tractogram or as tract density iumage (TDI))", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("surface", "", mitkDiffusionCommandLineParser::String, "Segmentation:", "input segmentation mesh (.vtp)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("out_dists", "", mitkDiffusionCommandLineParser::String, "Output:", "output text file", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("out_image", "", mitkDiffusionCommandLineParser::String, "Output:", "output image", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("thresholds", "", mitkDiffusionCommandLineParser::StringList, "Thresholds:", "distances for which voxels are counted", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string in_tracts = us::any_cast(parsedArgs["tracts"]); std::string in_seg = us::any_cast(parsedArgs["surface"]); std::string out_dists = ""; if (parsedArgs.count("out_dists")) out_dists = us::any_cast(parsedArgs["out_dists"]); std::string out_image = ""; if (parsedArgs.count("out_image")) out_image = us::any_cast(parsedArgs["out_image"]); std::vector< std::string > thresholds_str = {"0.0", "3.0", "5.0", "7.0"}; std::vector< float > thresholds; if (parsedArgs.count("thresholds")) thresholds_str = us::any_cast(parsedArgs["thresholds"]); for (auto s : thresholds_str) thresholds.push_back(boost::lexical_cast(s)); try { ItkFloatImgType::Pointer inputItkTDI; auto input= mitk::IOUtil::Load(in_tracts); if (dynamic_cast(input.GetPointer())) { mitk::CastToItkImage(dynamic_cast(input.GetPointer()), inputItkTDI); } else { mitk::FiberBundle::Pointer intput_tracts= dynamic_cast(input.GetPointer()); itk::TractDensityImageFilter< ItkFloatImgType >::Pointer tdi_filter = itk::TractDensityImageFilter< ItkFloatImgType >::New(); tdi_filter->SetFiberBundle(intput_tracts); - tdi_filter->SetBinaryOutput(false); + tdi_filter->SetMode(TDI_MODE::DENSITY); tdi_filter->SetOutputAbsoluteValues(false); tdi_filter->Update(); inputItkTDI = tdi_filter->GetOutput(); } mitk::Surface::Pointer inputSeg = mitk::IOUtil::Load(in_seg); typedef itk::DistanceFromSegmentationImageFilter< float > ImageGeneratorType; ImageGeneratorType::Pointer filter = ImageGeneratorType::New(); filter->SetInput(inputItkTDI); filter->SetSegmentationSurface(inputSeg); filter->SetThresholds(thresholds); filter->Update(); auto volume = inputItkTDI->GetSpacing()[0]*inputItkTDI->GetSpacing()[1]*inputItkTDI->GetSpacing()[2]; MITK_INFO << "Tractogram file: " << in_tracts; MITK_INFO << "Segmentation surface file: " << in_seg; MITK_INFO << "Reference image spacing: " << inputItkTDI->GetSpacing(); MITK_INFO << "Distance between surface and closest fiber-containing voxel-center: " << filter->GetMinDistance() << " mm"; for (unsigned int i=0; iGetThresholds().size(); ++i) MITK_INFO << "Voxels closer than " << filter->GetThresholds().at(i) << " mm: " << filter->GetCounts().at(i) << " (" << filter->GetCounts().at(i)*volume << " mm³)"; if (!out_dists.empty()) { MITK_INFO << "Saving data to text file: " << out_dists; std::ofstream statistics_file; statistics_file.open(out_dists, std::ios_base::out | std::ios_base::app); statistics_file << "MITK Diffusion git commit hash: " << MITKDIFFUSION_REVISION << std::endl; statistics_file << "MITK Diffusion branch name: " << MITKDIFFUSION_REVISION_NAME << std::endl; statistics_file << "MITK git commit hash: " << MITK_REVISION << std::endl; statistics_file << "MITK branch name: " << MITK_REVISION_NAME << std::endl; boost::posix_time::ptime timeLocal = boost::posix_time::second_clock::local_time(); statistics_file << "Current System Time = " << timeLocal << std::endl; statistics_file << "Tractogram file: " << in_tracts << std::endl; statistics_file << "Segmentation surface file: " << in_seg << std::endl; statistics_file << "Reference image spacing: " << inputItkTDI->GetSpacing() << std::endl; statistics_file << "Distance between surface and closest fiber-containing voxel-center: " << filter->GetMinDistance() << " mm" << std::endl; for (unsigned int i=0; iGetThresholds().size(); ++i) statistics_file << "Voxels closer than " << filter->GetThresholds().at(i) << " mm: " << filter->GetCounts().at(i) << " (" << filter->GetCounts().at(i)*volume << " mm³)" << std::endl; statistics_file << "-------------------------------------------------\n\n" << std::endl; statistics_file.close(); } // save image containing voxel-wise distances if (!out_image.empty()) { MITK_INFO << "Saving voxel-wise distances to image file: " << out_image; auto outImg = filter->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg); img->SetVolume(outImg->GetBufferPointer()); mitk::IOUtil::Save(img, out_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/FiberProcessing/TractDensity.cpp b/Modules/DiffusionCmdApps/FiberProcessing/TractDensity.cpp index e6b3853..3afd948 100644 --- a/Modules/DiffusionCmdApps/FiberProcessing/TractDensity.cpp +++ b/Modules/DiffusionCmdApps/FiberProcessing/TractDensity.cpp @@ -1,213 +1,213 @@ /*=================================================================== 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 "mitkDiffusionCommandLineParser.h" #include #include #include #include #include mitk::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } /*! \brief Modify input tractogram: fiber resampling, compression, pruning and transformation. */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Tract Density"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Generate tract density image, fiber envelope or fiber endpoints image."); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input fiber bundle", us::Any(), false); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output image", us::Any(), false); parser.addArgument("binary", "", mitkDiffusionCommandLineParser::Bool, "Binary output:", "calculate binary tract envelope", us::Any()); parser.addArgument("normalize", "", mitkDiffusionCommandLineParser::Bool, "Normalized output:", "normalize output to 0-1", us::Any()); parser.addArgument("endpoints", "", mitkDiffusionCommandLineParser::Bool, "Output endpoints image:", "calculate image of fiber endpoints instead of mask", us::Any()); parser.addArgument("reference_image", "", mitkDiffusionCommandLineParser::String, "Reference image:", "output image will have geometry of this reference image", us::Any()); parser.addArgument("upsampling", "", mitkDiffusionCommandLineParser::Float, "Upsampling:", "upsampling", 1.0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; bool binary = false; if (parsedArgs.count("binary")) binary = us::any_cast(parsedArgs["binary"]); bool endpoints = false; if (parsedArgs.count("endpoints")) endpoints = us::any_cast(parsedArgs["endpoints"]); bool normalize = false; if (parsedArgs.count("normalize")) normalize = us::any_cast(parsedArgs["normalize"]); float upsampling = 1.0; if (parsedArgs.count("upsampling")) upsampling = us::any_cast(parsedArgs["upsampling"]); MITK_INFO << "Upsampling: " << upsampling; std::string reference_image = ""; if (parsedArgs.count("reference_image")) reference_image = us::any_cast(parsedArgs["reference_image"]); std::string inFileName = us::any_cast(parsedArgs["i"]); std::string outFileName = us::any_cast(parsedArgs["o"]); try { mitk::FiberBundle::Pointer fib = LoadFib(inFileName); mitk::Image::Pointer ref_img; if (!reference_image.empty()) ref_img = mitk::IOUtil::Load(reference_image); if (endpoints) { typedef unsigned int OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetUpsamplingFactor(upsampling); if (ref_img.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(ref_img, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::IOUtil::Save(img, outFileName ); } else if (binary) { typedef unsigned char OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); - generator->SetBinaryOutput(binary); + generator->SetMode(TDI_MODE::BINARY); generator->SetOutputAbsoluteValues(!normalize); generator->SetUpsamplingFactor(upsampling); if (ref_img.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(ref_img, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::IOUtil::Save(img, outFileName ); } else { typedef float OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); - generator->SetBinaryOutput(binary); + generator->SetMode(TDI_MODE::DENSITY); generator->SetOutputAbsoluteValues(!normalize); generator->SetUpsamplingFactor(upsampling); if (ref_img.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(ref_img, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::IOUtil::Save(img, 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/FiberProcessing/TractDensityFilter.cpp b/Modules/DiffusionCmdApps/FiberProcessing/TractDensityFilter.cpp index 6271ee2..4912825 100644 --- a/Modules/DiffusionCmdApps/FiberProcessing/TractDensityFilter.cpp +++ b/Modules/DiffusionCmdApps/FiberProcessing/TractDensityFilter.cpp @@ -1,110 +1,110 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; /*! \brief Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Filter Outliers by Tract Density"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input tractogram (.fib/.trk/.tck/.dcm)", us::Any(), false); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output tractogram", us::Any(), false); parser.addArgument("threshold", "", mitkDiffusionCommandLineParser::Float, "Threshold:", "positive means ROI image value threshold", 0.05); parser.addArgument("overlap", "", mitkDiffusionCommandLineParser::Float, "Overlap:", "positive means ROI image value threshold", 0.5); parser.addArgument("min_fibers", "", mitkDiffusionCommandLineParser::Int, "Min. num. fibers:", "discard positive tracts with less fibers", 0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFib = us::any_cast(parsedArgs["i"]); std::string outFib = us::any_cast(parsedArgs["o"]); int min_fibers = 0; if (parsedArgs.count("min_fibers")) min_fibers = us::any_cast(parsedArgs["min_fibers"]); float overlap = 0.5; if (parsedArgs.count("overlap")) overlap = us::any_cast(parsedArgs["overlap"]); float threshold = 0.05f; if (parsedArgs.count("threshold")) threshold = us::any_cast(parsedArgs["threshold"]); try { mitk::FiberBundle::Pointer inputTractogram = mitk::IOUtil::Load(inFib); itk::TractDensityImageFilter< ItkFloatImgType >::Pointer generator = itk::TractDensityImageFilter< ItkFloatImgType >::New(); generator->SetFiberBundle(inputTractogram); - generator->SetBinaryOutput(false); + generator->SetMode(TDI_MODE::DENSITY); generator->SetOutputAbsoluteValues(false); generator->Update(); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetRoiImages({generator->GetOutput()}); extractor->SetInputFiberBundle(inputTractogram); extractor->SetOverlapFraction(overlap); extractor->SetInterpolate(true); extractor->SetThreshold(threshold); extractor->SetNoNegatives(true); extractor->Update(); if (extractor->GetPositives().at(0)->GetNumFibers() >= static_cast(min_fibers)) mitk::IOUtil::Save(extractor->GetPositives().at(0), outFib); } 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/TractographyEvaluation/AnchorConstrainedPlausibility.cpp b/Modules/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp index b227ea6..e2352de 100644 --- a/Modules/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp +++ b/Modules/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp @@ -1,447 +1,447 @@ /*=================================================================== 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 typedef itk::Point PointType4; typedef mitk::PeakImage::ItkPeakImageType PeakImgType; typedef itk::Image< unsigned char, 3 > ItkUcharImageType; /*! \brief Score input candidate tracts using ACP analysis */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Anchor Constrained Plausibility"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription("Score input candidate tracts using ACP analysis"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "a", mitkDiffusionCommandLineParser::String, "Anchor tractogram:", "anchor tracts in one tractogram file", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "p", mitkDiffusionCommandLineParser::String, "Input peaks:", "input peak image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "c", mitkDiffusionCommandLineParser::StringList, "Candidates:", "Folder(s) or file list of candidate tracts", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output folder:", "output folder", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("reference_mask_folders", "", mitkDiffusionCommandLineParser::StringList, "Reference Mask Folder(s):", "Folder(s) or file list containing reference tract masks for accuracy evaluation", true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("reference_peaks_folders", "", mitkDiffusionCommandLineParser::StringList, "Reference Peaks Folder(s):", "Folder(s) or file list containing reference peak images for accuracy evaluation", true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("mask", "", mitkDiffusionCommandLineParser::String, "Mask image:", "scoring is only performed inside the mask image", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("greedy_add", "", mitkDiffusionCommandLineParser::Bool, "Greedy:", "if enabled, the candidate tracts are not jointly fitted to the residual image but one after the other employing a greedy scheme", false); parser.addArgument("lambda", "", mitkDiffusionCommandLineParser::Float, "Lambda:", "modifier for regularization", 0.1); 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("regu", "", mitkDiffusionCommandLineParser::String, "Regularization:", "MSM; Variance; VoxelVariance; Lasso; GroupLasso; GroupVariance; NONE", std::string("NONE")); parser.addArgument("use_num_streamlines", "", mitkDiffusionCommandLineParser::Bool, "Use number of streamlines as score:", "Don't fit candidates, simply use number of streamlines per candidate as score", false); parser.addArgument("use_weights", "", mitkDiffusionCommandLineParser::Bool, "Use input weights as score:", "Don't fit candidates, simply use first input streamline weight per candidate as score", false); parser.addArgument("filter_zero_weights", "", mitkDiffusionCommandLineParser::Bool, "Filter zero-weights", "Remove streamlines with weight 0 from candidates", false); parser.addArgument("flipx", "", mitkDiffusionCommandLineParser::Bool, "Flip x", "flip along x-axis", false); parser.addArgument("flipy", "", mitkDiffusionCommandLineParser::Bool, "Flip y", "flip along y-axis", false); parser.addArgument("flipz", "", mitkDiffusionCommandLineParser::Bool, "Flip z", "flip along z-axis", false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string peak_file_name = us::any_cast(parsedArgs["p"]); std::string out_folder = us::any_cast(parsedArgs["o"]); mitkDiffusionCommandLineParser::StringContainerType candidate_tract_folders = us::any_cast(parsedArgs["c"]); if (!out_folder.empty() && out_folder.back() != '/') out_folder += "/"; bool greedy_add = false; if (parsedArgs.count("greedy_add")) greedy_add = us::any_cast(parsedArgs["greedy_add"]); float lambda = 0.1f; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); bool filter_outliers = false; if (parsedArgs.count("filter_outliers")) filter_outliers = us::any_cast(parsedArgs["filter_outliers"]); bool filter_zero_weights = false; if (parsedArgs.count("filter_zero_weights")) filter_zero_weights = us::any_cast(parsedArgs["filter_zero_weights"]); std::string mask_file = ""; if (parsedArgs.count("mask")) mask_file = us::any_cast(parsedArgs["mask"]); mitkDiffusionCommandLineParser::StringContainerType reference_mask_files_folders; if (parsedArgs.count("reference_mask_folders")) reference_mask_files_folders = us::any_cast(parsedArgs["reference_mask_folders"]); mitkDiffusionCommandLineParser::StringContainerType reference_peaks_files_folders; if (parsedArgs.count("reference_peaks_folders")) reference_peaks_files_folders = us::any_cast(parsedArgs["reference_peaks_folders"]); std::string regu = "NONE"; if (parsedArgs.count("regu")) regu = us::any_cast(parsedArgs["regu"]); bool use_weights = false; if (parsedArgs.count("use_weights")) use_weights = us::any_cast(parsedArgs["use_weights"]); bool use_num_streamlines = false; if (parsedArgs.count("use_num_streamlines")) use_num_streamlines = us::any_cast(parsedArgs["use_num_streamlines"]); bool flipx = false; if (parsedArgs.count("flipx")) flipx = us::any_cast(parsedArgs["flipx"]); bool flipy = false; if (parsedArgs.count("flipy")) flipy = us::any_cast(parsedArgs["flipy"]); bool flipz = false; if (parsedArgs.count("flipz")) flipz = us::any_cast(parsedArgs["flipz"]); try { itk::TimeProbe clock; clock.Start(); if (!ist::PathExists(out_folder)) { MITK_INFO << "Creating output directory"; ist::MakeDirectory(out_folder); } MITK_INFO << "Loading data"; // Load mask file. Fit is only performed inside the mask MITK_INFO << "Loading mask image"; auto mask = mitk::DiffusionDataIOHelper::load_itk_image(mask_file); // Load masks covering the true positives for evaluation purposes MITK_INFO << "Loading reference peaks and masks"; std::vector< std::string > anchor_mask_files; auto reference_masks = mitk::DiffusionDataIOHelper::load_itk_images(reference_mask_files_folders, &anchor_mask_files); auto reference_peaks = mitk::DiffusionDataIOHelper::load_itk_images(reference_peaks_files_folders); // Load peak image MITK_INFO << "Loading peak image"; auto peak_image = mitk::DiffusionDataIOHelper::load_itk_image(peak_file_name); // Load all candidate tracts MITK_INFO << "Loading candidate tracts"; std::vector< std::string > candidate_tract_files; auto input_candidates = mitk::DiffusionDataIOHelper::load_fibs(candidate_tract_folders, &candidate_tract_files); if (flipx || flipy || flipz) { itk::FlipPeaksFilter< float >::Pointer flipper = itk::FlipPeaksFilter< float >::New(); flipper->SetInput(peak_image); flipper->SetFlipX(flipx); flipper->SetFlipY(flipy); flipper->SetFlipZ(flipz); flipper->Update(); peak_image = flipper->GetOutput(); } mitk::LocaleSwitch localeSwitch("C"); itk::ImageFileWriter< PeakImgType >::Pointer peak_image_writer = itk::ImageFileWriter< PeakImgType >::New(); std::ofstream logfile; logfile.open (out_folder + "scores.txt"); double rmse = 0.0; int iteration = 0; std::string name = "NOANCHOR"; if (parsedArgs.count("a")) { // Load reference tractogram consisting of all known tracts std::string anchors_file = us::any_cast(parsedArgs["a"]); mitk::FiberBundle::Pointer anchor_tractogram = mitk::IOUtil::Load(anchors_file); if ( !(anchor_tractogram.IsNull() || anchor_tractogram->GetNumFibers()==0) ) { // Fit known tracts to peak image to obtain underexplained image MITK_INFO << "Fit anchor tracts"; itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetTractograms({anchor_tractogram}); fitter->SetLambda(static_cast(lambda)); fitter->SetFilterOutliers(filter_outliers); fitter->SetPeakImage(peak_image); fitter->SetVerbose(true); fitter->SetMaskImage(mask); fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); rmse = fitter->GetRMSE(); vnl_vector rms_diff = fitter->GetRmsDiffPerBundle(); name = ist::GetFilenameWithoutExtension(anchors_file); mitk::FiberBundle::Pointer anchor_tracts = fitter->GetTractograms().at(0); anchor_tracts->SetFiberColors(255,255,255); mitk::IOUtil::Save(anchor_tracts, out_folder + boost::lexical_cast(static_cast(100000*rms_diff[0])) + "_" + name + ".fib"); logfile << name << " " << setprecision(5) << rms_diff[0] << "\n"; peak_image = fitter->GetUnderexplainedImage(); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + "Residual_" + name + ".nii.gz"); peak_image_writer->Update(); } } if (use_weights || use_num_streamlines) { MITK_INFO << "Using tract weights as scores"; unsigned int c = 0; for (auto fib : input_candidates) { int mod = 1; float score = 0; if (use_weights) { score = fib->GetFiberWeight(0); mod = 100000; } else if (use_num_streamlines) score = fib->GetNumFibers(); fib->ColorFibersByOrientation(); std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast(static_cast(mod*score)) + "_" + bundle_name + ".fib"); unsigned int num_voxels = 0; { itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New(); masks_filter->SetInputImage(mask); - masks_filter->SetBinaryOutput(true); + masks_filter->SetMode(TDI_MODE::BINARY); masks_filter->SetFiberBundle(fib); masks_filter->SetUseImageGeometry(true); masks_filter->Update(); num_voxels = masks_filter->GetNumCoveredVoxels(); } float weight_sum = 0; for (unsigned int i=0; iGetNumFibers(); i++) weight_sum += fib->GetFiberWeight(i); std::cout.rdbuf (old); // <-- restore logfile << bundle_name << " " << setprecision(5) << score << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; ++c; } } else if (!greedy_add) { MITK_INFO << "Fit candidate tracts"; itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetLambda(static_cast(lambda)); fitter->SetFilterOutliers(filter_outliers); fitter->SetVerbose(true); fitter->SetPeakImage(peak_image); fitter->SetMaskImage(mask); fitter->SetTractograms(input_candidates); fitter->SetFitIndividualFibers(true); if (regu=="MSM") fitter->SetRegularization(VnlCostFunction::REGU::MSM); else if (regu=="Variance") fitter->SetRegularization(VnlCostFunction::REGU::VARIANCE); else if (regu=="Lasso") fitter->SetRegularization(VnlCostFunction::REGU::LASSO); else if (regu=="VoxelVariance") fitter->SetRegularization(VnlCostFunction::REGU::VOXEL_VARIANCE); else if (regu=="GroupLasso") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_LASSO); else if (regu=="GroupVariance") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_VARIANCE); else if (regu=="NONE") fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); vnl_vector rms_diff = fitter->GetRmsDiffPerBundle(); unsigned int c = 0; for (auto fib : input_candidates) { std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect if (filter_zero_weights) fib = fib->FilterByWeights(0); mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast((int)(100000*rms_diff[c])) + "_" + bundle_name + ".fib"); unsigned int num_voxels = 0; { itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New(); masks_filter->SetInputImage(mask); - masks_filter->SetBinaryOutput(true); + masks_filter->SetMode(TDI_MODE::BINARY); masks_filter->SetFiberBundle(fib); masks_filter->SetUseImageGeometry(true); masks_filter->Update(); num_voxels = masks_filter->GetNumCoveredVoxels(); } float weight_sum = 0; for (unsigned int i=0; iGetNumFibers(); i++) weight_sum += fib->GetFiberWeight(i); std::cout.rdbuf (old); // <-- restore logfile << bundle_name << " " << setprecision(5) << rms_diff[c] << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; ++c; } mitk::FiberBundle::Pointer out_fib = mitk::FiberBundle::New(); out_fib = out_fib->AddBundles(input_candidates); out_fib->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(out_fib, out_folder + "AllCandidates.fib"); peak_image = fitter->GetUnderexplainedImage(); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + "Residual_AllCandidates.nii.gz"); peak_image_writer->Update(); } else { MITK_INFO << "RMSE: " << setprecision(5) << rmse; // fitter->SetPeakImage(peak_image); // Iteratively add candidate bundles in a greedy manner while (!input_candidates.empty()) { double next_rmse = rmse; mitk::FiberBundle::Pointer best_candidate = nullptr; PeakImgType::Pointer best_candidate_peak_image = nullptr; for (unsigned int i=0; iSetLambda(static_cast(lambda)); fitter->SetFilterOutliers(filter_outliers); fitter->SetVerbose(false); fitter->SetPeakImage(peak_image); fitter->SetMaskImage(mask); // ****************************** fitter->SetTractograms({input_candidates.at(i)}); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect fitter->Update(); std::cout.rdbuf (old); // <-- restore double candidate_rmse = fitter->GetRMSE(); if (candidate_rmseGetTractograms().at(0); best_candidate_peak_image = fitter->GetUnderexplainedImage(); } } if (best_candidate.IsNull()) break; // fitter->SetPeakImage(peak_image); peak_image = best_candidate_peak_image; unsigned int i=0; std::vector< mitk::FiberBundle::Pointer > remaining_candidates; std::vector< std::string > remaining_candidate_files; for (auto fib : input_candidates) { if (fib!=best_candidate) { remaining_candidates.push_back(fib); remaining_candidate_files.push_back(candidate_tract_files.at(i)); } else name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(i)); ++i; } input_candidates = remaining_candidates; candidate_tract_files = remaining_candidate_files; iteration++; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect // Save winning candidate if (filter_zero_weights) best_candidate = best_candidate->FilterByWeights(0); mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast(iteration) + "_" + name + ".fib"); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + boost::lexical_cast(iteration) + "_" + name + ".nrrd"); peak_image_writer->Update(); std::cout.rdbuf (old); // <-- restore // logfile << name << " " << setprecision(5) << score << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; } } clock.Stop(); int h = static_cast(clock.GetTotal()/3600); int m = (static_cast(clock.GetTotal())%3600)/60; int s = static_cast(clock.GetTotal())%60; MITK_INFO << "Plausibility estimation took " << h << "h, " << m << "m and " << s << "s"; logfile.close(); } 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/TractographyEvaluation/GetOverlappingTracts.cpp b/Modules/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp index 67879d3..855da07 100644 --- a/Modules/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp +++ b/Modules/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp @@ -1,167 +1,167 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; /*! \brief Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Get Overlapping Tracts"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Find tracts that overlap with the reference masks or tracts"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::StringList, "Input:", "input tractograms (.fib/.trk/.tck/.dcm)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output Folder:", "move input tracts that do/don't overlap here", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("reference", "", mitkDiffusionCommandLineParser::StringList, "Reference:", "reference tractograms or mask images", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("overlap_fraction", "", mitkDiffusionCommandLineParser::Float, "Overlap fraction:", "", 0.9); parser.addArgument("use_any_overlap", "", mitkDiffusionCommandLineParser::Bool, "Use any overlap:", "Don't find maximum overlap but use first overlap larger threshold"); parser.addArgument("dont_save_tracts", "", mitkDiffusionCommandLineParser::Bool, "Don't save tracts:", "if true, only text files documenting the overlaps are saved and no tract files are copied"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::StringContainerType input = us::any_cast(parsedArgs["i"]); mitkDiffusionCommandLineParser::StringContainerType reference = us::any_cast(parsedArgs["reference"]); std::string out_folder = us::any_cast(parsedArgs["o"]); bool use_any_overlap = false; if (parsedArgs.count("use_any_overlap")) use_any_overlap = us::any_cast(parsedArgs["use_any_overlap"]); bool dont_save_tracts = false; if (parsedArgs.count("dont_save_tracts")) dont_save_tracts = us::any_cast(parsedArgs["dont_save_tracts"]); float overlap_threshold = 0.9; if (parsedArgs.count("overlap_fraction")) overlap_threshold = us::any_cast(parsedArgs["overlap_fraction"]); try { MITK_INFO << "Loading references"; std::vector< std::string > reference_names; auto masks = mitk::DiffusionDataIOHelper::load_itk_images(reference, &reference_names); auto reference_fibs = mitk::DiffusionDataIOHelper::load_fibs(reference, &reference_names); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect itk::TractDensityImageFilter< ItkFloatImgType >::Pointer filter = itk::TractDensityImageFilter< ItkFloatImgType >::New(); filter->SetUpsamplingFactor(0.25); - filter->SetBinaryOutput(true); + filter->SetMode(TDI_MODE::BINARY); for (auto fib : reference_fibs) { filter->SetFiberBundle(fib); filter->Update(); masks.push_back(filter->GetOutput()); } std::cout.rdbuf (old); // <-- restore MITK_INFO << "Loading input tractograms"; std::vector< std::string > input_names; auto input_fibs = mitk::DiffusionDataIOHelper::load_fibs(input, &input_names); MITK_INFO << "Finding overlaps"; std::ofstream logfile; logfile.open (out_folder + "Overlaps.txt"); std::ofstream logfile2; logfile2.open (out_folder + "AllOverlaps.txt"); boost::progress_display disp(input.size()); unsigned int c = 0; for (auto fib : input_fibs) { ++disp; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect bool is_overlapping = false; float overlap = 0; float max_overlap = 0; std::string max_ref = "-"; int i = 0; std::string overlap_string = ist::GetFilenameWithoutExtension(input_names.at(c)); for (auto m : masks) { overlap = fib->GetOverlap(m); if (overlap>max_overlap) { max_overlap = overlap; max_ref = ist::GetFilenameWithoutExtension(reference_names.at(i)); } if (use_any_overlap && overlap>=overlap_threshold) break; overlap_string += " " + ist::GetFilenameWithoutExtension(reference_names.at(i)) + " " + boost::lexical_cast(overlap); ++i; } if (overlap>=overlap_threshold) is_overlapping = true; logfile << ist::GetFilenameWithoutExtension(input_names.at(c)) << " - " << max_ref << ": " << boost::lexical_cast(max_overlap) << "\n"; logfile2 << overlap_string << "\n"; if (!dont_save_tracts && is_overlapping) ist::CopyAFile(input_names.at(c), out_folder + ist::GetFilenameName(input_names.at(c))); std::cout.rdbuf (old); // <-- restore ++c; } logfile.close(); logfile2.close(); } 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/TractographyEvaluation/MergeOverlappingTracts.cpp b/Modules/DiffusionCmdApps/TractographyEvaluation/MergeOverlappingTracts.cpp index cf440f6..ebe3fa0 100644 --- a/Modules/DiffusionCmdApps/TractographyEvaluation/MergeOverlappingTracts.cpp +++ b/Modules/DiffusionCmdApps/TractographyEvaluation/MergeOverlappingTracts.cpp @@ -1,214 +1,214 @@ /*=================================================================== 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 typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUIntImgType; /*! \brief */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Merge Overlapping Tracts"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::StringList, "Input:", "input tracts", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output Folder:", "output folder", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("overlap", "", mitkDiffusionCommandLineParser::Float, "Overlap threshold:", "Tracts with overlap larger than this threshold are merged", 0.8, false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::StringContainerType input_folder = us::any_cast(parsedArgs["i"]); std::string out_folder = us::any_cast(parsedArgs["o"]); float overlap = 0.8; if (parsedArgs.count("overlap")) overlap = us::any_cast(parsedArgs["overlap"]); try { if (!ist::PathExists(out_folder)) ist::MakeDirectory(out_folder); std::vector< mitk::FiberBundle::Pointer > fibs = mitk::DiffusionDataIOHelper::load_fibs(input_folder); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::FiberBundle::Pointer combined = mitk::FiberBundle::New(); combined = combined->AddBundles(fibs); itk::TractsToFiberEndingsImageFilter< ItkFloatImgType >::Pointer endings = itk::TractsToFiberEndingsImageFilter< ItkFloatImgType >::New(); endings->SetFiberBundle(combined); endings->SetUpsamplingFactor(0.25); endings->Update(); ItkFloatImgType::Pointer ref_image = endings->GetOutput(); std::cout.rdbuf (old); // <-- restore for (int its = 0; its<3; its++) { std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect std::vector< ItkFloatImgType::Pointer > mask_images; for (auto fib : fibs) { itk::TractDensityImageFilter< ItkFloatImgType >::Pointer masks = itk::TractDensityImageFilter< ItkFloatImgType >::New(); masks->SetInputImage(ref_image); - masks->SetBinaryOutput(true); + masks->SetMode(TDI_MODE::BINARY); masks->SetFiberBundle(fib); masks->SetUseImageGeometry(true); masks->Update(); mask_images.push_back(masks->GetOutput()); } int r=0; vnl_matrix< int > mat; mat.set_size(mask_images.size(), mask_images.size()); mat.fill(0); for (auto m1 : mask_images) { float max_overlap = overlap; int c = 0; for (auto m2 : mask_images) { if (c<=r) { ++c; continue; } itk::ImageRegionConstIterator it1(m1, m1->GetLargestPossibleRegion()); itk::ImageRegionConstIterator it2(m2, m2->GetLargestPossibleRegion()); unsigned int c1 = 0; unsigned int c2 = 0; unsigned int intersect = 0; while( !it1.IsAtEnd() ) { if( it1.Get()>0 && it2.Get()>0) ++intersect; if(it1.Get()>0) ++c1; if(it2.Get()>0) ++c2; ++it1; ++it2; } if ( (float)intersect/c1>max_overlap ) { max_overlap = (float)intersect/c1; mat.put(r,c, 1); } if ( (float)intersect/c2>max_overlap ) { max_overlap = (float)intersect/c2; mat.put(r,c, 1); } ++c; } ++r; } std::vector< mitk::FiberBundle::Pointer > out_fibs; std::vector< bool > used; for (unsigned int i=0; i0) { fib = fib->AddBundle(fibs.at(c)); used[c] = true; } } out_fibs.push_back(fib); } std::cout.rdbuf (old); // <-- restore MITK_INFO << fibs.size() << " --> " << out_fibs.size(); if (fibs.size()==out_fibs.size()) break; fibs = out_fibs; } int c = 0; for (auto fib : fibs) { std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::IOUtil::Save(fib, out_folder + "/bundle_" + boost::lexical_cast(c) + ".trk"); std::cout.rdbuf (old); // <-- restore ++c; } } catch (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/DiffusionCore/mitkDiffusionFunctionCollection.cpp b/Modules/DiffusionCore/mitkDiffusionFunctionCollection.cpp index 0308e98..541ef28 100644 --- a/Modules/DiffusionCore/mitkDiffusionFunctionCollection.cpp +++ b/Modules/DiffusionCore/mitkDiffusionFunctionCollection.cpp @@ -1,713 +1,713 @@ /*=================================================================== 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 "mitkDiffusionFunctionCollection.h" #include "mitkNumericTypes.h" #include #include #include #include #include #include "itkVectorContainer.h" #include "vnl/vnl_vector.h" #include #include #include #include #include // Intersect a finite line (with end points p0 and p1) with all of the // cells of a vtkImageData std::vector< std::pair< itk::Index<3>, double > > mitk::imv::IntersectImage(const itk::Vector& spacing, itk::Index<3>& si, itk::Index<3>& ei, itk::ContinuousIndex& sf, itk::ContinuousIndex& ef) { std::vector< std::pair< itk::Index<3>, double > > out; if (si == ei) { double d[3]; for (int i=0; i<3; ++i) d[i] = static_cast(sf[i]-ef[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); out.push_back( std::pair< itk::Index<3>, double >(si, len) ); return out; } double bounds[6]; double entrancePoint[3]; double exitPoint[3]; double startPoint[3]; double endPoint[3]; double t0, t1; for (unsigned int i=0; i<3; ++i) { startPoint[i] = static_cast(sf[i]); endPoint[i] = static_cast(ef[i]); if (si[i]>ei[i]) { auto t = si[i]; si[i] = ei[i]; ei[i] = t; } } for (auto x = si[0]; x<=ei[0]; ++x) for (auto y = si[1]; y<=ei[1]; ++y) for (auto z = si[2]; z<=ei[2]; ++z) { bounds[0] = static_cast(x) - 0.5; bounds[1] = static_cast(x) + 0.5; bounds[2] = static_cast(y) - 0.5; bounds[3] = static_cast(y) + 0.5; bounds[4] = static_cast(z) - 0.5; bounds[5] = static_cast(z) + 0.5; int entryPlane; int exitPlane; int hit = vtkBox::IntersectWithLine(bounds, startPoint, endPoint, t0, t1, entrancePoint, exitPoint, entryPlane, exitPlane); if (hit) { if (entryPlane>=0 && exitPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (exitPoint[i] - entrancePoint[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } else if (entryPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (static_cast(ef[i]) - entrancePoint[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } else if (exitPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (exitPoint[i]-static_cast(sf[i]))*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } } } return out; } //------------------------- SH-function ------------------------------------ double mitk::sh::factorial(int number) { if(number <= 1) return 1; double result = 1.0; for(int i=1; i<=number; i++) result *= i; return result; } void mitk::sh::Cart2Sph(double x, double y, double z, double *spherical) { double phi, th, rad; rad = sqrt(x*x+y*y+z*z); if( rad < mitk::eps ) { th = itk::Math::pi/2; phi = itk::Math::pi/2; } else { th = acos(z/rad); phi = atan2(y, x); } spherical[0] = phi; spherical[1] = th; spherical[2] = rad; } vnl_vector_fixed mitk::sh::Sph2Cart(const double& theta, const double& phi, const double& rad) { vnl_vector_fixed dir; dir[0] = rad * sin(theta) * cos(phi); dir[1] = rad * sin(theta) * sin(phi); dir[2] = rad * cos(theta); return dir; } double mitk::sh::legendre0(int l) { if( l%2 != 0 ) { return 0; } else { double prod1 = 1.0; for(int i=1;i(::boost::math::spherical_harmonic_r(static_cast(k), -m, theta, phi)); else if (m==0) return static_cast(::boost::math::spherical_harmonic_r(static_cast(k), m, theta, phi)); else return pow(-1.0,m)*sqrt(2.0)*static_cast(::boost::math::spherical_harmonic_i(static_cast(k), m, theta, phi)); } else { double plm = static_cast(::boost::math::legendre_p(k,abs(m),-cos(theta))); double mag = sqrt((2.0*k+1.0)/(4.0*itk::Math::pi)*::boost::math::factorial(k-abs(m))/::boost::math::factorial(k+abs(m)))*plm; if (m>0) return mag*static_cast(cos(m*phi)); else if (m==0) return mag; else return mag*static_cast(sin(-m*phi)); } return 0; } mitk::OdfImage::ItkOdfImageType::Pointer mitk::convert::GetItkOdfFromShImage(mitk::Image::Pointer mitkImage) { mitk::ShImage::Pointer mitkShImage = dynamic_cast(mitkImage.GetPointer()); if (mitkShImage.IsNull()) mitkThrow() << "Input image is not a SH image!"; mitk::OdfImage::ItkOdfImageType::Pointer output; switch (mitkShImage->ShOrder()) { case 2: { typedef itk::ShToOdfImageFilter< float, 2 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } case 4: { typedef itk::ShToOdfImageFilter< float, 4 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } case 6: { typedef itk::ShToOdfImageFilter< float, 6 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } case 8: { typedef itk::ShToOdfImageFilter< float, 8 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } case 10: { typedef itk::ShToOdfImageFilter< float, 10 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } case 12: { typedef itk::ShToOdfImageFilter< float, 12 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } default: mitkThrow() << "SH orders higher than 12 are not supported!"; } return output; } mitk::OdfImage::Pointer mitk::convert::GetOdfFromShImage(mitk::Image::Pointer mitkImage) { mitk::OdfImage::Pointer image = mitk::OdfImage::New(); auto img = GetItkOdfFromShImage(mitkImage); image->InitializeByItk( img.GetPointer() ); image->SetVolume( img->GetBufferPointer() ); return image; } mitk::OdfImage::ItkOdfImageType::Pointer mitk::convert::GetItkOdfFromTensorImage(mitk::Image::Pointer mitkImage) { typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( GetItkTensorFromTensorImage(mitkImage) ); filter->Update(); return filter->GetOutput(); } mitk::TensorImage::ItkTensorImageType::Pointer mitk::convert::GetItkTensorFromTensorImage(mitk::Image::Pointer mitkImage) { typedef mitk::ImageToItk< mitk::TensorImage::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImage); caster->Update(); return caster->GetOutput(); } mitk::PeakImage::ItkPeakImageType::Pointer mitk::convert::GetItkPeakFromPeakImage(mitk::Image::Pointer mitkImage) { typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImage); caster->SetCopyMemFlag(true); caster->Update(); return caster->GetOutput(); } mitk::OdfImage::Pointer mitk::convert::GetOdfFromTensorImage(mitk::Image::Pointer mitkImage) { mitk::OdfImage::Pointer image = mitk::OdfImage::New(); auto img = GetItkOdfFromTensorImage(mitkImage); image->InitializeByItk( img.GetPointer() ); image->SetVolume( img->GetBufferPointer() ); return image; } mitk::OdfImage::ItkOdfImageType::Pointer mitk::convert::GetItkOdfFromOdfImage(mitk::Image::Pointer mitkImage) { typedef mitk::ImageToItk< mitk::OdfImage::ItkOdfImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImage); caster->Update(); return caster->GetOutput(); } vnl_matrix mitk::sh::CalcShBasisForDirections(unsigned int sh_order, vnl_matrix U, bool mrtrix) { vnl_matrix sh_basis = vnl_matrix(U.cols(), (sh_order*sh_order + sh_order + 2)/2 + sh_order ); for(unsigned int i=0; i(sh_order); k+=2) { for(int m=-k; m<=k; m++) { int j = (k*k + k + 2)/2 + m - 1; double phi = U(0,i); double th = U(1,i); sh_basis(i,j) = mitk::sh::Yj(m,k,th,phi, mrtrix); } } } return sh_basis; } unsigned int mitk::sh::ShOrder(int num_coeffs) { int c=3, d=2-2*num_coeffs; int D = c*c-4*d; if (D>0) { int s = (-c+static_cast(sqrt(D)))/2; if (s<0) s = (-c-static_cast(sqrt(D)))/2; return static_cast(s); } else if (D==0) return static_cast(-c/2); return 0; } float mitk::sh::GetValue(const vnl_vector &coefficients, const int &sh_order, const double theta, const double phi, const bool mrtrix) { float val = 0; for(int k=0; k<=sh_order; k+=2) { for(int m=-k; m<=k; m++) { unsigned int j = static_cast((k*k + k + 2)/2 + m - 1); val += coefficients[j] * mitk::sh::Yj(m, k, theta, phi, mrtrix); } } return val; } float mitk::sh::GetValue(const vnl_vector &coefficients, const int &sh_order, const vnl_vector_fixed &dir, const bool mrtrix) { double spherical[3]; mitk::sh::Cart2Sph(dir[0], dir[1], dir[2], spherical); float val = 0; for(int k=0; k<=sh_order; k+=2) { for(int m=-k; m<=k; m++) { int j = (k*k + k + 2)/2 + m - 1; val += coefficients[j] * mitk::sh::Yj(m, k, spherical[1], spherical[0], mrtrix); } } return val; } //------------------------- gradients-function ------------------------------------ mitk::gradients::GradientDirectionContainerType::ConstPointer mitk::gradients::ReadBvalsBvecs(std::string bvals_file, std::string bvecs_file, double& reference_bval) { mitk::gradients::GradientDirectionContainerType::Pointer directioncontainer = mitk::gradients::GradientDirectionContainerType::New(); std::vector bvec_entries; if (!itksys::SystemTools::FileExists(bvecs_file)) mitkThrow() << "bvecs file not existing: " << bvecs_file; else { std::string line; std::ifstream myfile (bvecs_file.c_str()); if (myfile.is_open()) { while (std::getline(myfile, line)) { std::vector strs; - boost::split(strs,line,boost::is_any_of("\t \n")); + boost::split(strs,line,boost::is_any_of("\t \n\r")); for (auto token : strs) { if (!token.empty()) { try { bvec_entries.push_back(boost::lexical_cast(token)); } catch(...) { mitkThrow() << "Encountered invalid bvecs file entry >" << token << "<"; } } } } myfile.close(); } else { mitkThrow() << "bvecs file could not be opened: " << bvals_file; } } reference_bval = -1; std::vector bval_entries; if (!itksys::SystemTools::FileExists(bvals_file)) mitkThrow() << "bvals file not existing: " << bvals_file; else { std::string line; std::ifstream myfile (bvals_file.c_str()); if (myfile.is_open()) { while (std::getline(myfile, line)) { std::vector strs; - boost::split(strs,line,boost::is_any_of("\t \n")); + boost::split(strs,line,boost::is_any_of("\t \n\r")); for (auto token : strs) { if (!token.empty()) { try { bval_entries.push_back(boost::lexical_cast(token)); if (bval_entries.back()>reference_bval) reference_bval = bval_entries.back(); } catch(...) { mitkThrow() << "Encountered invalid bvals file entry >" << token << "<"; } } } } myfile.close(); } else { mitkThrow() << "bvals file could not be opened: " << bvals_file; } } for(unsigned int i=0; i0) { double factor = b_val/reference_bval; if(vec.magnitude() > 0) { vec.normalize(); vec[0] = sqrt(factor)*vec[0]; vec[1] = sqrt(factor)*vec[1]; vec[2] = sqrt(factor)*vec[2]; } } directioncontainer->InsertElement(i,vec); } return GradientDirectionContainerType::ConstPointer(directioncontainer); } void mitk::gradients::WriteBvalsBvecs(std::string bvals_file, std::string bvecs_file, GradientDirectionContainerType::ConstPointer gradients, double reference_bval) { std::ofstream myfile; myfile.open (bvals_file.c_str()); for(unsigned int i=0; iSize(); i++) { double twonorm = gradients->ElementAt(i).two_norm(); myfile << std::round(reference_bval*twonorm*twonorm) << " "; } myfile.close(); std::ofstream myfile2; myfile2.open (bvecs_file.c_str()); for(int j=0; j<3; j++) { for(unsigned int i=0; iSize(); i++) { GradientDirectionType direction = gradients->ElementAt(i); direction.normalize(); myfile2 << direction.get(j) << " "; } myfile2 << std::endl; } } std::vector mitk::gradients::GetAllUniqueDirections(const BValueMap & refBValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer ) { IndiciesVector directioncontainer; auto mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); auto containerIt = directioncontainer.begin(); bool directionExist = false; while(containerIt != directioncontainer.end()) { if (fabs(dot_product(refGradientsContainer->ElementAt(*containerIt), refGradientsContainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { directioncontainer.push_back(wntIndex); } } } return directioncontainer; } bool mitk::gradients::CheckForDifferingShellDirections(const BValueMap & refBValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer) { auto mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ auto mapIterator_2 = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator_2++; //skip bzero Values for( ; mapIterator_2 != refBValueMap.end(); mapIterator_2++){ if(mapIterator_2 == mapIterator) continue; IndiciesVector currentShell = mapIterator->second; IndiciesVector testShell = mapIterator_2->second; for (unsigned int i = 0; i< currentShell.size(); i++) if (fabs(dot_product(refGradientsContainer->ElementAt(currentShell[i]), refGradientsContainer->ElementAt(testShell[i]))) <= 0.9998) { return true; } } } return false; } vnl_matrix mitk::gradients::ComputeSphericalFromCartesian(const IndiciesVector & refShell, const GradientDirectionContainerType * refGradientsContainer) { vnl_matrix Q(3, refShell.size()); Q.fill(0.0); for(unsigned int i = 0; i < refShell.size(); i++) { GradientDirectionType dir = refGradientsContainer->ElementAt(refShell[i]); double x = dir.normalize().get(0); double y = dir.normalize().get(1); double z = dir.normalize().get(2); double cart[3]; mitk::sh::Cart2Sph(x,y,z,cart); Q(0,i) = cart[0]; Q(1,i) = cart[1]; Q(2,i) = cart[2]; } return Q; } vnl_matrix mitk::gradients::ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder) { vnl_matrix SHBasisOutput(QBallReference.cols(), (LOrder+1)*(LOrder+2)*0.5); SHBasisOutput.fill(0.0); for(unsigned int i=0; i 1){ mapIterator++; //skip bzero Values vnl_vector_fixed vec; vec.fill(0.0); directioncontainer->push_back(vec); } for( ; mapIterator != bValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); mitk::gradients::GradientDirectionContainerType::Iterator containerIt = directioncontainer->Begin(); bool directionExist = false; while(containerIt != directioncontainer->End()) { if (fabs(dot_product(containerIt.Value(), origninalGradentcontainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { GradientDirectionType dir(origninalGradentcontainer->ElementAt(wntIndex)); directioncontainer->push_back(dir.normalize()); } } } return directioncontainer; } diff --git a/Modules/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp b/Modules/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp index f26dd7f..d4b99d1 100644 --- a/Modules/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp +++ b/Modules/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp @@ -1,910 +1,910 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _TrackingForestHandler_cpp #define _TrackingForestHandler_cpp #include "mitkTrackingHandlerRandomForest.h" #include #include namespace mitk { template< int ShOrder, int NumberOfSignalFeatures > TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::TrackingHandlerRandomForest() : m_WmSampleDistance(-1) , m_NumTrees(30) , m_MaxTreeDepth(25) , m_SampleFraction(1.0) , m_NumberOfSamples(0) , m_GmSamplesPerVoxel(-1) , m_BidirectionalFiberSampling(false) , m_ZeroDirWmFeatures(true) , m_MaxNumWmSamples(-1) { vnl_vector_fixed ref; ref.fill(0); ref[0]=1; itk::OrientationDistributionFunction< float, 200 > odf; m_DirectionContainer.clear(); for (unsigned int i = 0; i odf_dir; odf_dir[0] = odf.GetDirection(i)[0]; odf_dir[1] = odf.GetDirection(i)[1]; odf_dir[2] = odf.GetDirection(i)[2]; if (dot_product(ref, odf_dir)>0) // only used directions on one hemisphere m_DirectionContainer.push_back(odf_dir); // store indices for later mapping the classifier output to the actual direction } m_OdfFloatDirs.set_size(m_DirectionContainer.size(), 3); for (unsigned int i=0; i TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::~TrackingHandlerRandomForest() { } template< int ShOrder, int NumberOfSignalFeatures > bool TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::WorldToIndex(itk::Point& pos, itk::Index<3>& index) { m_DwiFeatureImages.at(0)->TransformPhysicalPointToIndex(pos, index); return m_DwiFeatureImages.at(0)->GetLargestPossibleRegion().IsInside(index); } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InputDataValidForTracking() { if (m_InputDwis.empty()) mitkThrow() << "No diffusion-weighted images set!"; if (!IsForestValid()) mitkThrow() << "No or invalid random forest detected!"; } template< int ShOrder, int NumberOfSignalFeatures> template typename std::enable_if< NumberOfSignalFeatures <=99, T >::type TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InitDwiImageFeatures(mitk::Image::Pointer mitk_dwi) { MITK_INFO << "Calculating spherical harmonics features"; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter InterpolationFilterType; typename InterpolationFilterType::Pointer filter = InterpolationFilterType::New(); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(mitk_dwi)); filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(mitk_dwi), mitk::DiffusionPropertyHelper::GetItkVectorImage(mitk_dwi) ); filter->SetLambda(0.006); filter->SetNormalizationMethod(InterpolationFilterType::QBAR_RAW_SIGNAL); filter->Update(); m_DwiFeatureImages.push_back(filter->GetCoefficientImage()); return true; } template< int ShOrder, int NumberOfSignalFeatures> template typename std::enable_if< NumberOfSignalFeatures >=100, T >::type TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InitDwiImageFeatures(mitk::Image::Pointer mitk_dwi) { MITK_INFO << "Interpolating raw dwi signal features"; typedef itk::AnalyticalDiffusionQballReconstructionImageFilter InterpolationFilterType; typename InterpolationFilterType::Pointer filter = InterpolationFilterType::New(); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(mitk_dwi)); filter->SetGradientImage( mitk::DiffusionPropertyHelper::GetOriginalGradientContainer(mitk_dwi), mitk::DiffusionPropertyHelper::GetItkVectorImage(mitk_dwi) ); filter->SetLambda(0.006); filter->SetNormalizationMethod(InterpolationFilterType::QBAR_RAW_SIGNAL); filter->Update(); typename DwiFeatureImageType::Pointer dwiFeatureImage = DwiFeatureImageType::New(); dwiFeatureImage->SetSpacing(filter->GetOutput()->GetSpacing()); dwiFeatureImage->SetOrigin(filter->GetOutput()->GetOrigin()); dwiFeatureImage->SetDirection(filter->GetOutput()->GetDirection()); dwiFeatureImage->SetLargestPossibleRegion(filter->GetOutput()->GetLargestPossibleRegion()); dwiFeatureImage->SetBufferedRegion(filter->GetOutput()->GetLargestPossibleRegion()); dwiFeatureImage->SetRequestedRegion(filter->GetOutput()->GetLargestPossibleRegion()); dwiFeatureImage->Allocate(); // get signal values and store them in the feature image vnl_vector_fixed ref; ref.fill(0); ref[0]=1; itk::OrientationDistributionFunction< float, 2*NumberOfSignalFeatures > odf; itk::ImageRegionIterator< typename InterpolationFilterType::OutputImageType > it(filter->GetOutput(), filter->GetOutput()->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { typename DwiFeatureImageType::PixelType pix; int f = 0; for (unsigned int i = 0; i0) // only used directions on one hemisphere { pix[f] = it.Get()[i]; f++; } } dwiFeatureImage->SetPixel(it.GetIndex(), pix); ++it; } m_DwiFeatureImages.push_back(dwiFeatureImage); return true; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InitForTracking() { MITK_INFO << "Initializing random forest tracker."; if (m_NeedsDataInit) { InputDataValidForTracking(); m_DwiFeatureImages.clear(); InitDwiImageFeatures<>(m_InputDwis.at(0)); // initialize interpolators m_DwiFeatureImageInterpolator = DwiFeatureImageInterpolatorType::New(); m_DwiFeatureImageInterpolator->SetInputImage(m_DwiFeatureImages.at(0)); m_AdditionalFeatureImageInterpolators.clear(); for (auto afi_vec : m_AdditionalFeatureImages) { std::vector< FloatImageInterpolatorType::Pointer > v; for (auto img : afi_vec) { FloatImageInterpolatorType::Pointer interp = FloatImageInterpolatorType::New(); interp->SetInputImage(img); v.push_back(interp); } m_AdditionalFeatureImageInterpolators.push_back(v); } auto double_dir = m_DwiFeatureImages.at(0)->GetDirection().GetVnlMatrix(); for (int r=0; r<3; r++) for (int c=0; c<3; c++) { m_FloatImageRotation[r][c] = double_dir[r][c]; } m_NeedsDataInit = false; } this->CalculateMinVoxelSize(); } template< int ShOrder, int NumberOfSignalFeatures > vnl_vector_fixed TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { vnl_vector_fixed output_direction; output_direction.fill(0); itk::Index<3> idx; m_DwiFeatureImages.at(0)->TransformPhysicalPointToIndex(pos, idx); bool check_last_dir = false; vnl_vector_fixed last_dir; if (!olddirs.empty()) { last_dir = olddirs.back(); if (last_dir.magnitude()>0.5) check_last_dir = true; } if (!m_Parameters->m_InterpolateTractographyData && oldIndex==idx) return last_dir; // store feature pixel values in a vigra data type vigra::MultiArray<2, float> featureData = vigra::MultiArray<2, float>( vigra::Shape2(1,m_Forest->GetNumFeatures()) ); featureData.init(0.0); typename DwiFeatureImageType::PixelType dwiFeaturePixel = mitk::imv::GetImageValue< typename DwiFeatureImageType::PixelType >(pos, m_Parameters->m_InterpolateTractographyData, m_DwiFeatureImageInterpolator); for (unsigned int f=0; f direction_matrix = m_DwiFeatureImages.at(0)->GetDirection().GetVnlMatrix(); vnl_matrix_fixed inverse_direction_matrix = m_DwiFeatureImages.at(0)->GetInverseDirection().GetVnlMatrix(); // append normalized previous direction(s) to feature vector int i = 0; vnl_vector_fixed ref; ref.fill(0); ref[0]=1; for (auto d : olddirs) { vnl_vector_fixed tempD; tempD[0] = d[0]; tempD[1] = d[1]; tempD[2] = d[2]; if (m_Parameters->m_FlipX) tempD[0] *= -1; if (m_Parameters->m_FlipY) tempD[1] *= -1; if (m_Parameters->m_FlipZ) tempD[2] *= -1; tempD = inverse_direction_matrix * tempD; last_dir[0] = tempD[0]; last_dir[1] = tempD[1]; last_dir[2] = tempD[2]; int c = 0; for (int f=NumberOfSignalFeatures+3*i; f0) { int c = 0; for (auto interpolator : m_AdditionalFeatureImageInterpolators.at(0)) { float v = mitk::imv::GetImageValue(pos, false, interpolator); featureData(0,NumberOfSignalFeatures+m_Parameters->m_NumPreviousDirections*3+c) = v; c++; } } // perform classification vigra::MultiArray<2, float> probs(vigra::Shape2(1, m_Forest->GetNumClasses())); m_Forest->PredictProbabilities(featureData, probs); vnl_vector< float > angles = m_OdfFloatDirs*last_dir; vnl_vector< float > probs2; probs2.set_size(m_DirectionContainer.size()); probs2.fill(0.0); // used for probabilistic direction sampling float probs_sum = 0; float pNonFib = 0; // probability that we left the white matter float w = 0; // weight of the predicted direction for (int i=0; iGetNumClasses(); i++) // for each class (number of possible directions + out-of-wm class) { if (probs(0,i)>0) // if probability of respective class is 0, do nothing { // get label of class (does not correspond to the loop variable i) unsigned int classLabel = m_Forest->IndexToClassLabel(i); if (classLabelm_Mode==MODE::PROBABILISTIC) { probs2[classLabel] = probs(0,i); if (check_last_dir) probs2[classLabel] *= abs_angle; probs_sum += probs2[classLabel]; } else if (m_Parameters->m_Mode==MODE::DETERMINISTIC) { vnl_vector_fixed d = m_DirectionContainer.at(classLabel); // get direction vector assiciated with the respective direction index if (check_last_dir) // do we have a previous streamline direction or did we just start? { if (abs_angle>=m_Parameters->GetAngularThresholdDot()) // is angle between the directions smaller than our hard threshold? { if (angle<0) // make sure we don't walk backwards d *= -1; float w_i = probs(0,i)*abs_angle; output_direction += w_i*d; // weight contribution to output direction with its probability and the angular deviation from the previous direction w += w_i; // increase output weight of the final direction } } else { output_direction += probs(0,i)*d; w += probs(0,i); } } } else pNonFib += probs(0,i); // probability that we are not in the white matter anymore } } if (m_Parameters->m_Mode==MODE::PROBABILISTIC && pNonFib<0.5) { boost::random::discrete_distribution dist(probs2.begin(), probs2.end()); int sampled_idx = 0; for (int i=0; i<50; i++) // we allow 50 trials to exceed m_AngularThreshold { #pragma omp critical { boost::random::variate_generator> sampler(m_Rng, dist); sampled_idx = sampler(); } if ( probs2[sampled_idx]>0.1 && (!check_last_dir || (check_last_dir && fabs(angles[sampled_idx])>=m_Parameters->GetAngularThresholdDot())) ) break; } output_direction = m_DirectionContainer.at(sampled_idx); w = probs2[sampled_idx]; if (check_last_dir && angles[sampled_idx]<0) // make sure we don't walk backwards output_direction *= -1; } // if we did not find a suitable direction, make sure that we return (0,0,0) if (pNonFib>w && w>0) output_direction.fill(0.0); else { vnl_vector_fixed tempD; tempD[0] = output_direction[0]; tempD[1] = output_direction[1]; tempD[2] = output_direction[2]; tempD = direction_matrix * tempD; output_direction[0] = tempD[0]; output_direction[1] = tempD[1]; output_direction[2] = tempD[2]; if (m_Parameters->m_FlipX) output_direction[0] *= -1; if (m_Parameters->m_FlipY) output_direction[1] *= -1; if (m_Parameters->m_FlipZ) output_direction[2] *= -1; if (m_Parameters->m_ApplyDirectionMatrix) output_direction = m_FloatImageRotation*output_direction; } return output_direction * w; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::StartTraining() { m_StartTime = std::chrono::system_clock::now(); InputDataValidForTraining(); InitForTraining(); CalculateTrainingSamples(); MITK_INFO << "Maximum tree depths: " << m_MaxTreeDepth; MITK_INFO << "Sample fraction per tree: " << m_SampleFraction; MITK_INFO << "Number of trees: " << m_NumTrees; DefaultSplitType splitter; splitter.UsePointBasedWeights(true); splitter.SetWeights(m_Weights); splitter.UseRandomSplit(false); splitter.SetPrecision(mitk::eps); splitter.SetMaximumTreeDepth(m_MaxTreeDepth); std::vector< std::shared_ptr< vigra::RandomForest > > trees; int count = 0; #pragma omp parallel for for (int i = 0; i < m_NumTrees; ++i) { std::shared_ptr< vigra::RandomForest > lrf = std::make_shared< vigra::RandomForest >(); lrf->set_options().use_stratification(vigra::RF_NONE); // How the data should be made equal lrf->set_options().sample_with_replacement(true); // if sampled with replacement or not lrf->set_options().samples_per_tree(m_SampleFraction); // Fraction of samples that are used to train a tree lrf->set_options().tree_count(1); // Number of trees that are calculated; lrf->set_options().min_split_node_size(5); // Minimum number of datapoints that must be in a node lrf->ext_param_.max_tree_depth = m_MaxTreeDepth; lrf->learn(m_FeatureData, m_LabelData,vigra::rf::visitors::VisitorBase(),splitter); #pragma omp critical { count++; MITK_INFO << "Tree " << count << " finished training."; trees.push_back(lrf); } } for (int i = 1; i < m_NumTrees; ++i) trees.at(0)->trees_.push_back(trees.at(i)->trees_[0]); std::shared_ptr< vigra::RandomForest > forest = trees.at(0); forest->options_.tree_count_ = m_NumTrees; m_Forest = mitk::TractographyForest::New(forest); MITK_INFO << "Training finsihed"; m_EndTime = std::chrono::system_clock::now(); std::chrono::hours hh = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::minutes mm = std::chrono::duration_cast(m_EndTime - m_StartTime); mm %= 60; MITK_INFO << "Training took " << hh.count() << "h and " << mm.count() << "m"; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InputDataValidForTraining() { if (m_InputDwis.empty()) mitkThrow() << "No diffusion-weighted images set!"; if (m_Tractograms.empty()) mitkThrow() << "No tractograms set!"; if (m_InputDwis.size()!=m_Tractograms.size()) mitkThrow() << "Unequal number of diffusion-weighted images and tractograms detected!"; } template< int ShOrder, int NumberOfSignalFeatures > bool TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::IsForestValid() { int additional_features = 0; if (m_AdditionalFeatureImages.size()>0) additional_features = m_AdditionalFeatureImages.at(0).size(); if (!m_Forest) MITK_INFO << "No forest available!"; else { if (m_Forest->GetNumTrees() <= 0) MITK_ERROR << "Forest contains no trees!"; if ( m_Forest->GetNumFeatures() != static_cast(NumberOfSignalFeatures+3*m_Parameters->m_NumPreviousDirections+additional_features) ) MITK_ERROR << "Wrong number of features in forest: got " << m_Forest->GetNumFeatures() << ", expected " << (NumberOfSignalFeatures+3*m_Parameters->m_NumPreviousDirections+additional_features); } if(m_Forest && m_Forest->GetNumTrees()>0 && m_Forest->GetNumFeatures() == static_cast(NumberOfSignalFeatures+3*m_Parameters->m_NumPreviousDirections+additional_features)) return true; return false; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::InitForTraining() { MITK_INFO << "Spherical signal interpolation and sampling ..."; for (unsigned int i=0; i(m_InputDwis.at(i)); if (i>=m_AdditionalFeatureImages.size()) { m_AdditionalFeatureImages.push_back(std::vector< ItkFloatImgType::Pointer >()); } if (i>=m_FiberVolumeModImages.size()) { ItkFloatImgType::Pointer img = ItkFloatImgType::New(); img->SetSpacing( m_DwiFeatureImages.at(i)->GetSpacing() ); img->SetOrigin( m_DwiFeatureImages.at(i)->GetOrigin() ); img->SetDirection( m_DwiFeatureImages.at(i)->GetDirection() ); img->SetLargestPossibleRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); img->SetBufferedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); img->SetRequestedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); img->Allocate(); img->FillBuffer(1); m_FiberVolumeModImages.push_back(img); } if (m_FiberVolumeModImages.at(i)==nullptr) { m_FiberVolumeModImages.at(i) = ItkFloatImgType::New(); m_FiberVolumeModImages.at(i)->SetSpacing( m_DwiFeatureImages.at(i)->GetSpacing() ); m_FiberVolumeModImages.at(i)->SetOrigin( m_DwiFeatureImages.at(i)->GetOrigin() ); m_FiberVolumeModImages.at(i)->SetDirection( m_DwiFeatureImages.at(i)->GetDirection() ); m_FiberVolumeModImages.at(i)->SetLargestPossibleRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_FiberVolumeModImages.at(i)->SetBufferedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_FiberVolumeModImages.at(i)->SetRequestedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_FiberVolumeModImages.at(i)->Allocate(); m_FiberVolumeModImages.at(i)->FillBuffer(1); } if (i>=m_MaskImages.size()) { ItkUcharImgType::Pointer newMask = ItkUcharImgType::New(); newMask->SetSpacing( m_DwiFeatureImages.at(i)->GetSpacing() ); newMask->SetOrigin( m_DwiFeatureImages.at(i)->GetOrigin() ); newMask->SetDirection( m_DwiFeatureImages.at(i)->GetDirection() ); newMask->SetLargestPossibleRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); newMask->SetBufferedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); newMask->SetRequestedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); newMask->Allocate(); newMask->FillBuffer(1); m_MaskImages.push_back(newMask); } if (m_MaskImages.at(i)==nullptr) { m_MaskImages.at(i) = ItkUcharImgType::New(); m_MaskImages.at(i)->SetSpacing( m_DwiFeatureImages.at(i)->GetSpacing() ); m_MaskImages.at(i)->SetOrigin( m_DwiFeatureImages.at(i)->GetOrigin() ); m_MaskImages.at(i)->SetDirection( m_DwiFeatureImages.at(i)->GetDirection() ); m_MaskImages.at(i)->SetLargestPossibleRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_MaskImages.at(i)->SetBufferedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_MaskImages.at(i)->SetRequestedRegion( m_DwiFeatureImages.at(i)->GetLargestPossibleRegion() ); m_MaskImages.at(i)->Allocate(); m_MaskImages.at(i)->FillBuffer(1); } } // initialize interpolators m_AdditionalFeatureImageInterpolators.clear(); for (auto afi_vec : m_AdditionalFeatureImages) { std::vector< FloatImageInterpolatorType::Pointer > v; for (auto img : afi_vec) { FloatImageInterpolatorType::Pointer interp = FloatImageInterpolatorType::New(); interp->SetInputImage(img); v.push_back(interp); } m_AdditionalFeatureImageInterpolators.push_back(v); } MITK_INFO << "Resampling fibers and calculating number of samples ..."; m_NumberOfSamples = 0; m_SampleUsage.clear(); for (unsigned int t=0; t::Pointer env = itk::TractDensityImageFilter< ItkUcharImgType >::New(); env->SetFiberBundle(m_Tractograms.at(t)); env->SetInputImage(mask); - env->SetBinaryOutput(true); + env->SetMode(TDI_MODE::BINARY); env->SetUseImageGeometry(true); env->Update(); wmmask = env->GetOutput(); if (t>=m_WhiteMatterImages.size()) m_WhiteMatterImages.push_back(wmmask); else m_WhiteMatterImages.at(t) = wmmask; } // Calculate white-matter samples if (m_WmSampleDistance<0) { typename DwiFeatureImageType::Pointer image = m_DwiFeatureImages.at(t); float minSpacing = 1; if(image->GetSpacing()[0]GetSpacing()[1] && image->GetSpacing()[0]GetSpacing()[2]) minSpacing = image->GetSpacing()[0]; else if (image->GetSpacing()[1] < image->GetSpacing()[2]) minSpacing = image->GetSpacing()[1]; else minSpacing = image->GetSpacing()[2]; m_WmSampleDistance = minSpacing*0.5; } m_Tractograms.at(t)->ResampleLinear(m_WmSampleDistance); int wmSamples = m_Tractograms.at(t)->GetNumberOfPoints()-2*m_Tractograms.at(t)->GetNumFibers(); if (m_BidirectionalFiberSampling) wmSamples *= 2; if (m_ZeroDirWmFeatures) wmSamples *= (m_Parameters->m_NumPreviousDirections+1); MITK_INFO << "White matter samples available: " << wmSamples; // upper limit for samples if (m_MaxNumWmSamples>0 && wmSamples>m_MaxNumWmSamples) { if ((float)m_MaxNumWmSamples/wmSamples > 0.8) { m_SampleUsage.push_back(std::vector(wmSamples, true)); m_NumberOfSamples += wmSamples; } else { m_SampleUsage.push_back(std::vector(wmSamples, false)); m_NumberOfSamples += m_MaxNumWmSamples; wmSamples = m_MaxNumWmSamples; MITK_INFO << "Limiting white matter samples to: " << m_MaxNumWmSamples; itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer randgen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); randgen->SetSeed(); int c = 0; while (cGetIntegerVariate(m_MaxNumWmSamples-1); if (m_SampleUsage[t][idx]==false) { m_SampleUsage[t][idx]=true; ++c; } } } } else { m_SampleUsage.push_back(std::vector(wmSamples, true)); m_NumberOfSamples += wmSamples; } // calculate gray-matter samples itk::ImageRegionConstIterator it(wmmask, wmmask->GetLargestPossibleRegion()); int OUTOFWM = 0; while(!it.IsAtEnd()) { if (it.Get()==0 && mask->GetPixel(it.GetIndex())>0) OUTOFWM++; ++it; } MITK_INFO << "Non-white matter voxels: " << OUTOFWM; if (m_GmSamplesPerVoxel>0) { m_GmSamples.push_back(m_GmSamplesPerVoxel); m_NumberOfSamples += m_GmSamplesPerVoxel*OUTOFWM; } else if (OUTOFWM>0) { int gm_per_voxel = 0.5+(float)wmSamples/(float)OUTOFWM; if (gm_per_voxel<=0) gm_per_voxel = 1; m_GmSamples.push_back(gm_per_voxel); m_NumberOfSamples += m_GmSamples.back()*OUTOFWM; MITK_INFO << "Non-white matter samples per voxel: " << m_GmSamples.back(); } else { m_GmSamples.push_back(0); } MITK_INFO << "Non-white matter samples: " << m_GmSamples.back()*OUTOFWM; } MITK_INFO << "Number of samples: " << m_NumberOfSamples; } template< int ShOrder, int NumberOfSignalFeatures > void TrackingHandlerRandomForest< ShOrder, NumberOfSignalFeatures >::CalculateTrainingSamples() { vnl_vector_fixed ref; ref.fill(0); ref[0]=1; m_FeatureData.reshape( vigra::Shape2(m_NumberOfSamples, NumberOfSignalFeatures+m_Parameters->m_NumPreviousDirections*3+m_AdditionalFeatureImages.at(0).size()) ); m_LabelData.reshape( vigra::Shape2(m_NumberOfSamples,1) ); m_Weights.reshape( vigra::Shape2(m_NumberOfSamples,1) ); MITK_INFO << "Number of features: " << m_FeatureData.shape(1); itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); m_RandGen->SetSeed(); MITK_INFO << "Creating training data ..."; unsigned int sampleCounter = 0; for (unsigned int t=0; tSetInputImage(fiber_folume); typename DwiFeatureImageType::Pointer image = m_DwiFeatureImages.at(t); typename DwiFeatureImageInterpolatorType::Pointer dwi_interp = DwiFeatureImageInterpolatorType::New(); dwi_interp->SetInputImage(image); ItkUcharImgType::Pointer wmMask = m_WhiteMatterImages.at(t); ItkUcharImgType::Pointer mask; if (t it(wmMask, wmMask->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { if (it.Get()==0 && (mask.IsNull() || (mask.IsNotNull() && mask->GetPixel(it.GetIndex())>0))) { typename DwiFeatureImageType::PixelType pix = image->GetPixel(it.GetIndex()); // random directions for (unsigned int i=0; iGetIntegerVariate(m_Parameters->m_NumPreviousDirections); // how many directions should be zero? for (unsigned int i=0; im_NumPreviousDirections; i++) { int c=0; vnl_vector_fixed probe; if (static_cast(i)GetVariate()*2-1; probe[1] = m_RandGen->GetVariate()*2-1; probe[2] = m_RandGen->GetVariate()*2-1; probe.normalize(); if (dot_product(ref, probe)<0) probe *= -1; } for (unsigned int f=NumberOfSignalFeatures+3*i; f itkP; image->TransformIndexToPhysicalPoint(it.GetIndex(), itkP); float v = mitk::imv::GetImageValue(itkP, false, interpolator); m_FeatureData(sampleCounter,NumberOfSignalFeatures+m_Parameters->m_NumPreviousDirections*3+add_feat_c) = v; add_feat_c++; } // label and sample weight m_LabelData(sampleCounter,0) = m_DirectionContainer.size(); m_Weights(sampleCounter,0) = 1.0; sampleCounter++; } } ++it; } unsigned int num_gm_samples = sampleCounter; // white matter samples mitk::FiberBundle::Pointer fib = m_Tractograms.at(t); vtkSmartPointer< vtkPolyData > polyData = fib->GetFiberPolyData(); vnl_vector_fixed zero_dir; zero_dir.fill(0.0); for (unsigned int i=0; iGetNumFibers(); i++) { vtkCell* cell = polyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float fiber_weight = fib->GetFiberWeight(i); for (int n = 0; n <= static_cast(m_Parameters->m_NumPreviousDirections); ++n) { if (!m_ZeroDirWmFeatures) n = m_Parameters->m_NumPreviousDirections; for (bool reverse : {false, true}) { for (int j=1; j itkP1, itkP2; int num_nonzero_dirs = m_Parameters->m_NumPreviousDirections; if (!reverse) num_nonzero_dirs = std::min(n, j); else num_nonzero_dirs = std::min(n, numPoints-j-1); vnl_vector_fixed dir; // zero directions for (unsigned int k=0; km_NumPreviousDirections-num_nonzero_dirs; k++) { dir.fill(0.0); int c = 0; for (unsigned int f=NumberOfSignalFeatures+3*k; fGetPoint(j-n_idx); itkP1[0] = p[0]; itkP1[1] = p[1]; itkP1[2] = p[2]; p = points->GetPoint(j-n_idx+1); itkP2[0] = p[0]; itkP2[1] = p[1]; itkP2[2] = p[2]; } else { p = points->GetPoint(j+n_idx); itkP1[0] = p[0]; itkP1[1] = p[1]; itkP1[2] = p[2]; p = points->GetPoint(j+n_idx-1); itkP2[0] = p[0]; itkP2[1] = p[1]; itkP2[2] = p[2]; } dir[0]=itkP1[0]-itkP2[0]; dir[1]=itkP1[1]-itkP2[1]; dir[2]=itkP1[2]-itkP2[2]; if (dir.magnitude()<0.0001) mitkThrow() << "streamline error!"; dir.normalize(); if (dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2]) mitkThrow() << "ERROR: NaN direction!"; if (dot_product(ref, dir)<0) dir *= -1; int c = 0; for (unsigned int f=NumberOfSignalFeatures+3*(k+m_Parameters->m_NumPreviousDirections-num_nonzero_dirs); fm_NumPreviousDirections-num_nonzero_dirs); f++) { m_FeatureData(sampleCounter,f) = dir[c]; c++; } } // get target direction double* p = points->GetPoint(j); itkP1[0] = p[0]; itkP1[1] = p[1]; itkP1[2] = p[2]; if (reverse) { p = points->GetPoint(j-1); itkP2[0] = p[0]; itkP2[1] = p[1]; itkP2[2] = p[2]; } else { p = points->GetPoint(j+1); itkP2[0] = p[0]; itkP2[1] = p[1]; itkP2[2] = p[2]; } dir[0]=itkP2[0]-itkP1[0]; dir[1]=itkP2[1]-itkP1[1]; dir[2]=itkP2[2]-itkP1[2]; if (dir.magnitude()<0.0001) mitkThrow() << "streamline error!"; dir.normalize(); if (dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2]) mitkThrow() << "ERROR: NaN direction!"; if (dot_product(ref, dir)<0) dir *= -1; // image features float volume_mod = mitk::imv::GetImageValue(itkP1, false, volume_interpolator); // diffusion signal features typename DwiFeatureImageType::PixelType pix = mitk::imv::GetImageValue< typename DwiFeatureImageType::PixelType >(itkP1, m_Parameters->m_InterpolateTractographyData, dwi_interp); for (unsigned int f=0; f(itkP1, false, interpolator); add_feat_c++; m_FeatureData(sampleCounter,NumberOfSignalFeatures+2+add_feat_c) = v; } // set label values float angle = 0; float m = dir.magnitude(); if (m>0.0001) { int l = 0; for (auto d : m_DirectionContainer) { float a = fabs(dot_product(dir, d)); if (a>angle) { m_LabelData(sampleCounter,0) = l; m_Weights(sampleCounter,0) = fiber_weight*volume_mod; angle = a; } l++; } } sampleCounter++; } if (!m_BidirectionalFiberSampling) // don't sample fibers backward break; } } } } m_Tractograms.clear(); MITK_INFO << "done"; } } #endif diff --git a/Modules/FiberTracking/Algorithms/itkTractDensityImageFilter.cpp b/Modules/FiberTracking/Algorithms/itkTractDensityImageFilter.cpp index c5f54d9..b7e4601 100644 --- a/Modules/FiberTracking/Algorithms/itkTractDensityImageFilter.cpp +++ b/Modules/FiberTracking/Algorithms/itkTractDensityImageFilter.cpp @@ -1,179 +1,209 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Coindex[1]right (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 "itkTractDensityImageFilter.h" // VTK #include #include #include #include // misc #include #include #include #include namespace itk{ template< class OutputImageType > TractDensityImageFilter< OutputImageType >::TractDensityImageFilter() : m_UpsamplingFactor(1) , m_InvertImage(false) - , m_BinaryOutput(false) + , m_Mode(DENSITY) , m_UseImageGeometry(false) , m_OutputAbsoluteValues(false) , m_MaxDensity(0) , m_NumCoveredVoxels(0) { } template< class OutputImageType > TractDensityImageFilter< OutputImageType >::~TractDensityImageFilter() { } +template< class OutputImageType > +TDI_MODE TractDensityImageFilter< OutputImageType >::GetMode() const +{ + return m_Mode; +} + +template< class OutputImageType > +void TractDensityImageFilter< OutputImageType >::SetMode(const TDI_MODE &Mode) +{ + m_Mode = Mode; +} + template< class OutputImageType > void TractDensityImageFilter< OutputImageType >::GenerateData() { // generate upsampled image mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry(); typename OutputImageType::Pointer outImage = this->GetOutput(); // calculate new image parameters itk::Vector newSpacing; mitk::Point3D newOrigin; itk::Matrix newDirection; ImageRegion<3> upsampledRegion; if (m_UseImageGeometry && !m_InputImage.IsNull()) { MITK_INFO << "TractDensityImageFilter: using image geometry"; newSpacing = m_InputImage->GetSpacing()/m_UpsamplingFactor; upsampledRegion = m_InputImage->GetLargestPossibleRegion(); newOrigin = m_InputImage->GetOrigin(); typename OutputImageType::RegionType::SizeType size = upsampledRegion.GetSize(); size[0] *= m_UpsamplingFactor; size[1] *= m_UpsamplingFactor; size[2] *= m_UpsamplingFactor; upsampledRegion.SetSize(size); newDirection = m_InputImage->GetDirection(); } else { MITK_INFO << "TractDensityImageFilter: using fiber bundle geometry"; newSpacing = geometry->GetSpacing()/m_UpsamplingFactor; newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); // we retrieve the origin from a vtk-polydata (corner-based) and hance have to translate it to an image geometry // i.e. center-based newOrigin[0] += bounds.GetElement(0) + 0.5 * newSpacing[0]; newOrigin[1] += bounds.GetElement(2) + 0.5 * newSpacing[1]; newOrigin[2] += bounds.GetElement(4) + 0.5 * newSpacing[2]; for (int i=0; i<3; i++) for (int j=0; j<3; j++) newDirection[j][i] = geometry->GetMatrixColumn(i)[j]; upsampledRegion.SetSize(0, ceil( geometry->GetExtent(0)*m_UpsamplingFactor ) ); upsampledRegion.SetSize(1, ceil( geometry->GetExtent(1)*m_UpsamplingFactor ) ); upsampledRegion.SetSize(2, ceil( geometry->GetExtent(2)*m_UpsamplingFactor ) ); } typename OutputImageType::RegionType::SizeType upsampledSize = upsampledRegion.GetSize(); // apply new image parameters outImage->SetSpacing( newSpacing ); outImage->SetOrigin( newOrigin ); outImage->SetDirection( newDirection ); outImage->SetLargestPossibleRegion( upsampledRegion ); outImage->SetBufferedRegion( upsampledRegion ); outImage->SetRequestedRegion( upsampledRegion ); outImage->Allocate(); outImage->FillBuffer(0.0); int w = upsampledSize[0]; int h = upsampledSize[1]; int d = upsampledSize[2]; // set/initialize output OutPixelType* outImageBufferPointer = (OutPixelType*)outImage->GetBufferPointer(); MITK_INFO << "TractDensityImageFilter: starting image generation"; vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData(); int numFibers = m_FiberBundle->GetNumFibers(); boost::progress_display disp(numFibers); for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float weight = m_FiberBundle->GetFiberWeight(i); // fill output image for( int j=0; j startVertex = mitk::imv::GetItkPoint(points->GetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; outImage->TransformPhysicalPointToIndex(startVertex, startIndex); outImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); itk::Point endVertex = mitk::imv::GetItkPoint(points->GetPoint(j + 1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; outImage->TransformPhysicalPointToIndex(endVertex, endIndex); outImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(newSpacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > segment : segments) { if (!outImage->GetLargestPossibleRegion().IsInside(segment.first)) continue; if (outImage->GetPixel(segment.first)==0) m_NumCoveredVoxels++; - if (m_BinaryOutput) + switch (m_Mode) + { + case BINARY: + { outImage->SetPixel(segment.first, 1); - else + break; + } + case VISITATION_COUNT: + { + outImage->SetPixel(segment.first, outImage->GetPixel(segment.first) + 1); + break; + } + case DENSITY: + { + outImage->SetPixel(segment.first, outImage->GetPixel(segment.first)+segment.second * weight); + break; + } + default: + { outImage->SetPixel(segment.first, outImage->GetPixel(segment.first)+segment.second * weight); + } + } } } } m_MaxDensity = 0; for (int i=0; i0) for (int i=0; i #include #include #include #include +enum TDI_MODE : int { + BINARY, + VISITATION_COUNT, + DENSITY +}; + namespace itk{ /** * \brief Generates tract density images from input fiberbundles (Calamante 2010). */ template< class OutputImageType > class TractDensityImageFilter : public ImageSource< OutputImageType > { public: typedef TractDensityImageFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef typename OutputImageType::PixelType OutPixelType; itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkTypeMacro( TractDensityImageFilter, ImageSource ) itkSetMacro( UpsamplingFactor, float) ///< use higher resolution for ouput image itkGetMacro( UpsamplingFactor, float) ///< use higher resolution for ouput image itkSetMacro( InvertImage, bool) ///< voxelvalue = 1-voxelvalue itkGetMacro( InvertImage, bool) ///< voxelvalue = 1-voxelvalue - itkSetMacro( BinaryOutput, bool) ///< generate binary fiber envelope - itkGetMacro( BinaryOutput, bool) ///< generate binary fiber envelope itkSetMacro( OutputAbsoluteValues, bool) ///< output absolute values of the number of fibers per voxel itkGetMacro( OutputAbsoluteValues, bool) ///< output absolute values of the number of fibers per voxel itkSetMacro( UseImageGeometry, bool) ///< use input image geometry to initialize output image itkGetMacro( UseImageGeometry, bool) ///< use input image geometry to initialize output image itkSetMacro( FiberBundle, mitk::FiberBundle::Pointer) ///< input fiber bundle itkSetMacro( InputImage, typename OutputImageType::Pointer) ///< use input image geometry to initialize output image itkGetMacro( MaxDensity, OutPixelType) itkGetMacro( NumCoveredVoxels, unsigned int) void GenerateData() override; + TDI_MODE GetMode() const; + void SetMode(const TDI_MODE &Mode); + protected: TractDensityImageFilter(); ~TractDensityImageFilter() override; typename OutputImageType::Pointer m_InputImage; ///< use input image geometry to initialize output image mitk::FiberBundle::Pointer m_FiberBundle; ///< input fiber bundle float m_UpsamplingFactor; ///< use higher resolution for ouput image bool m_InvertImage; ///< voxelvalue = 1-voxelvalue - bool m_BinaryOutput; ///< generate binary fiber envelope + TDI_MODE m_Mode; ///< what should the output look like bool m_UseImageGeometry; ///< use input image geometry to initialize output image bool m_OutputAbsoluteValues; ///< do not normalize image values to 0-1 bool m_UseTrilinearInterpolation; bool m_DoFiberResampling; bool m_WorkOnFiberCopy; OutPixelType m_MaxDensity; unsigned int m_NumCoveredVoxels; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkTractDensityImageFilter.cpp" #endif #endif // __itkTractDensityImageFilter_h__ diff --git a/Modules/MriSimulation/Algorithms/itkTractsToDWIImageFilter.cpp b/Modules/MriSimulation/Algorithms/itkTractsToDWIImageFilter.cpp index d263f7a..2cd0113 100644 --- a/Modules/MriSimulation/Algorithms/itkTractsToDWIImageFilter.cpp +++ b/Modules/MriSimulation/Algorithms/itkTractsToDWIImageFilter.cpp @@ -1,1772 +1,1772 @@ /*=================================================================== 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 "itkTractsToDWIImageFilter.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 #include #include #include #include #include #include #include #include #include #include #include #include namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_StatusText("") , m_UseConstantRandSeed(false) , m_RandGen(itk::Statistics::MersenneTwisterRandomVariateGenerator::New()) { m_DoubleInterpolator = itk::LinearInterpolateImageFunction< ItkDoubleImgType, float >::New(); m_NullDir.Fill(0); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >:: SimulateKspaceAcquisition( std::vector< DoubleDwiType::Pointer >& compartment_images ) { unsigned int numFiberCompartments = m_Parameters.m_FiberModelList.size(); // create slice object ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]); sliceRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]); Vector< double, 2 > sliceSpacing; sliceSpacing[0] = m_WorkingSpacing[0]; sliceSpacing[1] = m_WorkingSpacing[1]; DoubleDwiType::PixelType nullPix; nullPix.SetSize(compartment_images.at(0)->GetVectorLength()); nullPix.Fill(0.0); auto magnitudeDwiImage = DoubleDwiType::New(); magnitudeDwiImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); magnitudeDwiImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); magnitudeDwiImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); magnitudeDwiImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); magnitudeDwiImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() ); magnitudeDwiImage->Allocate(); magnitudeDwiImage->FillBuffer(nullPix); m_PhaseImage = DoubleDwiType::New(); m_PhaseImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_PhaseImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_PhaseImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_PhaseImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_PhaseImage->SetVectorLength( compartment_images.at(0)->GetVectorLength() ); m_PhaseImage->Allocate(); m_PhaseImage->FillBuffer(nullPix); m_KspaceImage = DoubleDwiType::New(); m_KspaceImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_KspaceImage->SetOrigin( m_Parameters.m_SignalGen.m_ImageOrigin ); m_KspaceImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_KspaceImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_KspaceImage->SetVectorLength( m_Parameters.m_SignalGen.m_NumberOfCoils ); m_KspaceImage->Allocate(); m_KspaceImage->FillBuffer(nullPix); // calculate coil positions double a = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(0)*m_Parameters.m_SignalGen.m_ImageSpacing[0]; double b = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1)*m_Parameters.m_SignalGen.m_ImageSpacing[1]; double c = m_Parameters.m_SignalGen.m_ImageRegion.GetSize(2)*m_Parameters.m_SignalGen.m_ImageSpacing[2]; double diagonal = sqrt(a*a+b*b)/1000; // image diagonal in m m_CoilPointset = mitk::PointSet::New(); std::vector< itk::Vector > coilPositions; itk::Vector pos; pos.Fill(0.0); pos[1] = -diagonal/2; itk::Vector center; center[0] = a/2-m_Parameters.m_SignalGen.m_ImageSpacing[0]/2; center[1] = b/2-m_Parameters.m_SignalGen.m_ImageSpacing[2]/2; center[2] = c/2-m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; for (unsigned int c=0; c temp_p; temp_p[0] = temp_v[0]; temp_p[1] = temp_v[1]; temp_p[2] = temp_v[2]; m_CoilPointset->InsertPoint(c, temp_p); double rz = 360.0/m_Parameters.m_SignalGen.m_NumberOfCoils * itk::Math::pi/180; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; pos.SetVnlVector(rotZ*pos.GetVnlVector()); } auto num_slices = compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2); auto num_gradient_volumes = static_cast(compartment_images.at(0)->GetVectorLength()); auto max_threads = omp_get_max_threads(); int out_threads = Math::ceil(std::sqrt(max_threads)); int in_threads = Math::floor(std::sqrt(max_threads)); if (out_threads > num_gradient_volumes) { out_threads = num_gradient_volumes; in_threads = Math::floor(static_cast(max_threads/out_threads)); } PrintToLog("Parallel volumes: " + boost::lexical_cast(out_threads), false, true, true); PrintToLog("Threads per slice: " + boost::lexical_cast(in_threads), false, true, true); std::list< std::tuple > spikes; if (m_Parameters.m_Misc.m_DoAddSpikes) for (unsigned int i=0; i( m_RandGen->GetIntegerVariate()%num_gradient_volumes, m_RandGen->GetIntegerVariate()%num_slices, m_RandGen->GetIntegerVariate()%m_Parameters.m_SignalGen.m_NumberOfCoils); spikes.push_back(spike); } bool output_timing = m_Parameters.m_Misc.m_OutputAdditionalImages; PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); unsigned long lastTick = 0; boost::progress_display disp(static_cast(num_gradient_volumes)*compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)); #pragma omp parallel for num_threads(out_threads) for (int g=0; gGetAbortGenerateData()) continue; std::list< std::tuple > spikeSlice; #pragma omp critical { for (auto spike : spikes) if (std::get<0>(spike) == static_cast(g)) spikeSlice.push_back(std::tuple(std::get<1>(spike), std::get<2>(spike))); } for (unsigned int z=0; z compartment_slices; std::vector< float > t2Vector; std::vector< float > t1Vector; for (unsigned int i=0; i* signalModel; if (iSetLargestPossibleRegion( sliceRegion ); slice->SetBufferedRegion( sliceRegion ); slice->SetRequestedRegion( sliceRegion ); slice->SetSpacing(sliceSpacing); slice->Allocate(); slice->FillBuffer(0.0); // extract slice from channel g for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { Float2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; slice->SetPixel(index2D, compartment_images.at(i)->GetPixel(index3D)[g]); } compartment_slices.push_back(slice); t2Vector.push_back(signalModel->GetT2()); t1Vector.push_back(signalModel->GetT1()); } if (this->GetAbortGenerateData()) continue; for (unsigned int c=0; c(ss) == z && std::get<1>(ss) == c) ++numSpikes; // create k-sapce (inverse fourier transform slices) auto idft = itk::KspaceImageFilter< Float2DImageType::PixelType >::New(); idft->SetCompartmentImages(compartment_slices); idft->SetT2(t2Vector); idft->SetT1(t1Vector); if (m_UseConstantRandSeed) { int linear_seed = g + num_gradient_volumes*z + num_gradient_volumes*compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)*c; idft->SetRandSeed(linear_seed); } idft->SetParameters(&m_Parameters); idft->SetZ((float)z-(float)( compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2) -compartment_images.at(0)->GetLargestPossibleRegion().GetSize(2)%2 ) / 2.0); idft->SetZidx(z); idft->SetCoilPosition(coilPositions.at(c)); idft->SetFiberBundle(m_FiberBundle); idft->SetTranslation(m_Translations.at(g)); idft->SetRotationMatrix(m_RotationsInv.at(g)); idft->SetDiffusionGradientDirection(m_Parameters.m_SignalGen.GetGradientDirection(g)*m_Parameters.m_SignalGen.GetBvalue()/1000.0); idft->SetSpikesPerSlice(numSpikes); idft->SetNumberOfThreads(in_threads); #pragma omp critical if (output_timing) { idft->SetStoreTimings(true); output_timing = false; } idft->Update(); #pragma omp critical if (numSpikes>0) { m_SpikeLog += "Volume " + boost::lexical_cast(g) + " Coil " + boost::lexical_cast(c) + "\n"; m_SpikeLog += idft->GetSpikeLog(); } Complex2DImageType::Pointer fSlice; fSlice = idft->GetOutput(); if (idft->GetTickImage().IsNotNull()) m_TickImage = idft->GetTickImage(); if (idft->GetRfImage().IsNotNull()) m_RfImage = idft->GetRfImage(); // fourier transform slice Complex2DImageType::Pointer newSlice; auto dft = itk::DftImageFilter< Float2DImageType::PixelType >::New(); dft->SetInput(fSlice); dft->SetParameters(m_Parameters); dft->SetNumberOfThreads(in_threads); dft->Update(); newSlice = dft->GetOutput(); // put slice back into channel g for (unsigned int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (unsigned int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; Complex2DImageType::IndexType index2D; index2D[0]=x; index2D[1]=y; Complex2DImageType::PixelType cPix = newSlice->GetPixel(index2D); double magn = sqrt(cPix.real()*cPix.real()+cPix.imag()*cPix.imag()); double phase = 0; if (cPix.real()!=0) phase = atan( cPix.imag()/cPix.real() ); DoubleDwiType::PixelType real_pix = m_OutputImagesReal.at(c)->GetPixel(index3D); real_pix[g] = cPix.real(); m_OutputImagesReal.at(c)->SetPixel(index3D, real_pix); DoubleDwiType::PixelType imag_pix = m_OutputImagesImag.at(c)->GetPixel(index3D); imag_pix[g] = cPix.imag(); m_OutputImagesImag.at(c)->SetPixel(index3D, imag_pix); DoubleDwiType::PixelType dwiPix = magnitudeDwiImage->GetPixel(index3D); DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D); if (m_Parameters.m_SignalGen.m_NumberOfCoils>1) { dwiPix[g] += magn*magn; phasePix[g] += phase*phase; } else { dwiPix[g] = magn; phasePix[g] = phase; } //#pragma omp critical { magnitudeDwiImage->SetPixel(index3D, dwiPix); m_PhaseImage->SetPixel(index3D, phasePix); // k-space image if (g==0) { DoubleDwiType::PixelType kspacePix = m_KspaceImage->GetPixel(index3D); kspacePix[c] = idft->GetKSpaceImage()->GetPixel(index2D); m_KspaceImage->SetPixel(index3D, kspacePix); } } } } if (m_Parameters.m_SignalGen.m_NumberOfCoils>1) { for (int y=0; y(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(1)); y++) for (int x=0; x(magnitudeDwiImage->GetLargestPossibleRegion().GetSize(0)); x++) { DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; DoubleDwiType::PixelType magPix = magnitudeDwiImage->GetPixel(index3D); magPix[g] = sqrt(magPix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils); DoubleDwiType::PixelType phasePix = m_PhaseImage->GetPixel(index3D); phasePix[g] = sqrt(phasePix[g]/m_Parameters.m_SignalGen.m_NumberOfCoils); //#pragma omp critical { magnitudeDwiImage->SetPixel(index3D, magPix); m_PhaseImage->SetPixel(index3D, phasePix); } } } ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) PrintToLog("*", false, false, false); lastTick = newTick; } } PrintToLog("\n", false); return magnitudeDwiImage; } template< class PixelType > TractsToDWIImageFilter< PixelType >::ItkDoubleImgType::Pointer TractsToDWIImageFilter< PixelType >:: NormalizeInsideMask(ItkDoubleImgType::Pointer image) { double max = itk::NumericTraits< double >::min(); double min = itk::NumericTraits< double >::max(); itk::ImageRegionIterator< ItkDoubleImgType > it(image, image->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { if (m_Parameters.m_SignalGen.m_MaskImage.IsNotNull() && m_Parameters.m_SignalGen.m_MaskImage->GetPixel(it.GetIndex())<=0) { it.Set(0.0); ++it; continue; } if (it.Get()>max) max = it.Get(); if (it.Get()::New(); scaler->SetInput(image); scaler->SetShift(-min); scaler->SetScale(1.0/(max-min)); scaler->Update(); return scaler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::CheckVolumeFractionImages() { m_UseRelativeNonFiberVolumeFractions = false; // check for fiber volume fraction maps unsigned int fibVolImages = 0; for (std::size_t i=0; iGetVolumeFractionImage().IsNotNull()) { PrintToLog("Using volume fraction map for fiber compartment " + boost::lexical_cast(i+1), false); fibVolImages++; } } // check for non-fiber volume fraction maps unsigned int nonfibVolImages = 0; for (std::size_t i=0; iGetVolumeFractionImage().IsNotNull()) { PrintToLog("Using volume fraction map for non-fiber compartment " + boost::lexical_cast(i+1), false); nonfibVolImages++; } } // not all fiber compartments are using volume fraction maps // --> non-fiber volume fractions are assumed to be relative to the // non-fiber volume and not absolute voxel-volume fractions. // this means if two non-fiber compartments are used but only one of them // has an associated volume fraction map, the repesctive other volume fraction map // can be determined as inverse (1-val) of the present volume fraction map- if ( fibVolImages::New(); inverter->SetMaximum(1.0); if ( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNull() && m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNotNull() ) { // m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage( // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ) ); inverter->SetInput( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage() ); inverter->Update(); m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage(inverter->GetOutput()); } else if ( m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage().IsNull() && m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage().IsNotNull() ) { // m_Parameters.m_NonFiberModelList[0]->SetVolumeFractionImage( // NormalizeInsideMask( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ) ); inverter->SetInput( m_Parameters.m_NonFiberModelList[0]->GetVolumeFractionImage() ); inverter->Update(); m_Parameters.m_NonFiberModelList[1]->SetVolumeFractionImage(inverter->GetOutput()); } else { itkExceptionMacro("Something went wrong in automatically calculating the missing non-fiber volume fraction image!" " Did you use two non fiber compartments but only one volume fraction image?" " Then it should work and this error is really strange."); } m_UseRelativeNonFiberVolumeFractions = true; nonfibVolImages++; } // Up to two fiber compartments are allowed without volume fraction maps since the volume fractions can then be determined automatically if (m_Parameters.m_FiberModelList.size()>2 && fibVolImages!=m_Parameters.m_FiberModelList.size()) itkExceptionMacro("More than two fiber compartment selected but no corresponding volume fraction maps set!"); // One non-fiber compartment is allowed without volume fraction map since the volume fraction can then be determined automatically if (m_Parameters.m_NonFiberModelList.size()>1 && nonfibVolImages!=m_Parameters.m_NonFiberModelList.size()) itkExceptionMacro("More than one non-fiber compartment selected but no volume fraction maps set!"); if (fibVolImages0) { PrintToLog("Not all fiber compartments are using an associated volume fraction image.\n" "Assuming non-fiber volume fraction images to contain values relative to the" " remaining non-fiber volume, not absolute values.", false); m_UseRelativeNonFiberVolumeFractions = true; // mitk::LocaleSwitch localeSwitch("C"); // itk::ImageFileWriter::Pointer wr = itk::ImageFileWriter::New(); // wr->SetInput(m_Parameters.m_NonFiberModelList[1]->GetVolumeFractionImage()); // wr->SetFileName("/local/volumefraction.nrrd"); // wr->Update(); } // initialize the images that store the output volume fraction of each compartment m_VolumeFractions.clear(); for (std::size_t i=0; iSetSpacing( m_WorkingSpacing ); doubleImg->SetOrigin( m_WorkingOrigin ); doubleImg->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); doubleImg->SetLargestPossibleRegion( m_WorkingImageRegion ); doubleImg->SetBufferedRegion( m_WorkingImageRegion ); doubleImg->SetRequestedRegion( m_WorkingImageRegion ); doubleImg->Allocate(); doubleImg->FillBuffer(0); m_VolumeFractions.push_back(doubleImg); } } template< class PixelType > void TractsToDWIImageFilter< PixelType >::InitializeData() { m_Rotations.clear(); m_Translations.clear(); m_MotionLog = ""; m_SpikeLog = ""; m_TickImage = nullptr; m_RfImage = nullptr; // initialize output dwi image m_Parameters.m_SignalGen.m_CroppedRegion = m_Parameters.m_SignalGen.m_ImageRegion; if (m_Parameters.m_Misc.m_DoAddAliasing) m_Parameters.m_SignalGen.m_CroppedRegion.SetSize( 1, m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1) *m_Parameters.m_SignalGen.m_CroppingFactor); itk::Point shiftedOrigin = m_Parameters.m_SignalGen.m_ImageOrigin; shiftedOrigin[1] += (m_Parameters.m_SignalGen.m_ImageRegion.GetSize(1) -m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1))*m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_OutputImage = OutputImageType::New(); m_OutputImage->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); m_OutputImage->SetOrigin( shiftedOrigin ); m_OutputImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_OutputImage->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); m_OutputImage->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); m_OutputImage->Allocate(); typename OutputImageType::PixelType temp; temp.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); temp.Fill(0.0); m_OutputImage->FillBuffer(temp); PrintToLog("Output image spacing: [" + boost::lexical_cast(m_Parameters.m_SignalGen.m_ImageSpacing[0]) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_ImageSpacing[1]) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_ImageSpacing[2]) + "]", false); PrintToLog("Output image size: [" + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(0)) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(1)) + "," + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppedRegion.GetSize(2)) + "]", false); // images containing real and imaginary part of the dMRI signal for each coil m_OutputImagesReal.clear(); m_OutputImagesImag.clear(); for (unsigned int i=0; iSetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); outputImageReal->SetOrigin( shiftedOrigin ); outputImageReal->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); outputImageReal->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageReal->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageReal->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageReal->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); outputImageReal->Allocate(); outputImageReal->FillBuffer(temp); m_OutputImagesReal.push_back(outputImageReal); typename DoubleDwiType::Pointer outputImageImag = DoubleDwiType::New(); outputImageImag->SetSpacing( m_Parameters.m_SignalGen.m_ImageSpacing ); outputImageImag->SetOrigin( shiftedOrigin ); outputImageImag->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); outputImageImag->SetLargestPossibleRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageImag->SetBufferedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageImag->SetRequestedRegion( m_Parameters.m_SignalGen.m_CroppedRegion ); outputImageImag->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); outputImageImag->Allocate(); outputImageImag->FillBuffer(temp); m_OutputImagesImag.push_back(outputImageImag); } // Apply in-plane upsampling for Gibbs ringing artifact double upsampling = 1; if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging && m_Parameters.m_SignalGen.m_ZeroRinging==0) upsampling = 2; m_WorkingSpacing = m_Parameters.m_SignalGen.m_ImageSpacing; m_WorkingSpacing[0] /= upsampling; m_WorkingSpacing[1] /= upsampling; m_WorkingImageRegion = m_Parameters.m_SignalGen.m_ImageRegion; m_WorkingImageRegion.SetSize(0, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[0]*upsampling); m_WorkingImageRegion.SetSize(1, m_Parameters.m_SignalGen.m_ImageRegion.GetSize()[1]*upsampling); m_WorkingOrigin = m_Parameters.m_SignalGen.m_ImageOrigin; m_WorkingOrigin[0] -= m_Parameters.m_SignalGen.m_ImageSpacing[0]/2; m_WorkingOrigin[0] += m_WorkingSpacing[0]/2; m_WorkingOrigin[1] -= m_Parameters.m_SignalGen.m_ImageSpacing[1]/2; m_WorkingOrigin[1] += m_WorkingSpacing[1]/2; m_WorkingOrigin[2] -= m_Parameters.m_SignalGen.m_ImageSpacing[2]/2; m_WorkingOrigin[2] += m_WorkingSpacing[2]/2; m_VoxelVolume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2]; PrintToLog("Working image spacing: [" + boost::lexical_cast(m_WorkingSpacing[0]) + "," + boost::lexical_cast(m_WorkingSpacing[1]) + "," + boost::lexical_cast(m_WorkingSpacing[2]) + "]", false); PrintToLog("Working image size: [" + boost::lexical_cast(m_WorkingImageRegion.GetSize(0)) + "," + boost::lexical_cast(m_WorkingImageRegion.GetSize(1)) + "," + boost::lexical_cast(m_WorkingImageRegion.GetSize(2)) + "]", false); // generate double images to store the individual compartment signals m_CompartmentImages.clear(); int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); for (int i=0; iSetSpacing( m_WorkingSpacing ); doubleDwi->SetOrigin( m_WorkingOrigin ); doubleDwi->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); doubleDwi->SetLargestPossibleRegion( m_WorkingImageRegion ); doubleDwi->SetBufferedRegion( m_WorkingImageRegion ); doubleDwi->SetRequestedRegion( m_WorkingImageRegion ); doubleDwi->SetVectorLength( m_Parameters.m_SignalGen.GetNumVolumes() ); doubleDwi->Allocate(); DoubleDwiType::PixelType pix; pix.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); pix.Fill(0.0); doubleDwi->FillBuffer(pix); m_CompartmentImages.push_back(doubleDwi); } if (m_FiberBundle.IsNull() && m_InputImage.IsNotNull()) { m_CompartmentImages.clear(); m_Parameters.m_SignalGen.m_DoAddMotion = false; m_Parameters.m_SignalGen.m_DoSimulateRelaxation = false; PrintToLog("Simulating acquisition for input diffusion-weighted image.", false); auto caster = itk::CastImageFilter< OutputImageType, DoubleDwiType >::New(); caster->SetInput(m_InputImage); caster->Update(); if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging && m_Parameters.m_SignalGen.m_ZeroRinging==0) { PrintToLog("Upsampling input diffusion-weighted image for Gibbs ringing simulation.", false); auto resampler = itk::ResampleDwiImageFilter< double >::New(); resampler->SetInput(caster->GetOutput()); itk::Vector< double, 3 > samplingFactor; samplingFactor[0] = upsampling; samplingFactor[1] = upsampling; samplingFactor[2] = 1; resampler->SetSamplingFactor(samplingFactor); resampler->SetInterpolation(itk::ResampleDwiImageFilter< double >::Interpolate_WindowedSinc); resampler->Update(); m_CompartmentImages.push_back(resampler->GetOutput()); } else m_CompartmentImages.push_back(caster->GetOutput()); VectorType translation; translation.Fill(0.0); MatrixType rotation; rotation.SetIdentity(); for (unsigned int g=0; gGetLargestPossibleRegion()!=m_WorkingImageRegion) { PrintToLog("Resampling tissue mask", false); // rescale mask image (otherwise there are problems with the resampling) auto rescaler = itk::RescaleIntensityImageFilter::New(); rescaler->SetInput(0,m_Parameters.m_SignalGen.m_MaskImage); rescaler->SetOutputMaximum(100); rescaler->SetOutputMinimum(0); rescaler->Update(); // resample mask image auto resampler = itk::ResampleImageFilter::New(); resampler->SetInput(rescaler->GetOutput()); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); resampler->SetOutputDirection(m_Parameters.m_SignalGen.m_ImageDirection); resampler->SetOutputStartIndex ( m_WorkingImageRegion.GetIndex() ); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); m_Parameters.m_SignalGen.m_MaskImage = resampler->GetOutput(); } // resample frequency map if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && m_Parameters.m_SignalGen.m_FrequencyMap->GetLargestPossibleRegion()!=m_WorkingImageRegion) { PrintToLog("Resampling frequency map", false); auto resampler = itk::ResampleImageFilter::New(); resampler->SetInput(m_Parameters.m_SignalGen.m_FrequencyMap); resampler->SetSize(m_WorkingImageRegion.GetSize()); resampler->SetOutputSpacing(m_WorkingSpacing); resampler->SetOutputOrigin(m_WorkingOrigin); resampler->SetOutputDirection(m_Parameters.m_SignalGen.m_ImageDirection); resampler->SetOutputStartIndex ( m_WorkingImageRegion.GetIndex() ); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); m_Parameters.m_SignalGen.m_FrequencyMap = resampler->GetOutput(); } m_MaskImageSet = true; if (m_Parameters.m_SignalGen.m_MaskImage.IsNull()) { // no input tissue mask is set -> create default PrintToLog("No tissue mask set", false); m_Parameters.m_SignalGen.m_MaskImage = ItkUcharImgType::New(); m_Parameters.m_SignalGen.m_MaskImage->SetSpacing( m_WorkingSpacing ); m_Parameters.m_SignalGen.m_MaskImage->SetOrigin( m_WorkingOrigin ); m_Parameters.m_SignalGen.m_MaskImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); m_Parameters.m_SignalGen.m_MaskImage->SetLargestPossibleRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->SetBufferedRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->SetRequestedRegion( m_WorkingImageRegion ); m_Parameters.m_SignalGen.m_MaskImage->Allocate(); m_Parameters.m_SignalGen.m_MaskImage->FillBuffer(100); m_MaskImageSet = false; } else { PrintToLog("Using tissue mask", false); } if (m_Parameters.m_SignalGen.m_DoAddMotion) { if (m_Parameters.m_SignalGen.m_DoRandomizeMotion) { PrintToLog("Random motion artifacts:", false); PrintToLog("Maximum rotation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°", false); PrintToLog("Maximum translation: +/-" + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm", false); } else { PrintToLog("Linear motion artifacts:", false); PrintToLog("Maximum rotation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Rotation) + "°", false); PrintToLog("Maximum translation: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Translation) + "mm", false); } } if ( m_Parameters.m_SignalGen.m_MotionVolumes.empty() ) { // no motion in first volume m_Parameters.m_SignalGen.m_MotionVolumes.push_back(false); // motion in all other volumes while ( m_Parameters.m_SignalGen.m_MotionVolumes.size() < m_Parameters.m_SignalGen.GetNumVolumes() ) { m_Parameters.m_SignalGen.m_MotionVolumes.push_back(true); } } // we need to know for every volume if there is motion. if this information is missing, then set corresponding fal to false while ( m_Parameters.m_SignalGen.m_MotionVolumes.size()::New(); duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage); duplicator->Update(); m_TransformedMaskImage = duplicator->GetOutput(); // second upsampling needed for motion artifacts ImageRegion<3> upsampledImageRegion = m_WorkingImageRegion; VectorType upsampledSpacing = m_WorkingSpacing; upsampledSpacing[0] /= 4; upsampledSpacing[1] /= 4; upsampledSpacing[2] /= 4; upsampledImageRegion.SetSize(0, m_WorkingImageRegion.GetSize()[0]*4); upsampledImageRegion.SetSize(1, m_WorkingImageRegion.GetSize()[1]*4); upsampledImageRegion.SetSize(2, m_WorkingImageRegion.GetSize()[2]*4); itk::Point upsampledOrigin = m_WorkingOrigin; upsampledOrigin[0] -= m_WorkingSpacing[0]/2; upsampledOrigin[0] += upsampledSpacing[0]/2; upsampledOrigin[1] -= m_WorkingSpacing[1]/2; upsampledOrigin[1] += upsampledSpacing[1]/2; upsampledOrigin[2] -= m_WorkingSpacing[2]/2; upsampledOrigin[2] += upsampledSpacing[2]/2; m_UpsampledMaskImage = ItkUcharImgType::New(); auto upsampler = itk::ResampleImageFilter::New(); upsampler->SetInput(m_Parameters.m_SignalGen.m_MaskImage); upsampler->SetOutputParametersFromImage(m_Parameters.m_SignalGen.m_MaskImage); upsampler->SetSize(upsampledImageRegion.GetSize()); upsampler->SetOutputSpacing(upsampledSpacing); upsampler->SetOutputOrigin(upsampledOrigin); auto nn_interpolator = itk::NearestNeighborInterpolateImageFunction::New(); upsampler->SetInterpolator(nn_interpolator); upsampler->Update(); m_UpsampledMaskImage = upsampler->GetOutput(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::InitializeFiberData() { m_mmRadius = m_Parameters.m_SignalGen.m_AxonRadius/1000; auto caster = itk::CastImageFilter< itk::Image, itk::Image >::New(); caster->SetInput(m_TransformedMaskImage); caster->Update(); vtkSmartPointer weights = m_FiberBundle->GetFiberWeights(); float mean_weight = 0; for (int i=0; iGetSize(); i++) mean_weight += weights->GetValue(i); mean_weight /= weights->GetSize(); if (mean_weight>0.000001) for (int i=0; iGetSize(); i++) m_FiberBundle->SetFiberWeight(i, weights->GetValue(i)/mean_weight); else PrintToLog("\nWarning: streamlines have VERY low weights. Average weight: " + boost::lexical_cast(mean_weight) + ". Possible source of calculation errors.", false, true, true); auto density_calculator = itk::TractDensityImageFilter< itk::Image >::New(); density_calculator->SetFiberBundle(m_FiberBundle); density_calculator->SetInputImage(caster->GetOutput()); - density_calculator->SetBinaryOutput(false); + density_calculator->SetMode(TDI_MODE::DENSITY); density_calculator->SetUseImageGeometry(true); density_calculator->SetOutputAbsoluteValues(true); density_calculator->Update(); double max_density = density_calculator->GetMaxDensity(); double voxel_volume = m_WorkingSpacing[0]*m_WorkingSpacing[1]*m_WorkingSpacing[2]; if (m_mmRadius>0) { std::stringstream stream; stream << std::fixed << setprecision(2) << itk::Math::pi*m_mmRadius*m_mmRadius*max_density; std::string s = stream.str(); PrintToLog("\nMax. fiber volume: " + s + "mm².", false, true, true); { double full_radius = 1000*std::sqrt(voxel_volume/(max_density*itk::Math::pi)); std::stringstream stream; stream << std::fixed << setprecision(2) << full_radius; std::string s = stream.str(); PrintToLog("\nA full fiber voxel corresponds to a fiber radius of ~" + s + "µm, given the current fiber configuration.", false, true, true); } } else { m_mmRadius = std::sqrt(voxel_volume/(max_density*itk::Math::pi)); std::stringstream stream; stream << std::fixed << setprecision(2) << m_mmRadius*1000; std::string s = stream.str(); PrintToLog("\nSetting fiber radius to " + s + "µm to obtain full voxel.", false, true, true); } // a second fiber bundle is needed to store the transformed version of the m_FiberBundleWorkingCopy m_FiberBundleTransformed = m_FiberBundle->GetDeepCopy(); } template< class PixelType > bool TractsToDWIImageFilter< PixelType >::PrepareLogFile() { if(m_Logfile.is_open()) m_Logfile.close(); std::string filePath; std::string fileName; // Get directory name: if (m_Parameters.m_Misc.m_OutputPath.size() > 0) { filePath = m_Parameters.m_Misc.m_OutputPath; if( *(--(filePath.cend())) != '/') { filePath.push_back('/'); } } else return false; // Get file name: if( ! m_Parameters.m_Misc.m_ResultNode->GetName().empty() ) { fileName = m_Parameters.m_Misc.m_ResultNode->GetName(); } else { fileName = ""; } if( ! m_Parameters.m_Misc.m_OutputPrefix.empty() ) { fileName = m_Parameters.m_Misc.m_OutputPrefix + fileName; } try { m_Logfile.open( ( filePath + '/' + fileName + ".log" ).c_str() ); } catch (const std::ios_base::failure &fail) { MITK_ERROR << "itkTractsToDWIImageFilter.cpp: Exception " << fail.what() << " while trying to open file" << filePath << '/' << fileName << ".log"; return false; } if ( m_Logfile.is_open() ) { PrintToLog( "Logfile: " + filePath + '/' + fileName + ".log", false ); return true; } else return false; } template< class PixelType > void TractsToDWIImageFilter< PixelType >::GenerateData() { PrintToLog("\n**********************************************", false); // prepare logfile PrepareLogFile(); PrintToLog("Starting Fiberfox dMRI simulation"); m_TimeProbe.Start(); // check input data if (m_FiberBundle.IsNull() && m_InputImage.IsNull()) itkExceptionMacro("Input fiber bundle and input diffusion-weighted image is nullptr!"); if (m_Parameters.m_FiberModelList.empty() && m_InputImage.IsNull()) itkExceptionMacro("No diffusion model for fiber compartments defined and input diffusion-weighted" " image is nullptr! At least one fiber compartment is necessary to simulate diffusion."); if (m_Parameters.m_NonFiberModelList.empty() && m_InputImage.IsNull()) itkExceptionMacro("No diffusion model for non-fiber compartments defined and input diffusion-weighted" " image is nullptr! At least one non-fiber compartment is necessary to simulate diffusion."); if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) // no partial volume? remove all but first fiber compartment while (m_Parameters.m_FiberModelList.size()>1) m_Parameters.m_FiberModelList.pop_back(); if (!m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition) // No upsampling of input image needed if no k-space simulation is performed m_Parameters.m_SignalGen.m_DoAddGibbsRinging = false; if (m_UseConstantRandSeed) // always generate the same random numbers? m_RandGen->SetSeed(0); else m_RandGen->SetSeed(); InitializeData(); if ( m_FiberBundle.IsNotNull() ) // if no fiber bundle is found, we directly proceed to the k-space acquisition simulation { CheckVolumeFractionImages(); InitializeFiberData(); int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); double maxVolume = 0; unsigned long lastTick = 0; int signalModelSeed = m_RandGen->GetIntegerVariate(); PrintToLog("\n", false, false); PrintToLog("Generating " + boost::lexical_cast(numFiberCompartments+numNonFiberCompartments) + "-compartment diffusion-weighted signal."); std::vector< int > bVals = m_Parameters.m_SignalGen.GetBvalues(); PrintToLog("b-values: ", false, false, true); for (auto v : bVals) PrintToLog(boost::lexical_cast(v) + " ", false, false, true); PrintToLog("\nVolumes: " + boost::lexical_cast(m_Parameters.m_SignalGen.GetNumVolumes()), false, true, true); PrintToLog("\n", false, false, true); PrintToLog("\n", false, false, true); unsigned int image_size_x = m_WorkingImageRegion.GetSize(0); unsigned int region_size_y = m_WorkingImageRegion.GetSize(1); unsigned int num_gradients = m_Parameters.m_SignalGen.GetNumVolumes(); int numFibers = m_FiberBundle->GetNumFibers(); boost::progress_display disp(numFibers*num_gradients); if (m_FiberBundle->GetMeanFiberLength()<5.0) omp_set_num_threads(2); PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); for (unsigned int g=0; gSetSeed(signalModelSeed); for (std::size_t i=0; iSetSeed(signalModelSeed); // storing voxel-wise intra-axonal volume in mm³ auto intraAxonalVolumeImage = ItkDoubleImgType::New(); intraAxonalVolumeImage->SetSpacing( m_WorkingSpacing ); intraAxonalVolumeImage->SetOrigin( m_WorkingOrigin ); intraAxonalVolumeImage->SetDirection( m_Parameters.m_SignalGen.m_ImageDirection ); intraAxonalVolumeImage->SetLargestPossibleRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetBufferedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->SetRequestedRegion( m_WorkingImageRegion ); intraAxonalVolumeImage->Allocate(); intraAxonalVolumeImage->FillBuffer(0); maxVolume = 0; double* intraAxBuffer = intraAxonalVolumeImage->GetBufferPointer(); if (this->GetAbortGenerateData()) continue; vtkPolyData* fiberPolyData = m_FiberBundleTransformed->GetFiberPolyData(); // generate fiber signal (if there are any fiber models present) if (!m_Parameters.m_FiberModelList.empty()) { std::vector< double* > buffers; for (unsigned int i=0; iGetBufferPointer()); #pragma omp parallel for for( int i=0; iGetAbortGenerateData()) continue; float fiberWeight = m_FiberBundleTransformed->GetFiberWeight(i); if (fiberWeight == 0) continue; int numPoints = -1; std::vector< itk::Vector > points_copy; #pragma omp critical { vtkCell* cell = fiberPolyData->GetCell(i); numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j))); } if (numPoints<2) continue; double seg_volume = fiberWeight*itk::Math::pi*m_mmRadius*m_mmRadius; for( int j=0; jGetAbortGenerateData()) { j=numPoints; continue; } itk::Vector v = points_copy.at(j); itk::Vector dir = points_copy.at(j+1)-v; if ( dir.GetSquaredNorm()<0.0001 || dir[0]!=dir[0] || dir[1]!=dir[1] || dir[2]!=dir[2] ) continue; dir.Normalize(); itk::Point startVertex = points_copy.at(j); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; m_TransformedMaskImage->TransformPhysicalPointToIndex(startVertex, startIndex); m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); itk::Point endVertex = points_copy.at(j+1); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; m_TransformedMaskImage->TransformPhysicalPointToIndex(endVertex, endIndex); m_TransformedMaskImage->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(m_WorkingSpacing, startIndex, endIndex, startIndexCont, endIndexCont); // generate signal for each fiber compartment for (int k=0; kSimulateMeasurement(g, dir)*seg_volume; for (std::pair< itk::Index<3>, double > seg : segments) { if (!m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(seg.first) || m_TransformedMaskImage->GetPixel(seg.first)<=0) continue; double seg_signal = seg.second*signal_add; unsigned int linear_index = g + num_gradients*seg.first[0] + num_gradients*image_size_x*seg.first[1] + num_gradients*image_size_x*region_size_y*seg.first[2]; // update dMRI volume #pragma omp atomic buffers[k][linear_index] += seg_signal; // update fiber volume image if (k==0) { linear_index = seg.first[0] + image_size_x*seg.first[1] + image_size_x*region_size_y*seg.first[2]; #pragma omp atomic intraAxBuffer[linear_index] += seg.second*seg_volume; double vol = intraAxBuffer[linear_index]; if (vol>maxVolume) { maxVolume = vol; } } } } } #pragma omp critical { // progress report ++disp; unsigned long newTick = 50*disp.count()/disp.expected_count(); for (unsigned int tick = 0; tick<(newTick-lastTick); ++tick) PrintToLog("*", false, false, false); lastTick = newTick; } } } // axon radius not manually defined --> set fullest voxel (maxVolume) to full fiber voxel double density_correctiony_global = 1.0; if (m_Parameters.m_SignalGen.m_AxonRadius<0.0001) density_correctiony_global = m_VoxelVolume/maxVolume; // generate non-fiber signal ImageRegionIterator it3(m_TransformedMaskImage, m_TransformedMaskImage->GetLargestPossibleRegion()); while(!it3.IsAtEnd()) { if (it3.Get()>0) { DoubleDwiType::IndexType index = it3.GetIndex(); double iAxVolume = intraAxonalVolumeImage->GetPixel(index); // get non-transformed point (remove headmotion tranformation) // this point lives in the volume fraction image space itk::Point volume_fraction_point; if ( m_Parameters.m_SignalGen.m_DoAddMotion ) volume_fraction_point = GetMovedPoint(index, false); else m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, volume_fraction_point); if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) { if (iAxVolume>0.0001) // scale fiber compartment to voxel { DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); pix[g] *= m_VoxelVolume/iAxVolume; m_CompartmentImages.at(0)->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(0)->SetPixel(index, 1); } else { DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); pix[g] = 0; m_CompartmentImages.at(0)->SetPixel(index, pix); SimulateExtraAxonalSignal(index, volume_fraction_point, 0, g); } } else { // manually defined axon radius and voxel overflow --> rescale to voxel volume if ( m_Parameters.m_SignalGen.m_AxonRadius>=0.0001 && iAxVolume>m_VoxelVolume ) { for (int i=0; iGetPixel(index); pix[g] *= m_VoxelVolume/iAxVolume; m_CompartmentImages.at(i)->SetPixel(index, pix); } iAxVolume = m_VoxelVolume; } // if volume fraction image is set use it, otherwise use global scaling factor double density_correction_voxel = density_correctiony_global; if ( m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr && iAxVolume>0.0001 ) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()); double volume_fraction = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator); if (volume_fraction<0) mitkThrow() << "Volume fraction image (index 1) contains negative values (intra-axonal compartment)!"; density_correction_voxel = m_VoxelVolume*volume_fraction/iAxVolume; // remove iAxVolume sclaing and scale to volume_fraction } else if (m_Parameters.m_FiberModelList[0]->GetVolumeFractionImage()!=nullptr) density_correction_voxel = 0.0; // adjust intra-axonal compartment volume by density correction factor DoubleDwiType::PixelType pix = m_CompartmentImages.at(0)->GetPixel(index); pix[g] *= density_correction_voxel; m_CompartmentImages.at(0)->SetPixel(index, pix); // normalize remaining fiber volume fractions (they are rescaled in SimulateExtraAxonalSignal) if (iAxVolume>0.0001) { for (int i=1; iGetPixel(index); pix[g] /= iAxVolume; m_CompartmentImages.at(i)->SetPixel(index, pix); } } else { for (int i=1; iGetPixel(index); pix[g] = 0; m_CompartmentImages.at(i)->SetPixel(index, pix); } } iAxVolume = density_correction_voxel*iAxVolume; // new intra-axonal volume = old intra-axonal volume * correction factor // simulate other compartments SimulateExtraAxonalSignal(index, volume_fraction_point, iAxVolume, g); } } ++it3; } } PrintToLog("\n", false); } if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } DoubleDwiType::Pointer doubleOutImage; double signalScale = m_Parameters.m_SignalGen.m_SignalScale; if ( m_Parameters.m_SignalGen.m_SimulateKspaceAcquisition ) // do k-space stuff { PrintToLog("\n", false, false); PrintToLog("Simulating k-space acquisition using " +boost::lexical_cast(m_Parameters.m_SignalGen.m_NumberOfCoils) +" coil(s)"); switch (m_Parameters.m_SignalGen.m_AcquisitionType) { case SignalGenerationParameters::SingleShotEpi: { PrintToLog("Acquisition type: single shot EPI", false); break; } case SignalGenerationParameters::ConventionalSpinEcho: { PrintToLog("Acquisition type: conventional spin echo (one RF pulse per line) with cartesian k-space trajectory", false); break; } case SignalGenerationParameters::FastSpinEcho: { PrintToLog("Acquisition type: fast spin echo (one RF pulse per ETL lines) with cartesian k-space trajectory (ETL: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_EchoTrainLength) + ")", false); break; } default: { PrintToLog("Acquisition type: single shot EPI", false); break; } } if(m_Parameters.m_SignalGen.m_tInv>0) PrintToLog("Using inversion pulse with TI " + boost::lexical_cast(m_Parameters.m_SignalGen.m_tInv) + "ms", false); if (m_Parameters.m_SignalGen.m_DoSimulateRelaxation) PrintToLog("Simulating signal relaxation", false); if (m_Parameters.m_SignalGen.m_NoiseVariance>0 && m_Parameters.m_Misc.m_DoAddNoise) PrintToLog("Simulating complex Gaussian noise: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_NoiseVariance), false); if (m_Parameters.m_SignalGen.m_FrequencyMap.IsNotNull() && m_Parameters.m_Misc.m_DoAddDistortions) PrintToLog("Simulating distortions", false); if (m_Parameters.m_SignalGen.m_DoAddGibbsRinging) { if (m_Parameters.m_SignalGen.m_ZeroRinging > 0) PrintToLog("Simulating ringing artifacts by zeroing " + boost::lexical_cast(m_Parameters.m_SignalGen.m_ZeroRinging) + "% of k-space frequencies", false); else PrintToLog("Simulating ringing artifacts by cropping high resolution inputs during k-space simulation", false); } if (m_Parameters.m_Misc.m_DoAddEddyCurrents && m_Parameters.m_SignalGen.m_EddyStrength>0) PrintToLog("Simulating eddy currents: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_EddyStrength), false); if (m_Parameters.m_Misc.m_DoAddSpikes && m_Parameters.m_SignalGen.m_Spikes>0) PrintToLog("Simulating spikes: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Spikes), false); if (m_Parameters.m_Misc.m_DoAddAliasing && m_Parameters.m_SignalGen.m_CroppingFactor<1.0) PrintToLog("Simulating aliasing: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_CroppingFactor), false); if (m_Parameters.m_Misc.m_DoAddGhosts && m_Parameters.m_SignalGen.m_KspaceLineOffset>0) PrintToLog("Simulating ghosts: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_KspaceLineOffset), false); doubleOutImage = SimulateKspaceAcquisition(m_CompartmentImages); signalScale = 1; // already scaled in SimulateKspaceAcquisition() } else // don't do k-space stuff, just sum compartments { PrintToLog("Summing compartments"); doubleOutImage = m_CompartmentImages.at(0); for (unsigned int i=1; i::New(); adder->SetInput1(doubleOutImage); adder->SetInput2(m_CompartmentImages.at(i)); adder->Update(); doubleOutImage = adder->GetOutput(); } } if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } PrintToLog("Finalizing image"); if (m_Parameters.m_SignalGen.m_DoAddDrift && m_Parameters.m_SignalGen.m_Drift>0.0) PrintToLog("Adding signal drift: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_Drift), false); if (signalScale>1) PrintToLog("Scaling signal", false); if (m_Parameters.m_NoiseModel) PrintToLog("Adding noise: " + boost::lexical_cast(m_Parameters.m_SignalGen.m_NoiseVariance), false); ImageRegionIterator it4 (m_OutputImage, m_OutputImage->GetLargestPossibleRegion()); DoubleDwiType::PixelType signal; signal.SetSize(m_Parameters.m_SignalGen.GetNumVolumes()); boost::progress_display disp2(m_OutputImage->GetLargestPossibleRegion().GetNumberOfPixels()); PrintToLog("0% 10 20 30 40 50 60 70 80 90 100%", false, true, false); PrintToLog("|----|----|----|----|----|----|----|----|----|----|\n*", false, false, false); int lastTick = 0; while(!it4.IsAtEnd()) { if (this->GetAbortGenerateData()) { PrintToLog("\n", false, false); PrintToLog("Simulation aborted"); return; } ++disp2; unsigned long newTick = 50*disp2.count()/disp2.expected_count(); for (unsigned long tick = 0; tick<(newTick-lastTick); tick++) PrintToLog("*", false, false, false); lastTick = newTick; typename OutputImageType::IndexType index = it4.GetIndex(); signal = doubleOutImage->GetPixel(index)*signalScale; for (unsigned int i=0; iAddNoise(signal); for (unsigned int i=0; i0) signal[i] = floor(signal[i]+0.5); else signal[i] = ceil(signal[i]-0.5); } it4.Set(signal); ++it4; } this->SetNthOutput(0, m_OutputImage); PrintToLog("\n", false); PrintToLog("Finished simulation"); m_TimeProbe.Stop(); if (m_Parameters.m_SignalGen.m_DoAddMotion) { PrintToLog("\nHead motion log:", false); PrintToLog(m_MotionLog, false, false); } if (m_Parameters.m_Misc.m_DoAddSpikes && m_Parameters.m_SignalGen.m_Spikes>0) { PrintToLog("\nSpike log:", false); PrintToLog(m_SpikeLog, false, false); } if (m_Logfile.is_open()) m_Logfile.close(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::PrintToLog(std::string m, bool addTime, bool linebreak, bool stdOut) { // timestamp if (addTime) { if ( m_Logfile.is_open() ) m_Logfile << this->GetTime() << " > "; m_StatusText += this->GetTime() + " > "; if (stdOut) std::cout << this->GetTime() << " > "; } // message if (m_Logfile.is_open()) m_Logfile << m; m_StatusText += m; if (stdOut) std::cout << m; // new line if (linebreak) { if (m_Logfile.is_open()) m_Logfile << "\n"; m_StatusText += "\n"; if (stdOut) std::cout << "\n"; } if ( m_Logfile.is_open() ) m_Logfile.flush(); } template< class PixelType > void TractsToDWIImageFilter< PixelType >::SimulateMotion(int g) { if ( m_Parameters.m_SignalGen.m_DoAddMotion && m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0 && m_Parameters.m_SignalGen.m_MotionVolumes[g-1]) { // The last volume was randomly moved, so we have to reset to fiberbundle and the mask. // Without motion or with linear motion, we keep the last position --> no reset. m_FiberBundleTransformed = m_FiberBundle->GetDeepCopy(); if (m_MaskImageSet) { auto duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(m_Parameters.m_SignalGen.m_MaskImage); duplicator->Update(); m_TransformedMaskImage = duplicator->GetOutput(); } } VectorType rotation; VectorType translation; // is motion artifact enabled? // is the current volume g affected by motion? if ( m_Parameters.m_SignalGen.m_DoAddMotion && m_Parameters.m_SignalGen.m_MotionVolumes[g] && g(m_Parameters.m_SignalGen.GetNumVolumes()) ) { // adjust motion transforms if ( m_Parameters.m_SignalGen.m_DoRandomizeMotion ) { // randomly rotation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[0]*2) -m_Parameters.m_SignalGen.m_Rotation[0]; rotation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[1]*2) -m_Parameters.m_SignalGen.m_Rotation[1]; rotation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Rotation[2]*2) -m_Parameters.m_SignalGen.m_Rotation[2]; translation[0] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[0]*2) -m_Parameters.m_SignalGen.m_Translation[0]; translation[1] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[1]*2) -m_Parameters.m_SignalGen.m_Translation[1]; translation[2] = m_RandGen->GetVariateWithClosedRange(m_Parameters.m_SignalGen.m_Translation[2]*2) -m_Parameters.m_SignalGen.m_Translation[2]; m_FiberBundleTransformed->TransformFibers(rotation[0], rotation[1], rotation[2], translation[0], translation[1], translation[2]); } else { // linearly rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes; translation = m_Parameters.m_SignalGen.m_Translation / m_NumMotionVolumes; m_MotionCounter++; m_FiberBundleTransformed->TransformFibers(rotation[0], rotation[1], rotation[2], translation[0], translation[1], translation[2]); rotation *= m_MotionCounter; translation *= m_MotionCounter; } MatrixType rotationMatrix = mitk::imv::GetRotationMatrixItk(rotation[0], rotation[1], rotation[2]); MatrixType rotationMatrixInv = mitk::imv::GetRotationMatrixItk(-rotation[0], -rotation[1], -rotation[2]); m_Rotations.push_back(rotationMatrix); m_RotationsInv.push_back(rotationMatrixInv); m_Translations.push_back(translation); // move mask image accoring to new transform if (m_MaskImageSet) { ImageRegionIterator maskIt(m_UpsampledMaskImage, m_UpsampledMaskImage->GetLargestPossibleRegion()); m_TransformedMaskImage->FillBuffer(0); while(!maskIt.IsAtEnd()) { if (maskIt.Get()<=0) { ++maskIt; continue; } DoubleDwiType::IndexType index = maskIt.GetIndex(); m_TransformedMaskImage->TransformPhysicalPointToIndex(GetMovedPoint(index, true), index); if (m_TransformedMaskImage->GetLargestPossibleRegion().IsInside(index)) m_TransformedMaskImage->SetPixel(index, 100); ++maskIt; } } } else { if (m_Parameters.m_SignalGen.m_DoAddMotion && !m_Parameters.m_SignalGen.m_DoRandomizeMotion && g>0) { rotation = m_Parameters.m_SignalGen.m_Rotation / m_NumMotionVolumes; rotation *= m_MotionCounter; m_Rotations.push_back(m_Rotations.back()); m_RotationsInv.push_back(m_RotationsInv.back()); m_Translations.push_back(m_Translations.back()); } else { rotation.Fill(0.0); VectorType translation; translation.Fill(0.0); MatrixType rotation_matrix; rotation_matrix.SetIdentity(); m_Rotations.push_back(rotation_matrix); m_RotationsInv.push_back(rotation_matrix); m_Translations.push_back(translation); } } if (m_Parameters.m_SignalGen.m_DoAddMotion) { m_MotionLog += boost::lexical_cast(g) + " rotation: " + boost::lexical_cast(rotation[0]) + "," + boost::lexical_cast(rotation[1]) + "," + boost::lexical_cast(rotation[2]) + ";"; m_MotionLog += " translation: " + boost::lexical_cast(m_Translations.back()[0]) + "," + boost::lexical_cast(m_Translations.back()[1]) + "," + boost::lexical_cast(m_Translations.back()[2]) + "\n"; } } template< class PixelType > itk::Point TractsToDWIImageFilter< PixelType >::GetMovedPoint(itk::Index<3>& index, bool forward) { itk::Point transformed_point; float tx = m_Translations.back()[0]; float ty = m_Translations.back()[1]; float tz = m_Translations.back()[2]; if (forward) { m_UpsampledMaskImage->TransformIndexToPhysicalPoint(index, transformed_point); m_FiberBundle->TransformPoint<>(transformed_point, m_Rotations.back(), tx, ty, tz); } else { tx *= -1; ty *= -1; tz *= -1; m_TransformedMaskImage->TransformIndexToPhysicalPoint(index, transformed_point); m_FiberBundle->TransformPoint<>(transformed_point, m_RotationsInv.back(), tx, ty, tz); } return transformed_point; } template< class PixelType > void TractsToDWIImageFilter< PixelType >:: SimulateExtraAxonalSignal(ItkUcharImgType::IndexType& index, itk::Point& volume_fraction_point, double intraAxonalVolume, int g) { int numFiberCompartments = m_Parameters.m_FiberModelList.size(); int numNonFiberCompartments = m_Parameters.m_NonFiberModelList.size(); if (m_Parameters.m_SignalGen.m_DoDisablePartialVolume) { // simulate signal for largest non-fiber compartment int max_compartment_index = 0; double max_fraction = 0; if (numNonFiberCompartments>1) { for (int i=0; iSetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()); double compartment_fraction = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator); if (compartment_fraction<0) mitkThrow() << "Volume fraction image (index " << i << ") contains values less than zero!"; if (compartment_fraction>max_fraction) { max_fraction = compartment_fraction; max_compartment_index = i; } } } DoubleDwiType::Pointer doubleDwi = m_CompartmentImages.at(max_compartment_index+numFiberCompartments); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(index); pix[g] += m_Parameters.m_NonFiberModelList[max_compartment_index]->SimulateMeasurement(g, m_NullDir)*m_VoxelVolume; doubleDwi->SetPixel(index, pix); if (g==0) m_VolumeFractions.at(max_compartment_index+numFiberCompartments)->SetPixel(index, 1); } else { std::vector< double > fractions; if (g==0) m_VolumeFractions.at(0)->SetPixel(index, intraAxonalVolume/m_VoxelVolume); double extraAxonalVolume = m_VoxelVolume-intraAxonalVolume; // non-fiber volume if (extraAxonalVolume<0) { if (extraAxonalVolume<-0.001) MITK_ERROR << "Corrupted intra-axonal signal voxel detected. Fiber volume larger voxel volume! " << m_VoxelVolume << "<" << intraAxonalVolume; extraAxonalVolume = 0; } double interAxonalVolume = 0; if (numFiberCompartments>1) interAxonalVolume = extraAxonalVolume * intraAxonalVolume/m_VoxelVolume; // inter-axonal fraction of non fiber compartment double nonFiberVolume = extraAxonalVolume - interAxonalVolume; // rest of compartment if (nonFiberVolume<0) { if (nonFiberVolume<-0.001) MITK_ERROR << "Corrupted signal voxel detected. Fiber volume larger voxel volume!"; nonFiberVolume = 0; interAxonalVolume = extraAxonalVolume; } double compartmentSum = intraAxonalVolume; fractions.push_back(intraAxonalVolume/m_VoxelVolume); // rescale extra-axonal fiber signal for (int i=1; iGetVolumeFractionImage()!=nullptr) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_FiberModelList[i]->GetVolumeFractionImage()); interAxonalVolume = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator)*m_VoxelVolume; if (interAxonalVolume<0) mitkThrow() << "Volume fraction image (index " << i+1 << ") contains negative values!"; } DoubleDwiType::PixelType pix = m_CompartmentImages.at(i)->GetPixel(index); pix[g] *= interAxonalVolume; m_CompartmentImages.at(i)->SetPixel(index, pix); compartmentSum += interAxonalVolume; fractions.push_back(interAxonalVolume/m_VoxelVolume); if (g==0) m_VolumeFractions.at(i)->SetPixel(index, interAxonalVolume/m_VoxelVolume); } for (int i=0; iGetVolumeFractionImage()!=nullptr) { m_DoubleInterpolator->SetInputImage(m_Parameters.m_NonFiberModelList[i]->GetVolumeFractionImage()); volume = mitk::imv::GetImageValue(volume_fraction_point, true, m_DoubleInterpolator)*m_VoxelVolume; if (volume<0) mitkThrow() << "Volume fraction image (index " << numFiberCompartments+i+1 << ") contains negative values (non-fiber compartment)!"; if (m_UseRelativeNonFiberVolumeFractions) volume *= nonFiberVolume/m_VoxelVolume; } DoubleDwiType::PixelType pix = m_CompartmentImages.at(i+numFiberCompartments)->GetPixel(index); pix[g] += m_Parameters.m_NonFiberModelList[i]->SimulateMeasurement(g, m_NullDir)*volume; m_CompartmentImages.at(i+numFiberCompartments)->SetPixel(index, pix); compartmentSum += volume; fractions.push_back(volume/m_VoxelVolume); if (g==0) m_VolumeFractions.at(i+numFiberCompartments)->SetPixel(index, volume/m_VoxelVolume); } if (compartmentSum/m_VoxelVolume>1.05) { MITK_ERROR << "Compartments do not sum to 1 in voxel " << index << " (" << compartmentSum/m_VoxelVolume << ")"; for (auto val : fractions) MITK_ERROR << val; } } } template< class PixelType > itk::Vector TractsToDWIImageFilter< PixelType >::GetItkVector(double point[3]) { itk::Vector itkVector; itkVector[0] = point[0]; itkVector[1] = point[1]; itkVector[2] = point[2]; return itkVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(double point[3]) { vnl_vector_fixed vnlVector; vnlVector[0] = point[0]; vnlVector[1] = point[1]; vnlVector[2] = point[2]; return vnlVector; } template< class PixelType > vnl_vector_fixed TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector& vector) { vnl_vector_fixed vnlVector; vnlVector[0] = vector[0]; vnlVector[1] = vector[1]; vnlVector[2] = vector[2]; return vnlVector; } template< class PixelType > double TractsToDWIImageFilter< PixelType >::RoundToNearest(double num) { return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5); } template< class PixelType > std::string TractsToDWIImageFilter< PixelType >::GetTime() { m_TimeProbe.Stop(); unsigned long total = RoundToNearest(m_TimeProbe.GetTotal()); unsigned long hours = total/3600; unsigned long minutes = (total%3600)/60; unsigned long seconds = total%60; std::string out = ""; out.append(boost::lexical_cast(hours)); out.append(":"); out.append(boost::lexical_cast(minutes)); out.append(":"); out.append(boost::lexical_cast(seconds)); m_TimeProbe.Start(); return out; } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp index 706cba4..25ee540 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp @@ -1,1741 +1,1741 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkFiberProcessingView.h" #include #include #include #include #include #include #include #include #include #include #include #include "usModuleRegistry.h" #include #include "mitkNodePredicateDataType.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include const std::string QmitkFiberProcessingView::VIEW_ID = "org.mitk.views.fiberprocessing"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace mitk; QmitkFiberProcessingView::QmitkFiberProcessingView() : QmitkAbstractView() , m_Controls( 0 ) , m_CircleCounter(0) , m_PolygonCounter(0) , m_UpsamplingFactor(1) { } // Destructor QmitkFiberProcessingView::~QmitkFiberProcessingView() { RemoveObservers(); } void QmitkFiberProcessingView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberProcessingViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( OnDrawCircle() ) ); connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( OnDrawPolygon() ) ); connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(GenerateAndComposite()) ); connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) ); connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) ); connect(m_Controls->m_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) ); connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) ); connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubtractBundles()) ); connect(m_Controls->m_CopyBundle, SIGNAL(clicked()), this, SLOT(CopyBundles()) ); connect(m_Controls->m_ExtractFibersButton, SIGNAL(clicked()), this, SLOT(Extract())); connect(m_Controls->m_RemoveButton, SIGNAL(clicked()), this, SLOT(Remove())); connect(m_Controls->m_ModifyButton, SIGNAL(clicked()), this, SLOT(Modify())); connect(m_Controls->m_ExtractionMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_RemovalMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ModificationMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ExtractionBoxMask, SIGNAL(currentIndexChanged(int)), this, SLOT(OnMaskExtractionChanged())); m_Controls->m_ColorMapBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isOdf = mitk::NodePredicateDataType::New("OdfImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isOdf); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); m_Controls->m_ColorMapBox->SetPredicate(finalPredicate); } UpdateGui(); OnMaskExtractionChanged(); } void QmitkFiberProcessingView::OnMaskExtractionChanged() { m_Controls->m_FiberExtractionFractionLabel->setVisible(false); m_Controls->m_FiberExtractionFractionBox->setVisible(false); m_Controls->m_FiberExtractionThresholdLabel->setVisible(false); m_Controls->m_FiberExtractionThresholdBox->setVisible(false); m_Controls->m_InterpolateRoiBox->setVisible(false); m_Controls->m_BothEnds->setVisible(false); m_Controls->m_LabelsBox->setVisible(false); m_Controls->m_LabelsLabel->setVisible(false); if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2 || m_Controls->m_ExtractionBoxMask->currentIndex() == 3) { m_Controls->m_FiberExtractionFractionLabel->setVisible(true); m_Controls->m_FiberExtractionFractionBox->setVisible(true); m_Controls->m_FiberExtractionThresholdLabel->setVisible(true); m_Controls->m_FiberExtractionThresholdBox->setVisible(true); m_Controls->m_InterpolateRoiBox->setVisible(true); } else if (m_Controls->m_ExtractionBoxMask->currentIndex() == 0 || m_Controls->m_ExtractionBoxMask->currentIndex() == 1) { if (m_Controls->m_ExtractionBoxMask->currentIndex() != 3) m_Controls->m_BothEnds->setVisible(true); m_Controls->m_InterpolateRoiBox->setVisible(true); m_Controls->m_FiberExtractionThresholdLabel->setVisible(true); m_Controls->m_FiberExtractionThresholdBox->setVisible(true); } else if (m_Controls->m_ExtractionBoxMask->currentIndex() == 4) { m_Controls->m_BothEnds->setVisible(true); m_Controls->m_LabelsBox->setVisible(true); m_Controls->m_LabelsLabel->setVisible(true); } } void QmitkFiberProcessingView::SetFocus() { m_Controls->toolBoxx->setFocus(); } void QmitkFiberProcessingView::Modify() { switch (m_Controls->m_ModificationMethodBox->currentIndex()) { case 0: { ResampleSelectedBundlesSpline(); break; } case 1: { ResampleSelectedBundlesLinear(); break; } case 2: { CompressSelectedBundles(); break; } case 3: { MirrorFibers(); break; } case 4: { WeightFibers(); break; } case 5: { DoWeightColorCoding(); break; } case 6: { DoCurvatureColorCoding(); break; } case 7: { DoLengthColorCoding(); break; } case 8: { DoImageColorCoding(); break; } } } void QmitkFiberProcessingView::WeightFibers() { float weight = this->m_Controls->m_BundleWeightBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->SetFiberWeights(weight); } } void QmitkFiberProcessingView::Remove() { switch (m_Controls->m_RemovalMethodBox->currentIndex()) { case 0: { RemoveDir(); break; } case 1: { PruneBundle(); break; } case 2: { ApplyCurvatureThreshold(); break; } case 3: { RemoveWithMask(false); break; } case 4: { RemoveWithMask(true); break; } case 5: { ApplyWeightThreshold(); break; } case 6: { ApplyDensityThreshold(); break; } } } void QmitkFiberProcessingView::Extract() { switch (m_Controls->m_ExtractionMethodBox->currentIndex()) { case 0: { ExtractWithPlanarFigure(); break; } case 1: { switch (m_Controls->m_ExtractionBoxMask->currentIndex()) { { case 0: ExtractWithMask(true, false, false); break; } { case 1: ExtractWithMask(true, true, false); break; } { case 2: ExtractWithMask(false, false, false); break; } { case 3: ExtractWithMask(false, true, false); break; } { case 4: ExtractWithMask(true, false, true); break; } } break; } } } void QmitkFiberProcessingView::PruneBundle() { int minLength = this->m_Controls->m_PruneFibersMinBox->value(); int maxLength = this->m_Controls->m_PruneFibersMaxBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); if (!fib->RemoveShortFibers(minLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); else if (!fib->RemoveLongFibers(maxLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyWeightThreshold() { float thr = this->m_Controls->m_WeightThresholdBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->FilterByWeights(thr); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByFiberWeights(false, true); node->SetData(newFib); } else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyDensityThreshold() { float thr = this->m_Controls->m_DensityThresholdBox->value(); float ol = this->m_Controls->m_DensityOverlapBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); itk::TractDensityImageFilter< ItkFloatImageType >::Pointer generator = itk::TractDensityImageFilter< ItkFloatImageType >::New(); generator->SetFiberBundle(fib); - generator->SetBinaryOutput(false); + generator->SetMode(TDI_MODE::DENSITY); generator->SetOutputAbsoluteValues(false); generator->Update(); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetRoiImages({generator->GetOutput()}); extractor->SetInputFiberBundle(fib); extractor->SetOverlapFraction(ol); extractor->SetInterpolate(true); extractor->SetThreshold(thr); extractor->SetNoNegatives(true); extractor->Update(); if (extractor->GetPositives().empty()) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } mitk::FiberBundle::Pointer newFib = extractor->GetPositives().at(0); if (newFib->GetNumFibers()>0) node->SetData(newFib); else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyCurvatureThreshold() { int angle = this->m_Controls->m_CurvSpinBox->value(); int dist = this->m_Controls->m_CurvDistanceSpinBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); itk::FiberCurvatureFilter::Pointer filter = itk::FiberCurvatureFilter::New(); filter->SetInputFiberBundle(fib); filter->SetAngularDeviation(angle); filter->SetDistance(dist); filter->SetRemoveFibers(m_Controls->m_RemoveCurvedFibersBox->isChecked()); filter->Update(); mitk::FiberBundle::Pointer newFib = filter->GetOutputFiberBundle(); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByOrientation(); node->SetData(newFib); } else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveDir() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); vnl_vector_fixed dir; dir[0] = m_Controls->m_ExtractDirX->value(); dir[1] = m_Controls->m_ExtractDirY->value(); dir[2] = m_Controls->m_ExtractDirZ->value(); fib->RemoveDir(dir,cos((float)m_Controls->m_ExtractAngle->value()*itk::Math::pi/180)); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveWithMask(bool removeInside) { if (m_RoiImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_RoiImageNode->GetData()); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(mitkMask, mask); mitk::FiberBundle::Pointer newFib = fib->RemoveFibersOutside(mask, removeInside); if (newFib->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } node->SetData(newFib); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert, bool labelmap) { if (m_RoiImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_RoiImageNode->GetData()); for (auto node : m_SelectedFB) { std::string roi_name = m_RoiImageNode->GetName(); mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(mitkMask, mask); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(fib); extractor->SetRoiImages({mask}); extractor->SetRoiImageNames({roi_name}); extractor->SetThreshold(m_Controls->m_FiberExtractionThresholdBox->value()); extractor->SetOverlapFraction(m_Controls->m_FiberExtractionFractionBox->value()); extractor->SetBothEnds(m_Controls->m_BothEnds->isChecked()); extractor->SetInterpolate(m_Controls->m_InterpolateRoiBox->isChecked()); extractor->SetMinFibersPerTract(m_Controls->m_MinExtractedFibersBox->value()); extractor->SetSplitByRoi(true); if (invert) extractor->SetNoPositives(true); else extractor->SetNoNegatives(true); if (labelmap) { std::string labels_string = m_Controls->m_LabelsBox->text().toStdString(); if (labels_string!="ALL") { std::vector strs; boost::split(strs,labels_string,boost::is_any_of(" ,;\t")); std::vector< unsigned short > labels_vector; for (auto v : strs) { try{ unsigned short l = boost::lexical_cast(v); labels_vector.push_back(l); } catch(...) { } } extractor->SetLabels(labels_vector); } extractor->SetInterpolate(false); extractor->SetInputType(itk::FiberExtractionFilter::INPUT::LABEL_MAP); extractor->SetSplitLabels(true); onlyEnds = true; } if (onlyEnds) extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); extractor->Update(); std::vector< mitk::FiberBundle::Pointer > newFibs; if (invert) newFibs = extractor->GetNegatives(); else newFibs = extractor->GetPositives(); if (newFibs.empty()) { QMessageBox::information(nullptr, "No output generated:", "No fibers could be extracted."); continue; } auto labels = extractor->GetPositiveLabels(); for (unsigned int i=0; iSetData(fib); std::string name = roi_name; if (iGetFloatProperty("Fiber2DSliceThickness", currentThickness); newNode->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(currentThickness)); GetDataStorage()->Add(newNode, node); } node->SetVisibility(false); } } void QmitkFiberProcessingView::GenerateRoiImage() { if (m_SelectedPF.empty()) return; mitk::BaseGeometry::Pointer geometry; if (!m_SelectedFB.empty()) { mitk::FiberBundle::Pointer fib = dynamic_cast(m_SelectedFB.front()->GetData()); geometry = fib->GetGeometry(); } else if (m_SelectedImage) geometry = m_SelectedImage->GetGeometry(); else return; itk::Vector spacing = geometry->GetSpacing(); spacing /= m_UpsamplingFactor; mitk::Point3D newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); newOrigin[0] += bounds.GetElement(0); newOrigin[1] += bounds.GetElement(2); newOrigin[2] += bounds.GetElement(4); itk::Matrix direction; itk::ImageRegion<3> imageRegion; for (int i=0; i<3; i++) for (int j=0; j<3; j++) direction[j][i] = geometry->GetMatrixColumn(i)[j]/spacing[j]; imageRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); imageRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); imageRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); m_PlanarFigureImage = ItkUCharImageType::New(); m_PlanarFigureImage->SetSpacing( spacing ); // Set the image spacing m_PlanarFigureImage->SetOrigin( newOrigin ); // Set the image origin m_PlanarFigureImage->SetDirection( direction ); // Set the image direction m_PlanarFigureImage->SetRegions( imageRegion ); m_PlanarFigureImage->Allocate(); m_PlanarFigureImage->FillBuffer( 0 ); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); std::string name = m_SelectedPF.at(0)->GetName(); WritePfToImage(m_SelectedPF.at(0), tmpImage); for (unsigned int i=1; iGetName(); WritePfToImage(m_SelectedPF.at(i), tmpImage); } DataNode::Pointer node = DataNode::New(); tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); node->SetData(tmpImage); node->SetName(name); this->GetDataStorage()->Add(node); } void QmitkFiberProcessingView::WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image) { if (dynamic_cast(node->GetData())) { m_PlanarFigure = dynamic_cast(node->GetData()); AccessFixedDimensionByItk_2( image, InternalReorientImagePlane, 3, m_PlanarFigure->GetGeometry(), -1); AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 3, 2, node->GetName() ); } else if (dynamic_cast(node->GetData())) { DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node); for (unsigned int i=0; iSize(); i++) { WritePfToImage(children->at(i), image); } } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< float, VImageDimension > FloatImageType; typedef itk::ResampleImageFilter ResamplerType; typename ResamplerType::Pointer resampler = ResamplerType::New(); mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); float upsamp = m_UpsamplingFactor; float gausssigma = 0.5; // Spacing typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); spacing[0] = image->GetSpacing()[0] / upsamp; spacing[1] = image->GetSpacing()[1] / upsamp; spacing[2] = image->GetSpacing()[2]; resampler->SetOutputSpacing( spacing ); // Size typename ResamplerType::SizeType size; size[0] = planegeo->GetExtentInMM(0) / spacing[0]; size[1] = planegeo->GetExtentInMM(1) / spacing[1]; size[2] = 1; resampler->SetSize( size ); // Origin typename mitk::Point3D orig = planegeo->GetOrigin(); typename mitk::Point3D corrorig; planegeo3D->WorldToIndex(orig,corrorig); corrorig[0] += 0.5/upsamp; corrorig[1] += 0.5/upsamp; corrorig[2] += 0; planegeo3D->IndexToWorld(corrorig,corrorig); resampler->SetOutputOrigin(corrorig ); // Direction typename ResamplerType::DirectionType direction; typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); for(unsigned int c=0; cSetOutputDirection( direction ); // Gaussian interpolation if(gausssigma != 0) { double sigma[3]; for( unsigned int d = 0; d < 3; d++ ) sigma[d] = gausssigma * image->GetSpacing()[d]; double alpha = 2.0; typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New(); interpolator->SetInputImage( image ); interpolator->SetParameters( sigma, alpha ); resampler->SetInterpolator( interpolator ); } else { typedef typename itk::LinearInterpolateImageFunction InterpolatorType; typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); } resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); resampler->Update(); if(additionalIndex < 0) { this->m_InternalImage = mitk::Image::New(); this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string ) { typedef itk::Image< TPixel, VImageDimension > ImageType; // Generate mask image as new image with same header as input image and // initialize with "1". ItkUCharImageType::Pointer newMaskImage = ItkUCharImageType::New(); newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); newMaskImage->Allocate(); newMaskImage->FillBuffer( 1 ); // Generate VTK polygon from (closed) PlanarFigure polyline // (The polyline points are shifted by -0.5 in z-direction to make sure // that the extrusion filter, which afterwards elevates all points by +0.5 // in z-direction, creates a 3D object which is cut by the the plane z=0) const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const BaseGeometry *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); vtkPolyData *polyline = vtkPolyData::New(); polyline->Allocate( 1, 1 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // Create VTK polydata object of polyline contour vtkPoints *points = vtkPoints::New(); PlanarFigure::PolyLineType::const_iterator it; unsigned int numberOfPoints = 0; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected image Point2D point2D = *it; planarFigurePlaneGeometry->WorldToIndex(point2D, point2D); point2D[0] -= 0.5/m_UpsamplingFactor; point2D[1] -= 0.5/m_UpsamplingFactor; planarFigurePlaneGeometry->IndexToWorld(point2D, point2D); planarFigurePlaneGeometry->Map( point2D, point3D ); // Polygons (partially) outside of the image bounds can not be processed further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { float bounds[2] = {0,0}; bounds[0] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i0); bounds[1] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i1); imageGeometry3D->WorldToIndex( point3D, point3D ); if (point3D[i0]<0) point3D[i0] = 0.0; else if (point3D[i0]>bounds[0]) point3D[i0] = bounds[0]-0.001; if (point3D[i1]<0) point3D[i1] = 0.0; else if (point3D[i1]>bounds[1]) point3D[i1] = bounds[1]-0.001; points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } else { imageGeometry3D->WorldToIndex( point3D, point3D ); // Add point to polyline array points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } } polyline->SetPoints( points ); points->Delete(); vtkIdType *ptIds = new vtkIdType[numberOfPoints]; for ( vtkIdType i = 0; i < numberOfPoints; ++i ) ptIds[i] = i; polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); // Extrude the generated contour polygon vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); extrudeFilter->SetInputData( polyline ); extrudeFilter->SetScaleFactor( 1 ); extrudeFilter->SetExtrusionTypeToNormalExtrusion(); extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); // Make a stencil from the extruded polygon vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); polyDataToImageStencil->SetInputConnection( extrudeFilter->GetOutputPort() ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< ItkUCharImageType > ImageImportType; typedef itk::VTKImageExport< ItkUCharImageType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( newMaskImage ); vtkImageImport *vtkImporter = vtkImageImport::New(); this->ConnectPipelines( itkExporter, vtkImporter ); vtkImporter->Update(); // Apply the generated image stencil to the input image vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(polyDataToImageStencil->GetOutputPort() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkImageExport *vtkExporter = vtkImageExport::New(); vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // calculate cropping bounding box m_InternalImageMask3D = itkImporter->GetOutput(); m_InternalImageMask3D->SetDirection(image->GetDirection()); itk::ImageRegionConstIterator itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); itk::ImageRegionIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); typename ImageType::SizeType lowersize = {{itk::NumericTraits::max(),itk::NumericTraits::max(),itk::NumericTraits::max()}}; typename ImageType::SizeType uppersize = {{0,0,0}}; while( !itmask.IsAtEnd() ) { if(itmask.Get() == 0) itimage.Set(0); else { typename ImageType::IndexType index = itimage.GetIndex(); typename ImageType::SizeType signedindex; signedindex[0] = index[0]; signedindex[1] = index[1]; signedindex[2] = index[2]; lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; } ++itmask; ++itimage; } typename ImageType::IndexType index; index[0] = lowersize[0]; index[1] = lowersize[1]; index[2] = lowersize[2]; typename ImageType::SizeType size; size[0] = uppersize[0] - lowersize[0] + 1; size[1] = uppersize[1] - lowersize[1] + 1; size[2] = uppersize[2] - lowersize[2] + 1; itk::ImageRegion<3> cropRegion = itk::ImageRegion<3>(index, size); // crop internal mask typedef itk::RegionOfInterestImageFilter< ItkUCharImageType, ItkUCharImageType > ROIMaskFilterType; typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); roi2->SetRegionOfInterest(cropRegion); roi2->SetInput(m_InternalImageMask3D); roi2->Update(); m_InternalImageMask3D = roi2->GetOutput(); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer()); tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer()); Image::Pointer tmpImage2 = Image::New(); tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer()); const BaseGeometry *pfImageGeometry3D = tmpImage2->GetGeometry( 0 ); const BaseGeometry *intImageGeometry3D = tmpImage->GetGeometry( 0 ); typedef itk::ImageRegionIteratorWithIndex IteratorType; IteratorType imageIterator (m_InternalImageMask3D, m_InternalImageMask3D->GetRequestedRegion()); imageIterator.GoToBegin(); while ( !imageIterator.IsAtEnd() ) { unsigned char val = imageIterator.Value(); if (val>0) { itk::Index<3> index = imageIterator.GetIndex(); Point3D point; point[0] = index[0]; point[1] = index[1]; point[2] = index[2]; intImageGeometry3D->IndexToWorld(point, point); pfImageGeometry3D->WorldToIndex(point, point); point[i0] += 0.5; point[i1] += 0.5; index[0] = point[0]; index[1] = point[1]; index[2] = point[2]; if (pfImageGeometry3D->IsIndexInside(index)) m_PlanarFigureImage->SetPixel(index, 1); } ++imageIterator; } // Clean up VTK objects polyline->Delete(); extrudeFilter->Delete(); polyDataToImageStencil->Delete(); vtkImporter->Delete(); imageStencilFilter->Delete(); //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? delete[] ptIds; } void QmitkFiberProcessingView::UpdateGui() { m_Controls->m_FibLabel->setText("mandatory"); m_Controls->m_PfLabel->setText("needed for extraction"); m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->m_RemoveButton->setEnabled(false); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); m_Controls->PFCompoANDButton->setEnabled(false); m_Controls->PFCompoORButton->setEnabled(false); m_Controls->PFCompoNOTButton->setEnabled(false); m_Controls->m_GenerateRoiImage->setEnabled(false); m_Controls->m_ExtractFibersButton->setEnabled(false); m_Controls->m_ModifyButton->setEnabled(false); m_Controls->m_CopyBundle->setEnabled(false); m_Controls->m_JoinBundles->setEnabled(false); m_Controls->m_SubstractBundles->setEnabled(false); // disable alle frames m_Controls->m_BundleWeightFrame->setVisible(false); m_Controls->m_ExtactionFramePF->setVisible(false); m_Controls->m_RemoveDirectionFrame->setVisible(false); m_Controls->m_RemoveLengthFrame->setVisible(false); m_Controls->m_RemoveCurvatureFrame->setVisible(false); m_Controls->m_RemoveByWeightFrame->setVisible(false); m_Controls->m_RemoveByDensityFrame->setVisible(false); m_Controls->m_SmoothFibersFrame->setVisible(false); m_Controls->m_CompressFibersFrame->setVisible(false); m_Controls->m_ColorFibersFrame->setVisible(false); m_Controls->m_MirrorFibersFrame->setVisible(false); m_Controls->m_MaskExtractionFrame->setVisible(false); m_Controls->m_ColorMapBox->setVisible(false); m_Controls->m_ValueAsWeightBox->setVisible(false); bool pfSelected = !m_SelectedPF.empty(); bool fibSelected = !m_SelectedFB.empty(); bool multipleFibsSelected = (m_SelectedFB.size()>1); bool maskSelected = m_RoiImageNode.IsNotNull(); bool imageSelected = m_SelectedImage.IsNotNull(); // toggle visibility of elements according to selected method switch ( m_Controls->m_ExtractionMethodBox->currentIndex() ) { case 0: m_Controls->m_ExtactionFramePF->setVisible(true); break; case 1: m_Controls->m_MaskExtractionFrame->setVisible(true); break; } switch ( m_Controls->m_RemovalMethodBox->currentIndex() ) { case 0: m_Controls->m_RemoveDirectionFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 1: m_Controls->m_RemoveLengthFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 2: m_Controls->m_RemoveCurvatureFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 3: break; case 4: break; case 5: m_Controls->m_RemoveByWeightFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 6: m_Controls->m_RemoveByDensityFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; } switch ( m_Controls->m_ModificationMethodBox->currentIndex() ) { case 0: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 1: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 2: m_Controls->m_CompressFibersFrame->setVisible(true); break; case 3: m_Controls->m_MirrorFibersFrame->setVisible(true); if (m_SelectedSurfaces.size()>0) m_Controls->m_ModifyButton->setEnabled(true); break; case 4: m_Controls->m_BundleWeightFrame->setVisible(true); break; case 5: m_Controls->m_ColorFibersFrame->setVisible(true); break; case 6: m_Controls->m_ValueAsWeightBox->setVisible(true); m_Controls->m_ColorFibersFrame->setVisible(true); break; case 7: m_Controls->m_ValueAsWeightBox->setVisible(true); m_Controls->m_ColorFibersFrame->setVisible(true); break; case 8: m_Controls->m_ValueAsWeightBox->setVisible(true); m_Controls->m_ColorFibersFrame->setVisible(true); m_Controls->m_ColorMapBox->setVisible(true); break; } // are fiber bundles selected? if ( fibSelected ) { m_Controls->m_CopyBundle->setEnabled(true); m_Controls->m_ModifyButton->setEnabled(true); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); m_Controls->m_FibLabel->setText(QString(m_SelectedFB.at(0)->GetName().c_str())); // one bundle and one planar figure needed to extract fibers if (pfSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==0) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_SelectedPF.at(0)->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } // more than two bundles needed to join/subtract if (multipleFibsSelected) { m_Controls->m_FibLabel->setText("multiple bundles selected"); m_Controls->m_JoinBundles->setEnabled(true); m_Controls->m_SubstractBundles->setEnabled(true); } if (maskSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==1) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_RoiImageNode->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } if (maskSelected && (m_Controls->m_RemovalMethodBox->currentIndex()==3 || m_Controls->m_RemovalMethodBox->currentIndex()==4) ) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_RoiImageNode->GetName().c_str())); m_Controls->m_RemoveButton->setEnabled(true); } } // are planar figures selected? if (pfSelected) { if ( fibSelected || m_SelectedImage.IsNotNull()) m_Controls->m_GenerateRoiImage->setEnabled(true); if (m_SelectedPF.size() > 1) { m_Controls->PFCompoANDButton->setEnabled(true); m_Controls->PFCompoORButton->setEnabled(true); } else m_Controls->PFCompoNOTButton->setEnabled(true); } // is image selected if (imageSelected || maskSelected) { m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); } } void QmitkFiberProcessingView::NodeRemoved(const mitk::DataNode* node ) { for (auto fnode: m_SelectedFB) if (node == fnode) { m_SelectedFB.clear(); break; } berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } void QmitkFiberProcessingView::NodeAdded(const mitk::DataNode* ) { if (!m_Controls->m_InteractiveBox->isChecked()) { berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } } void QmitkFiberProcessingView::OnEndInteraction() { if (m_Controls->m_InteractiveBox->isChecked()) ExtractWithPlanarFigure(true); } void QmitkFiberProcessingView::AddObservers() { typedef itk::SimpleMemberCommand< QmitkFiberProcessingView > SimpleCommandType; for (auto node : m_SelectedPF) { mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); if (figure!=nullptr) { figure->RemoveAllObservers(); // add observer for event when interaction with figure starts SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); endInteractionCommand->SetCallbackFunction( this, &QmitkFiberProcessingView::OnEndInteraction); m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); } } } void QmitkFiberProcessingView::RemoveObservers() { for (auto node : m_SelectedPF) { mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); if (figure!=nullptr) figure->RemoveAllObservers(); } } void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { RemoveObservers(); //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection std::vector lastSelectedFB = m_SelectedFB; m_SelectedFB.clear(); m_SelectedPF.clear(); m_SelectedSurfaces.clear(); m_SelectedImage = nullptr; m_RoiImageNode = nullptr; for (auto node: nodes) { if ( dynamic_cast(node->GetData()) ) m_SelectedFB.push_back(node); else if (dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData())) m_SelectedPF.push_back(node); else if (dynamic_cast(node->GetData())) { m_SelectedImage = dynamic_cast(node->GetData()); if (m_SelectedImage->GetDimension()==3) m_RoiImageNode = node; } else if (dynamic_cast(node->GetData())) m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); } // if we perform interactive fiber extraction, we want to avoid auto-selection of the extracted bundle if (m_SelectedFB.empty() && m_Controls->m_InteractiveBox->isChecked()) m_SelectedFB = lastSelectedFB; // if no fibers or surfaces are selected, select topmost if (m_SelectedFB.empty() && m_SelectedSurfaces.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedFB.clear(); m_SelectedFB.push_back(nodes->at(i)); } } } // if no plar figure is selected, select topmost if (m_SelectedPF.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedPF.clear(); m_SelectedPF.push_back(nodes->at(i)); } } } AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::OnDrawPolygon() { mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); figure->ClosedOn(); this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); } void QmitkFiberProcessingView::OnDrawCircle() { mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter)); } void QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *, mitk::BaseProperty* ) { // initialize figure's geometry with empty geometry mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New(); figure->SetPlaneGeometry( emptygeometry ); //set desired data to DataNode where Planarfigure is stored mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetName(name.toStdString()); newNode->SetData(figure); newNode->SetBoolProperty("planarfigure.3drendering", true); newNode->SetBoolProperty("planarfigure.3drendering.fill", true); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(newNode->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode(newNode); } // figure drawn on the topmost layer / image GetDataStorage()->Add(newNode ); RemoveObservers(); for(unsigned int i = 0; i < m_SelectedPF.size(); i++) m_SelectedPF[i]->SetSelected(false); newNode->SetSelected(true); m_SelectedPF.clear(); m_SelectedPF.push_back(newNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::ExtractWithPlanarFigure(bool interactive) { if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){ QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); return; } try { std::vector fiberBundles = m_SelectedFB; mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0); for (unsigned int i=0; i(fiberBundles.at(i)->GetData()); mitk::FiberBundle::Pointer extFB = fib->ExtractFiberSubset(planarFigure, GetDataStorage()); float currentThickness = 0; fiberBundles.at(i)->GetFloatProperty("Fiber2DSliceThickness", currentThickness); if (interactive && m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); m_InteractiveNode->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(currentThickness)); GetDataStorage()->Add(m_InteractiveNode); } float op = 5.0/sqrt(fib->GetNumFibers()); float currentOp = 0; fiberBundles.at(i)->GetFloatProperty("opacity", currentOp); if (currentOp!=op) { fib->SetFiberColors(255, 255, 255); fiberBundles.at(i)->SetFloatProperty("opacity", op); fiberBundles.at(i)->SetBoolProperty("Fiber2DfadeEFX", false); } m_InteractiveNode->SetData(extFB); } else { if (extFB->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } mitk::DataNode::Pointer node; node = mitk::DataNode::New(); node->SetData(extFB); QString name(fiberBundles.at(i)->GetName().c_str()); name += "*"; node->SetName(name.toStdString()); node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(currentThickness)); fiberBundles.at(i)->SetVisibility(false); GetDataStorage()->Add(node); } } } catch(const std::out_of_range& ) { QMessageBox::warning( nullptr, "Fiber extraction failed", "Did you only create the planar figure, using the circle or polygon button, but forgot to actually place it in the image afterwards? \nAfter creating a planar figure, simply left-click at the desired position in the image or on the tractogram to place it."); } } void QmitkFiberProcessingView::GenerateAndComposite() { mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); PFCAnd->setOperationType(mitk::PlanarFigureComposite::AND); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("AND"); newPFCNode->SetData(PFCAnd); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); RemoveObservers(); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::GenerateOrComposite() { mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); PFCOr->setOperationType(mitk::PlanarFigureComposite::OR); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("OR"); newPFCNode->SetData(PFCOr); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); UpdateGui(); } void QmitkFiberProcessingView::GenerateNotComposite() { mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); PFCNot->setOperationType(mitk::PlanarFigureComposite::NOT); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("NOT"); newPFCNode->SetData(PFCNot); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode ) { pfc->SetSelected(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(pfc, parentNode); else GetDataStorage()->Add(pfc); for (auto child : children) { if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); newChild->SetBoolProperty("planarfigure.3drendering", true); newChild->SetBoolProperty("planarfigure.3drendering.fill", true); GetDataStorage()->Add(newChild, pfc); GetDataStorage()->Remove(child); } else if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); std::vector< mitk::DataNode::Pointer > grandChildVector; mitk::DataStorage::SetOfObjects::ConstPointer grandchildren = GetDataStorage()->GetDerivations(child); for( mitk::DataStorage::SetOfObjects::const_iterator it = grandchildren->begin(); it != grandchildren->end(); ++it ) grandChildVector.push_back(*it); AddCompositeToDatastorage(newChild, grandChildVector, pfc); GetDataStorage()->Remove(child); } } UpdateGui(); } void QmitkFiberProcessingView::CopyBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( nullptr, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->GetDeepCopy(); node->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newFib); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); } UpdateGui(); } void QmitkFiberProcessingView::JoinBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } m_SelectedFB.at(0)->SetVisibility(false); mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); std::vector< mitk::FiberBundle::Pointer > tractograms; for (unsigned int i=1; iSetVisibility(false); tractograms.push_back(dynamic_cast(m_SelectedFB.at(i)->GetData())); } newBundle = newBundle->AddBundles(tractograms); mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName("Joined_Tractograms"); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } void QmitkFiberProcessingView::SubtractBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); m_SelectedFB.at(0)->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); for (unsigned int i=1; iSubtractBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); if (newBundle.IsNull()) break; name += "-"+QString(m_SelectedFB.at(i)->GetName().c_str()); m_SelectedFB.at(i)->SetVisibility(false); } if (newBundle.IsNull()) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers. Did you select the fiber bundles in the correct order? X-Y is not equal to Y-X!"); return; } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } void QmitkFiberProcessingView::ResampleSelectedBundlesSpline() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleSpline(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ResampleSelectedBundlesLinear() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleLinear(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::CompressSelectedBundles() { double factor = this->m_Controls->m_ErrorThresholdBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->Compress(factor); fib->ColorFibersByOrientation(); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::DoImageColorCoding() { if (m_Controls->m_ColorMapBox->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Bundle coloring aborted:", "No image providing the scalar values for coloring the selected bundle available."); return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByScalarMap(dynamic_cast(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked(), m_Controls->m_ValueAsWeightBox->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoCurvatureColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByCurvature(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked(), m_Controls->m_ValueAsWeightBox->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoLengthColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByLength(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked(), m_Controls->m_ValueAsWeightBox->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoWeightColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByFiberWeights(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::MirrorFibers() { unsigned int axis = this->m_Controls->m_MirrorFibersBox->currentIndex(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->MirrorFibers(axis); } for (auto surf : m_SelectedSurfaces) { vtkSmartPointer poly = surf->GetVtkPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); for (int i=0; iGetNumberOfPoints(); i++) { double* point = poly->GetPoint(i); point[axis] *= -1; vtkNewPoints->InsertNextPoint(point); } poly->SetPoints(vtkNewPoints); surf->CalculateBoundingBox(); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp index 6b4ef2f..32a7f27 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.cpp @@ -1,539 +1,543 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkFiberQuantificationView.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include // ITK #include #include #include #include #include #include #include const std::string QmitkFiberQuantificationView::VIEW_ID = "org.mitk.views.fiberquantification"; using namespace mitk; QmitkFiberQuantificationView::QmitkFiberQuantificationView() : QmitkAbstractView() , m_Controls( 0 ) , m_UpsamplingFactor(5) , m_Visible(false) { } // Destructor QmitkFiberQuantificationView::~QmitkFiberQuantificationView() { } void QmitkFiberQuantificationView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberQuantificationViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_ProcessFiberBundleButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedBundles()) ); connect( m_Controls->m_ExtractFiberPeaks, SIGNAL(clicked()), this, SLOT(CalculateFiberDirections()) ); m_Controls->m_TractBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isFib = mitk::TNodePredicateDataType::New(); m_Controls->m_TractBox->SetPredicate( isFib ); m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ImageBox->SetZeroEntryText("--"); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDimension::Pointer is3D = mitk::NodePredicateDimension::New(3); m_Controls->m_ImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, is3D) ); connect( (QObject*)(m_Controls->m_TractBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect( (QObject*)(m_Controls->m_ImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); } } void QmitkFiberQuantificationView::Activated() { } void QmitkFiberQuantificationView::Deactivated() { } void QmitkFiberQuantificationView::Visible() { m_Visible = true; } void QmitkFiberQuantificationView::Hidden() { m_Visible = false; } void QmitkFiberQuantificationView::SetFocus() { m_Controls->m_ProcessFiberBundleButton->setFocus(); } void QmitkFiberQuantificationView::CalculateFiberDirections() { typedef itk::Image ItkUcharImgType; // load fiber bundle mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(m_SelectedFB.back()->GetData()); itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); if (m_SelectedImage.IsNotNull()) { ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); mitk::CastToItkImage(m_SelectedImage, itkMaskImage); fOdfFilter->SetMaskImage(itkMaskImage); } // extract directions from fiber bundle fOdfFilter->SetFiberBundle(inputTractogram); fOdfFilter->SetAngularThreshold(cos(m_Controls->m_AngularThreshold->value()*itk::Math::pi/180)); switch (m_Controls->m_FiberDirNormBox->currentIndex()) { case 0: fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::GLOBAL_MAX); break; case 1: fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::SINGLE_VEC_NORM); break; case 2: fOdfFilter->SetNormalizationMethod(itk::TractsToVectorImageFilter::NormalizationMethods::MAX_VEC_NORM); break; } fOdfFilter->SetSizeThreshold(m_Controls->m_PeakThreshold->value()); fOdfFilter->SetMaxNumDirections(m_Controls->m_MaxNumDirections->value()); fOdfFilter->Update(); QString name = m_SelectedFB.back()->GetName().c_str(); if (m_Controls->m_NumDirectionsBox->isChecked()) { mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( fOdfFilter->GetNumDirectionsImage().GetPointer() ); mitkImage->SetVolume( fOdfFilter->GetNumDirectionsImage()->GetBufferPointer() ); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(mitkImage); node->SetName((name+"_NUM_DIRECTIONS").toStdString().c_str()); GetDataStorage()->Add(node, m_SelectedFB.back()); } Image::Pointer mitkImage = dynamic_cast(PeakImage::New().GetPointer()); mitk::CastToMitkImage(fOdfFilter->GetDirectionImage(), mitkImage); mitkImage->SetVolume(fOdfFilter->GetDirectionImage()->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(mitkImage); node->SetName( (name+"_DIRECTIONS").toStdString().c_str()); GetDataStorage()->Add(node, m_SelectedFB.back()); } void QmitkFiberQuantificationView::UpdateGui() { m_SelectedFB.clear(); if (m_Controls->m_TractBox->GetSelectedNode().IsNotNull()) m_SelectedFB.push_back(m_Controls->m_TractBox->GetSelectedNode()); m_SelectedImage = nullptr; if (m_Controls->m_ImageBox->GetSelectedNode().IsNotNull()) m_SelectedImage = dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData()); m_Controls->m_ProcessFiberBundleButton->setEnabled(!m_SelectedFB.empty()); m_Controls->m_ExtractFiberPeaks->setEnabled(!m_SelectedFB.empty()); } void QmitkFiberQuantificationView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { UpdateGui(); } void QmitkFiberQuantificationView::ProcessSelectedBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); MITK_WARN("QmitkFiberQuantificationView") << "no fibe bundle selected"; return; } int generationMethod = m_Controls->m_GenerationBox->currentIndex(); for( unsigned int i=0; i(node->GetData())) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); QString name(node->GetName().c_str()); DataNode::Pointer newNode = nullptr; switch(generationMethod){ case 0: - newNode = GenerateTractDensityImage(fib, false, true, node->GetName()); + newNode = GenerateTractDensityImage(fib, TDI_MODE::DENSITY, true, node->GetName()); name += "_TDI"; break; case 1: - newNode = GenerateTractDensityImage(fib, false, false, node->GetName()); + newNode = GenerateTractDensityImage(fib, TDI_MODE::DENSITY, false, node->GetName()); name += "_TDI"; break; case 2: - newNode = GenerateTractDensityImage(fib, true, false, node->GetName()); + newNode = GenerateTractDensityImage(fib, TDI_MODE::BINARY, false, node->GetName()); name += "_envelope"; break; case 3: newNode = GenerateColorHeatmap(fib); break; case 4: newNode = GenerateFiberEndingsImage(fib); name += "_fiber_endings"; break; case 5: newNode = GenerateFiberEndingsPointSet(fib); name += "_fiber_endings"; break; case 6: newNode = GenerateDistanceMap(fib); name += "_distance_map"; break; case 7: newNode = GenerateBinarySkeleton(fib); name += "_skeleton"; break; + case 8: + newNode = GenerateTractDensityImage(fib, TDI_MODE::VISITATION_COUNT, true, node->GetName()); + name += "_visitations"; + break; } if (newNode.IsNotNull()) { newNode->SetName(name.toStdString()); GetDataStorage()->Add(newNode); } } } } // generate pointset displaying the fiber endings mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib) { mitk::PointSet::Pointer pointSet = mitk::PointSet::New(); vtkSmartPointer fiberPolyData = fib->GetFiberPolyData(); int count = 0; int numFibers = fib->GetNumFibers(); for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints>0) { double* point = points->GetPoint(0); itk::Point itkPoint = mitk::imv::GetItkPoint(point); pointSet->InsertPoint(count, itkPoint); count++; } if (numPoints>2) { double* point = points->GetPoint(numPoints-1); itk::Point itkPoint = mitk::imv::GetItkPoint(point); pointSet->InsertPoint(count, itkPoint); count++; } } mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( pointSet ); return node; } // generate image displaying the fiber endings mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib) { typedef unsigned int OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image OutImageType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } // generate rgba heatmap from fiber bundle mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateColorHeatmap(mitk::FiberBundle::Pointer fib) { typedef itk::RGBAPixel OutPixType; typedef itk::Image OutImageType; typedef itk::TractsToRgbaImageFilter< OutImageType > ImageGeneratorType; ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); generator->SetFiberBundle(fib); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { itk::Image::Pointer itkImage = itk::Image::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateBinarySkeleton(mitk::FiberBundle::Pointer fib) { typedef itk::Image UcharImageType; itk::TractDensityImageFilter< UcharImageType >::Pointer envelope_generator = itk::TractDensityImageFilter< UcharImageType >::New(); envelope_generator->SetFiberBundle(fib); - envelope_generator->SetBinaryOutput(true); + envelope_generator->SetMode(TDI_MODE::BINARY); envelope_generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { UcharImageType::Pointer itkImage = UcharImageType::New(); CastToItkImage(m_SelectedImage, itkImage); envelope_generator->SetInputImage(itkImage); envelope_generator->SetUseImageGeometry(true); } envelope_generator->Update(); itk::BinaryThinningImageFilter::Pointer map_generator = itk::BinaryThinningImageFilter::New(); map_generator->SetInput(envelope_generator->GetOutput()); map_generator->Update(); UcharImageType::Pointer outImg = map_generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateDistanceMap(mitk::FiberBundle::Pointer fib) { typedef itk::Image UcharImageType; typedef itk::Image FloatImageType; itk::TractDensityImageFilter< UcharImageType >::Pointer envelope_generator = itk::TractDensityImageFilter< UcharImageType >::New(); envelope_generator->SetFiberBundle(fib); - envelope_generator->SetBinaryOutput(true); + envelope_generator->SetMode(TDI_MODE::BINARY); envelope_generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { UcharImageType::Pointer itkImage = UcharImageType::New(); CastToItkImage(m_SelectedImage, itkImage); envelope_generator->SetInputImage(itkImage); envelope_generator->SetUseImageGeometry(true); } envelope_generator->Update(); itk::SignedMaurerDistanceMapImageFilter::Pointer map_generator = itk::SignedMaurerDistanceMapImageFilter::New(); map_generator->SetInput(envelope_generator->GetOutput()); map_generator->SetUseImageSpacing(true); map_generator->SetSquaredDistance(false); map_generator->SetInsideIsPositive(true); map_generator->Update(); FloatImageType::Pointer outImg = map_generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); return node; } // generate tract density image from fiber bundle -mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute, std::string name) +mitk::DataNode::Pointer QmitkFiberQuantificationView::GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, TDI_MODE mode, bool absolute, std::string name) { mitk::DataNode::Pointer node = mitk::DataNode::New(); - if (binary) + if (mode==TDI_MODE::BINARY) { typedef unsigned char OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); - generator->SetBinaryOutput(binary); + generator->SetMode(mode); generator->SetOutputAbsoluteValues(absolute); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (m_SelectedImage.IsNotNull()) { mitk::LabelSetImage::Pointer multilabelImage = mitk::LabelSetImage::New(); multilabelImage->InitializeByLabeledImage(img); multilabelImage->GetActiveLabelSet()->SetActiveLabel(1); mitk::Label::Pointer label = multilabelImage->GetActiveLabel(); label->SetName("Tractogram"); // Add Segmented Property Category Code Sequence tags (0062, 0003): Sequence defining the general category of this // segment. // (0008,0100) Code Value label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_VALUE_PATH()).c_str(), TemporoSpatialStringProperty::New("T-D000A")); // (0008,0102) Coding Scheme Designator label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_SCHEME_PATH()).c_str(), TemporoSpatialStringProperty::New("SRT")); // (0008,0104) Code Meaning label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_CATEGORY_CODE_MEANING_PATH()).c_str(), TemporoSpatialStringProperty::New("Anatomical Structure")); //------------------------------------------------------------ // Add Segmented Property Type Code Sequence (0062, 000F): Sequence defining the specific property type of this // segment. // (0008,0100) Code Value label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_VALUE_PATH()).c_str(), TemporoSpatialStringProperty::New("DUMMY")); // (0008,0102) Coding Scheme Designator label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_SCHEME_PATH()).c_str(), TemporoSpatialStringProperty::New("SRT")); // (0008,0104) Code Meaning label->SetProperty( DICOMTagPathToPropertyName(DICOMSegmentationConstants::SEGMENT_TYPE_CODE_MEANING_PATH()).c_str(), TemporoSpatialStringProperty::New(name)); //Error: is undeclared// mitk::DICOMQIPropertyHandler::DeriveDICOMSourceProperties(m_SelectedImage, multilabelImage); // init data node node->SetData(multilabelImage); } else { // init data node node->SetData(img); } } else { typedef float OutPixType; typedef itk::Image OutImageType; itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); generator->SetFiberBundle(fib); - generator->SetBinaryOutput(binary); + generator->SetMode(mode); generator->SetOutputAbsoluteValues(absolute); generator->SetUpsamplingFactor(m_Controls->m_UpsamplingSpinBox->value()); if (m_SelectedImage.IsNotNull()) { OutImageType::Pointer itkImage = OutImageType::New(); CastToItkImage(m_SelectedImage, itkImage); generator->SetInputImage(itkImage); generator->SetUseImageGeometry(true); } //generator->SetDoFiberResampling(false); generator->Update(); // get output image typedef itk::Image OutType; OutType::Pointer outImg = generator->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node node->SetData(img); } return node; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.h index bb3d808..53d0c61 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationView.h @@ -1,89 +1,90 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef QmitkFiberQuantificationView_h #define QmitkFiberQuantificationView_h #include #include "ui_QmitkFiberQuantificationViewControls.h" #include #include #include #include +#include /*! \brief Generation of images from fiber bundles (TDI, envelopes, endpoint distribution) and extraction of principal fiber directions from tractograms. */ class QmitkFiberQuantificationView : public QmitkAbstractView, public mitk::ILifecycleAwarePart { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: typedef itk::Image< unsigned char, 3 > itkUCharImageType; static const std::string VIEW_ID; QmitkFiberQuantificationView(); virtual ~QmitkFiberQuantificationView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; virtual void Activated() override; virtual void Deactivated() override; virtual void Visible() override; virtual void Hidden() override; protected slots: void ProcessSelectedBundles(); ///< start selected operation on fiber bundle (e.g. tract density image generation) void CalculateFiberDirections(); ///< Calculate main fiber directions from tractogram void UpdateGui(); ///< update button activity etc. dpending on current datamanager selection protected: /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkFiberQuantificationViewControls* m_Controls; std::vector m_SelectedFB; ///< selected fiber bundle nodes mitk::Image::Pointer m_SelectedImage; float m_UpsamplingFactor; ///< upsampling factor for all image generations - mitk::DataNode::Pointer GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, bool binary, bool absolute, std::string name); + mitk::DataNode::Pointer GenerateTractDensityImage(mitk::FiberBundle::Pointer fib, TDI_MODE mode, bool absolute, std::string name); mitk::DataNode::Pointer GenerateColorHeatmap(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsImage(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateFiberEndingsPointSet(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateDistanceMap(mitk::FiberBundle::Pointer fib); mitk::DataNode::Pointer GenerateBinarySkeleton(mitk::FiberBundle::Pointer fib); bool m_Visible; }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui index c427969..c2cbb6e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberQuantificationViewControls.ui @@ -1,430 +1,435 @@ QmitkFiberQuantificationViewControls 0 0 365 581 Form QCommandLinkButton:disabled { border: none; } QGroupBox { background-color: transparent; } 25 Fiber-derived images 6 6 6 6 false 0 0 200 16777215 11 Perform selected operation on all selected fiber bundles. Generate Image 0 0 Upsampling factor 1 0.100000000000000 10.000000000000000 0.100000000000000 1.000000000000000 0 0 Tract Density Image (TDI) Normalized TDI Binary Envelope Fiber Bundle Image Fiber Endings Image Fiber Endings Pointset Distance Map Binary Skeleton + + + Streamline Visitation Count + + Principal Fiber Directions 6 6 6 6 QFrame::NoFrame QFrame::Raised 0 0 0 0 0 0 Fiber directions with an angle smaller than the defined threshold are clustered. 2 0.000000000000000 90.000000000000000 1.000000000000000 30.000000000000000 0 0 <html><head/><body><p>Directions shorter than the defined threshold are discarded.</p></body></html> 3 1.000000000000000 0.100000000000000 0.300000000000000 Angular Threshold: Max. Peaks: Size Threshold: 0 0 Maximum number of fiber directions per voxel. 100 3 Normalization: 0 0 0 Global maximum Single vector Voxel-wise maximum 0 0 Image containing the number of distinct fiber clusters per voxel. Output #Directions per Voxel false false Generate Directions Input Data 6 6 6 6 Tractogram: Reference Image: Qt::Vertical 20 40 QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h
 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp index 0726d1f..7dc0143 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp @@ -1,577 +1,689 @@ /*=================================================================== 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 "QmitkTractometryView.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include +#include const std::string QmitkTractometryView::VIEW_ID = "org.mitk.views.tractometry"; using namespace mitk; QmitkTractometryView::QmitkTractometryView() : QmitkAbstractView() , m_Controls( nullptr ) , m_Visible(false) { } // Destructor QmitkTractometryView::~QmitkTractometryView() { } void QmitkTractometryView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkTractometryViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_MethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_StartButton, SIGNAL(clicked()), this, SLOT(StartTractometry()) ); mitk::TNodePredicateDataType::Pointer imageP = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDimension::Pointer dimP = mitk::NodePredicateDimension::New(3); m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ImageBox->SetPredicate(mitk::NodePredicateAnd::New(imageP, dimP)); m_Controls->m_ChartWidget->SetXAxisLabel("Tract position"); m_Controls->m_ChartWidget->SetYAxisLabel("Image Value"); m_Controls->m_ChartWidget->SetTheme(QmitkChartWidget::ColorTheme::darkstyle); } } void QmitkTractometryView::SetFocus() { } void QmitkTractometryView::UpdateGui() { berry::IWorkbenchPart::Pointer nullPart; OnSelectionChanged(nullPart, QList(m_CurrentSelection)); } bool QmitkTractometryView::Flip(vtkSmartPointer< vtkPolyData > polydata1, int i, vtkSmartPointer< vtkPolyData > ref_poly) { double d_direct = 0; double d_flipped = 0; vtkCell* cell1 = polydata1->GetCell(0); if (ref_poly!=nullptr) cell1 = ref_poly->GetCell(0); auto numPoints1 = cell1->GetNumberOfPoints(); vtkPoints* points1 = cell1->GetPoints(); std::vector> ref_points; for (int j=0; jGetPoint(j); itk::Point itk_p; itk_p[0] = p1[0]; itk_p[1] = p1[1]; itk_p[2] = p1[2]; ref_points.push_back(itk_p); } vtkCell* cell2 = polydata1->GetCell(i); vtkPoints* points2 = cell2->GetPoints(); for (int j=0; jGetPoint(j); d_direct = (p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) + (p1[2]-p2[2])*(p1[2]-p2[2]); double* p3 = points2->GetPoint(numPoints1-j-1); d_flipped = (p1[0]-p3[0])*(p1[0]-p3[0]) + (p1[1]-p3[1])*(p1[1]-p3[1]) + (p1[2]-p3[2])*(p1[2]-p3[2]); } if (d_direct>d_flipped) return true; return false; } template void QmitkTractometryView::StaticResamplingTractometry(const mitk::PixelType, mitk::Image::Pointer image, mitk::DataNode::Pointer node, std::vector > &data, std::string& clipboard_string) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); unsigned int num_points = m_Controls->m_SamplingPointsBox->value(); mitk::ImagePixelReadAccessor readimage(image, image->GetVolumeData(0)); mitk::FiberBundle::Pointer working_fib = fib->GetDeepCopy(); working_fib->ResampleToNumPoints(num_points); vtkSmartPointer< vtkPolyData > polydata = working_fib->GetFiberPolyData(); double rgb[3] = {0,0,0}; mitk::LookupTable::Pointer lookupTable = mitk::LookupTable::New(); lookupTable->SetType(mitk::LookupTable::MULTILABEL); std::vector > all_values; std::vector< double > mean_values; for (unsigned int i=0; iGetNumFibers(); ++i) { vtkCell* cell = polydata->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); std::vector< double > fib_vals; bool flip = false; if (i>0) flip = Flip(polydata, i); else if (m_ReferencePolyData!=nullptr) flip = Flip(polydata, 0, m_ReferencePolyData); for (int j=0; jGetTableValue(j, rgb); double* p; if (flip) { auto p_idx = numPoints - j - 1; p = points->GetPoint(p_idx); working_fib->ColorSinglePoint(i, p_idx, rgb); } else { p = points->GetPoint(j); working_fib->ColorSinglePoint(i, j, rgb); } Point3D px; px[0] = p[0]; px[1] = p[1]; px[2] = p[2]; double pixelValue = static_cast(readimage.GetPixelByWorldCoordinates(px)); fib_vals.push_back(pixelValue); mean += pixelValue; if (pixelValuemax) max = pixelValue; mean_values.at(j) += pixelValue; } all_values.push_back(fib_vals); } if (m_ReferencePolyData==nullptr) m_ReferencePolyData = polydata; std::vector< double > std_values1; std::vector< double > std_values2; for (unsigned int i=0; iGetNumFibers(); double stdev = 0; for (unsigned int j=0; j(mean_values.at(i)); clipboard_string += " "; clipboard_string += boost::lexical_cast(stdev); clipboard_string += "\n"; } clipboard_string += "\n"; data.push_back(mean_values); data.push_back(std_values1); data.push_back(std_values2); MITK_INFO << "Min: " << min; MITK_INFO << "Max: " << max; MITK_INFO << "Mean: " << mean/working_fib->GetNumberOfPoints(); if (m_Controls->m_ShowBinned->isChecked()) { mitk::DataNode::Pointer new_node = mitk::DataNode::New(); auto children = GetDataStorage()->GetDerivations(node); for (unsigned int i=0; isize(); ++i) { if (children->at(i)->GetName() == "binned_static") { new_node = children->at(i); new_node->SetData(working_fib); new_node->SetVisibility(true); node->SetVisibility(false); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); return; } } new_node->SetData(working_fib); new_node->SetName("binned_static"); new_node->SetVisibility(true); node->SetVisibility(false); GetDataStorage()->Add(new_node, node); } } template void QmitkTractometryView::NearestCentroidPointTractometry(const mitk::PixelType, mitk::Image::Pointer image, mitk::DataNode::Pointer node, std::vector< std::vector< double > >& data, std::string& clipboard_string) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); unsigned int num_points = m_Controls->m_SamplingPointsBox->value(); mitk::ImagePixelReadAccessor readimage(image, image->GetVolumeData(0)); mitk::FiberBundle::Pointer working_fib = fib->GetDeepCopy(); working_fib->ResampleSpline(1.0); vtkSmartPointer< vtkPolyData > working_polydata = working_fib->GetFiberPolyData(); // clustering std::vector< mitk::ClusteringMetric* > metrics; metrics.push_back({new mitk::ClusteringMetricEuclideanStd()}); mitk::FiberBundle::Pointer fib_static_resampled = fib->GetDeepCopy(); fib_static_resampled->ResampleToNumPoints(num_points); vtkSmartPointer< vtkPolyData > polydata_static_resampled = fib_static_resampled->GetFiberPolyData(); std::vector centroids; std::shared_ptr< mitk::TractClusteringFilter > clusterer = std::make_shared(); int c=0; while (c<30 && (centroids.empty() || centroids.size()>static_cast(m_Controls->m_MaxCentroids->value()))) { float cluster_size = m_Controls->m_ClusterSize->value() + m_Controls->m_ClusterSize->value()*c*0.2; float max_d = 0; int i=1; std::vector< float > distances; while (max_d < working_fib->GetGeometry()->GetDiagonalLength()/2) { distances.push_back(cluster_size*i); max_d = cluster_size*i; ++i; } clusterer->SetDistances(distances); clusterer->SetTractogram(fib_static_resampled); clusterer->SetMetrics(metrics); clusterer->SetMergeDuplicateThreshold(cluster_size); clusterer->SetDoResampling(false); clusterer->SetNumPoints(num_points); // clusterer->SetMaxClusters(m_Controls->m_MaxCentroids->value()); clusterer->SetMinClusterSize(1); clusterer->Update(); centroids = clusterer->GetOutCentroids(); ++c; } double rgb[3] = {0,0,0}; mitk::LookupTable::Pointer lookupTable = mitk::LookupTable::New(); lookupTable->SetType(mitk::LookupTable::MULTILABEL); std::vector > all_values; std::vector< double > mean_values; std::vector< unsigned int > value_count; for (unsigned int i=0; iGetNumFibers(); ++i) { vtkCell* cell = working_polydata->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); std::vector< double > fib_vals; for (int j=0; jGetPoint(j); int min_bin = 0; float d=999999; for (auto centroid : centroids) { auto centroid_polydata = centroid->GetFiberPolyData(); vtkCell* centroid_cell = centroid_polydata->GetCell(0); auto centroid_numPoints = centroid_cell->GetNumberOfPoints(); vtkPoints* centroid_points = centroid_cell->GetPoints(); bool centroid_flip = Flip(centroid_polydata, 0, centroids.at(0)->GetFiberPolyData()); for (int bin=0; binGetPoint(bin); float temp_d = std::sqrt((p[0]-centroid_p[0])*(p[0]-centroid_p[0]) + (p[1]-centroid_p[1])*(p[1]-centroid_p[1]) + (p[2]-centroid_p[2])*(p[2]-centroid_p[2])); if (temp_dGetTableValue(min_bin, rgb); working_fib->ColorSinglePoint(i, j, rgb); Point3D px; px[0] = p[0]; px[1] = p[1]; px[2] = p[2]; double pixelValue = static_cast(readimage.GetPixelByWorldCoordinates(px)); fib_vals.push_back(pixelValue); mean += pixelValue; if (pixelValuemax) max = pixelValue; mean_values.at(min_bin) += pixelValue; value_count.at(min_bin) += 1; } all_values.push_back(fib_vals); } if (m_ReferencePolyData==nullptr) m_ReferencePolyData = working_polydata; std::vector< double > std_values1; std::vector< double > std_values2; for (unsigned int i=0; i(mean_values.at(i)); clipboard_string += " "; clipboard_string += boost::lexical_cast(stdev); clipboard_string += "\n"; } clipboard_string += "\n"; data.push_back(mean_values); data.push_back(std_values1); data.push_back(std_values2); MITK_INFO << "Min: " << min; MITK_INFO << "Max: " << max; MITK_INFO << "Mean: " << mean/working_fib->GetNumberOfPoints(); if (m_Controls->m_ShowBinned->isChecked()) { mitk::DataNode::Pointer new_node = mitk::DataNode::New(); // mitk::DataNode::Pointer new_node2; auto children = GetDataStorage()->GetDerivations(node); for (unsigned int i=0; isize(); ++i) { if (children->at(i)->GetName() == "binned_centroid") { new_node = children->at(i); new_node->SetData(working_fib); new_node->SetVisibility(true); node->SetVisibility(false); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); return; } } new_node->SetData(working_fib); new_node->SetName("binned_centroid"); new_node->SetVisibility(true); node->SetVisibility(false); GetDataStorage()->Add(new_node, node); } } void QmitkTractometryView::Activated() { } void QmitkTractometryView::Deactivated() { } void QmitkTractometryView::Visible() { m_Visible = true; QList selection = GetDataManagerSelection(); berry::IWorkbenchPart::Pointer nullPart; OnSelectionChanged(nullPart, selection); } void QmitkTractometryView::Hidden() { m_Visible = false; } std::string QmitkTractometryView::RGBToHexString(double *rgb) { std::ostringstream os; for (int i = 0; i < 3; ++i) { os << std::setw(2) << std::setfill('0') << std::hex << static_cast(rgb[i] * 255); } return os.str(); } +void QmitkTractometryView::AlongTractRadiomicsPreprocessing(mitk::Image::Pointer image, mitk::DataNode::Pointer node) +{ + mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); + + // calculate mask + typedef unsigned int OutPixType; + typedef itk::Image OutImageType; + itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); + generator->SetFiberBundle(fib); + generator->SetMode(TDI_MODE::BINARY); + OutImageType::Pointer itkImage = OutImageType::New(); + CastToItkImage(image, itkImage); + generator->SetInputImage(itkImage); + generator->SetUseImageGeometry(true); + generator->Update(); + OutImageType::Pointer count_map = generator->GetOutput(); + + unsigned int num_points = m_Controls->m_SamplingPointsBox->value(); + mitk::FiberBundle::Pointer working_fib = fib->GetDeepCopy(); + working_fib->ResampleToNumPoints(num_points); + vtkSmartPointer< vtkPolyData > polydata = working_fib->GetFiberPolyData(); + + + itk::ImageRegionIterator< OutImageType > it(count_map, count_map->GetLargestPossibleRegion()); + while( !it.IsAtEnd() ) + { + if (it.Get()>0) + { + std::vector seg_vote; seg_vote.resize(num_points, 0); + typename OutImageType::PointType image_point; + count_map->TransformIndexToPhysicalPoint(it.GetIndex(), image_point); + + for (unsigned int i=0; iGetNumFibers(); ++i) + { + vtkCell* cell = polydata->GetCell(i); + auto numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + bool flip = false; + if (i>0) + flip = Flip(polydata, i); + else if (m_ReferencePolyData!=nullptr) + flip = Flip(polydata, 0, m_ReferencePolyData); + + float local_d = 99999999; + int local_closest_seg = -1; + + for (int j=0; jGetPoint(segment_id); + } + else + { + p = points->GetPoint(j); + segment_id = j; + } + + float d = std::fabs( (p[0]-image_point[0]) ) + std::fabs( (p[1]-image_point[1]) ) + std::fabs( (p[2]-image_point[2]) ); + if (d=max_count) + { + final_seg_id = i; + max_count = seg_vote.at(i); + } + } + + it.Set(final_seg_id + 1); + } + ++it; + } + + + mitk::Image::Pointer seg_img = mitk::Image::New(); + seg_img->InitializeByItk(count_map.GetPointer()); + seg_img->SetVolume(count_map->GetBufferPointer()); + + mitk::DataNode::Pointer new_node = mitk::DataNode::New(); + new_node->SetData(seg_img); + new_node->SetName("segment image"); + new_node->SetVisibility(true); + node->SetVisibility(false); + GetDataStorage()->Add(new_node, node); +} + void QmitkTractometryView::StartTractometry() { m_ReferencePolyData = nullptr; double color[3] = {0,0,0}; mitk::LookupTable::Pointer lookupTable = mitk::LookupTable::New(); lookupTable->SetType(mitk::LookupTable::MULTILABEL); mitk::Image::Pointer image = dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData()); this->m_Controls->m_ChartWidget->Clear(); std::string clipboardString = ""; int c = 1; for (auto node : m_CurrentSelection) { clipboardString += node->GetName() + "\n"; clipboardString += "mean stdev\n"; std::vector< std::vector< double > > data; - if (m_Controls->m_MethodBox->currentIndex()==0) + switch (m_Controls->m_MethodBox->currentIndex()) + { + case 1: { mitkPixelTypeMultiplex4( StaticResamplingTractometry, image->GetPixelType(), image, node, data, clipboardString ); + break; + } + case 2: + { + AlongTractRadiomicsPreprocessing(image, node); + return; } - else + default: { mitkPixelTypeMultiplex4( NearestCentroidPointTractometry, image->GetPixelType(), image, node, data, clipboardString ); } - + } m_Controls->m_ChartWidget->AddData1D(data.at(0), node->GetName() + " Mean", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->SetLineStyle(node->GetName() + " Mean", QmitkChartWidget::LineStyle::solid); if (m_Controls->m_StDevBox->isChecked()) { m_Controls->m_ChartWidget->AddData1D(data.at(1), node->GetName() + " +STDEV", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->AddData1D(data.at(2), node->GetName() + " -STDEV", QmitkChartWidget::ChartType::line); m_Controls->m_ChartWidget->SetLineStyle(node->GetName() + " +STDEV", QmitkChartWidget::LineStyle::dashed); m_Controls->m_ChartWidget->SetLineStyle(node->GetName() + " -STDEV", QmitkChartWidget::LineStyle::dashed); } lookupTable->GetTableValue(c, color); this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " Mean", RGBToHexString(color)); if (m_Controls->m_StDevBox->isChecked()) { color[0] *= 0.5; color[1] *= 0.5; color[2] *= 0.5; this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " +STDEV", RGBToHexString(color)); this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " -STDEV", RGBToHexString(color)); } this->m_Controls->m_ChartWidget->Show(true); this->m_Controls->m_ChartWidget->SetShowDataPoints(false); ++c; } QApplication::clipboard()->setText(clipboardString.c_str(), QClipboard::Clipboard); } void QmitkTractometryView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { m_Controls->m_StartButton->setEnabled(false); if (!m_Visible) return; if (m_Controls->m_MethodBox->currentIndex()==0) m_Controls->m_ClusterFrame->setVisible(false); else m_Controls->m_ClusterFrame->setVisible(true); m_CurrentSelection.clear(); if(m_Controls->m_ImageBox->GetSelectedNode().IsNull()) return; for (auto node: nodes) if ( dynamic_cast(node->GetData()) ) m_CurrentSelection.push_back(node); if (!m_CurrentSelection.empty()) m_Controls->m_StartButton->setEnabled(true); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.h index e85a6ed..e610038 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.h @@ -1,84 +1,86 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef QmitkTractometryView_h #define QmitkTractometryView_h #include #include "ui_QmitkTractometryViewControls.h" #include #include #include #include #include #include /*! \brief Weight fibers by linearly fitting them to the image data. */ class QmitkTractometryView : public QmitkAbstractView, public mitk::ILifecycleAwarePart { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; QmitkTractometryView(); virtual ~QmitkTractometryView(); virtual void CreateQtPartControl(QWidget *parent) override; template void StaticResamplingTractometry(const mitk::PixelType, mitk::Image::Pointer image, mitk::DataNode::Pointer node, std::vector< std::vector< double > >& data, std::string& clipboard_string); template void NearestCentroidPointTractometry(const mitk::PixelType, mitk::Image::Pointer image, mitk::DataNode::Pointer node, std::vector< std::vector< double > >& data, std::string& clipboard_string); + void AlongTractRadiomicsPreprocessing(mitk::Image::Pointer image, mitk::DataNode::Pointer node); + virtual void SetFocus() override; virtual void Activated() override; virtual void Deactivated() override; virtual void Visible() override; virtual void Hidden() override; protected slots: void UpdateGui(); void StartTractometry(); protected: /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkTractometryViewControls* m_Controls; bool Flip(vtkSmartPointer< vtkPolyData > polydata1, int i, vtkSmartPointer ref_poly=nullptr); std::string RGBToHexString(double *rgb); vtkSmartPointer< vtkPolyData > m_ReferencePolyData; QList m_CurrentSelection; bool m_Visible; }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryViewControls.ui index 5123bbe..23a578c 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryViewControls.ui @@ -1,217 +1,222 @@ QmitkTractometryViewControls 0 0 484 - 744 + 922 Form QCommandLinkButton:disabled { border: none; } QGroupBox { background-color: transparent; } Qt::Vertical 20 40 Show STDEV true Show binned tract true 0 600 Centroid Options Cluster Size: Max. Number of Centroids: 1 99 1 15.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 6 3 99999 20 Input Image: Sampling Points: Static Resampling Centroid Based + + + Segment voting + + Binning Method: Start Tractometry QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkChartWidget QWidget
QmitkChartWidget.h