diff --git a/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp b/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp index ec0fc3d..d25f851 100644 --- a/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp @@ -1,584 +1,582 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include const int numOdfSamples = 200; typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType; /*! \brief */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform streamline tractography"); parser.setContributor("MIC"); // parameters fo all methods parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::StringList, "Input:", "input image (multiple possible for 'DetTensor' algorithm)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output fiberbundle/probability map", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("type", "", mitkDiffusionCommandLineParser::String, "Type:", "which tracker to use (Peaks; Tensor; ODF; ODF-DIPY/FSL; RF)", us::Any(), false); parser.addArgument("probabilistic", "", mitkDiffusionCommandLineParser::Bool, "Probabilistic:", "Probabilistic tractography", us::Any(false)); parser.endGroup(); parser.beginGroup("2. Seeding:"); parser.addArgument("seeds", "", mitkDiffusionCommandLineParser::Int, "Seeds per voxel:", "number of seed points per voxel", 1); parser.addArgument("seed_image", "", mitkDiffusionCommandLineParser::String, "Seed image:", "mask image defining seed voxels", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("trials_per_seed", "", mitkDiffusionCommandLineParser::Int, "Max. trials per seed:", "try each seed N times until a valid streamline is obtained (only for probabilistic tractography)", 10); parser.addArgument("max_tracts", "", mitkDiffusionCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded", -1); parser.endGroup(); parser.beginGroup("3. Tractography constraints:"); parser.addArgument("tracking_mask", "", mitkDiffusionCommandLineParser::String, "Mask image:", "streamlines leaving the mask will stop immediately", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("stop_image", "", mitkDiffusionCommandLineParser::String, "Stop ROI image:", "streamlines entering the mask will stop immediately", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("exclusion_image", "", mitkDiffusionCommandLineParser::String, "Exclusion ROI image:", "streamlines entering the mask will be discarded", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("ep_constraint", "", mitkDiffusionCommandLineParser::String, "Endpoint constraint:", "determines which fibers are accepted based on their endpoint location - options are NONE, EPS_IN_TARGET, EPS_IN_TARGET_LABELDIFF, EPS_IN_SEED_AND_TARGET, MIN_ONE_EP_IN_TARGET, ONE_EP_IN_TARGET and NO_EP_IN_TARGET", us::Any()); parser.addArgument("target_image", "", mitkDiffusionCommandLineParser::String, "Target ROI image:", "effact depends on the chosen endpoint constraint (option ep_constraint)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.endGroup(); parser.beginGroup("4. Streamline integration parameters:"); parser.addArgument("sharpen_odfs", "", mitkDiffusionCommandLineParser::Bool, "SHarpen ODFs:", "if you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). this is not necessary for CSD fODFs, since they are narurally much sharper."); parser.addArgument("cutoff", "", mitkDiffusionCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 0.1); parser.addArgument("odf_cutoff", "", mitkDiffusionCommandLineParser::Float, "ODF Cutoff:", "threshold on the ODF magnitude. this is useful in case of CSD fODF tractography.", 0.0); parser.addArgument("step_size", "", mitkDiffusionCommandLineParser::Float, "Step size:", "step size (in voxels)", 0.5); parser.addArgument("min_tract_length", "", mitkDiffusionCommandLineParser::Float, "Min. tract length:", "minimum fiber length (in mm)", 20); parser.addArgument("angular_threshold", "", mitkDiffusionCommandLineParser::Float, "Angular threshold:", "angular threshold between two successive steps, (default: 90° * step_size, minimum 15°)"); parser.addArgument("loop_check", "", mitkDiffusionCommandLineParser::Float, "Check for loops:", "threshold on angular stdev over the last 4 voxel lengths"); parser.addArgument("peak_jitter", "", mitkDiffusionCommandLineParser::Float, "Peak jitter:", "important for probabilistic peak tractography and peak prior. actual jitter is drawn from a normal distribution with peak_jitter*fabs(direction_value) as standard deviation.", 0.01); parser.endGroup(); parser.beginGroup("5. Tractography prior:"); parser.addArgument("prior_image", "", mitkDiffusionCommandLineParser::String, "Peak prior:", "tractography prior in thr for of a peak image", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("prior_weight", "", mitkDiffusionCommandLineParser::Float, "Prior weight", "weighting factor between prior and data.", 0.5); parser.addArgument("dont_restrict_to_prior", "", mitkDiffusionCommandLineParser::Bool, "Don't restrict to prior:", "don't restrict tractography to regions where the prior is valid.", us::Any(false)); parser.addArgument("no_new_directions_from_prior", "", mitkDiffusionCommandLineParser::Bool, "No new directios from prior:", "the prior cannot create directions where there are none in the data.", us::Any(false)); parser.addArgument("prior_flip_x", "", mitkDiffusionCommandLineParser::Bool, "Prior Flip X:", "multiply x-coordinate of prior direction by -1"); parser.addArgument("prior_flip_y", "", mitkDiffusionCommandLineParser::Bool, "Prior Flip Y:", "multiply y-coordinate of prior direction by -1"); parser.addArgument("prior_flip_z", "", mitkDiffusionCommandLineParser::Bool, "Prior Flip Z:", "multiply z-coordinate of prior direction by -1"); parser.endGroup(); parser.beginGroup("6. Neighborhood sampling:"); parser.addArgument("num_samples", "", mitkDiffusionCommandLineParser::Int, "Num. neighborhood samples:", "number of neighborhood samples that are use to determine the next progression direction", 0); parser.addArgument("sampling_distance", "", mitkDiffusionCommandLineParser::Float, "Sampling distance:", "distance of neighborhood sampling points (in voxels)", 0.25); parser.addArgument("use_stop_votes", "", mitkDiffusionCommandLineParser::Bool, "Use stop votes:", "use stop votes"); parser.addArgument("use_only_forward_samples", "", mitkDiffusionCommandLineParser::Bool, "Use only forward samples:", "use only forward samples"); parser.endGroup(); parser.beginGroup("7. Tensor tractography specific:"); parser.addArgument("tend_f", "", mitkDiffusionCommandLineParser::Float, "Weight f", "weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0); parser.addArgument("tend_g", "", mitkDiffusionCommandLineParser::Float, "Weight g", "weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0); parser.endGroup(); parser.beginGroup("8. Random forest tractography specific:"); parser.addArgument("forest", "", mitkDiffusionCommandLineParser::String, "Forest:", "input random forest (HDF5 file)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("use_sh_features", "", mitkDiffusionCommandLineParser::Bool, "Use SH features:", "use SH features"); parser.endGroup(); parser.beginGroup("9. Additional input:"); parser.addArgument("additional_images", "", mitkDiffusionCommandLineParser::StringList, "Additional images:", "specify a list of float images that hold additional information (FA, GFA, additional features for RF tractography)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.endGroup(); parser.beginGroup("10. Misc:"); parser.addArgument("flip_x", "", mitkDiffusionCommandLineParser::Bool, "Flip X:", "multiply x-coordinate of direction proposal by -1"); parser.addArgument("flip_y", "", mitkDiffusionCommandLineParser::Bool, "Flip Y:", "multiply y-coordinate of direction proposal by -1"); parser.addArgument("flip_z", "", mitkDiffusionCommandLineParser::Bool, "Flip Z:", "multiply z-coordinate of direction proposal by -1"); parser.addArgument("no_data_interpolation", "", mitkDiffusionCommandLineParser::Bool, "Don't interpolate input data:", "don't interpolate input image values"); parser.addArgument("no_mask_interpolation", "", mitkDiffusionCommandLineParser::Bool, "Don't interpolate masks:", "don't interpolate mask image values"); parser.addArgument("compress", "", mitkDiffusionCommandLineParser::Bool, "Compress:", "compress output fibers (lossy)"); parser.addArgument("fix_seed", "", mitkDiffusionCommandLineParser::Bool, "Fix Random Seed:", "always use the same random numbers"); parser.addArgument("parameter_file", "", mitkDiffusionCommandLineParser::String, "Parameter File:", "load parameters from json file (svae using MITK Diffusion GUI). the parameters loaded form this file are overwritten by the manually set parameters.", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["i"]); std::string outFile = us::any_cast(parsedArgs["o"]); std::string type = us::any_cast(parsedArgs["type"]); std::shared_ptr< mitk::StreamlineTractographyParameters > params = std::make_shared(); if (parsedArgs.count("parameter_file")) { auto parameter_file = us::any_cast(parsedArgs["parameter_file"]); params->LoadParameters(parameter_file); } if (parsedArgs.count("probabilistic")) params->m_Mode = mitk::StreamlineTractographyParameters::MODE::PROBABILISTIC; else { params->m_Mode = mitk::StreamlineTractographyParameters::MODE::DETERMINISTIC; } std::string prior_image = ""; if (parsedArgs.count("prior_image")) prior_image = us::any_cast(parsedArgs["prior_image"]); if (parsedArgs.count("prior_weight")) params->m_Weight = us::any_cast(parsedArgs["prior_weight"]); if (parsedArgs.count("fix_seed")) params->m_FixRandomSeed = us::any_cast(parsedArgs["fix_seed"]); if (parsedArgs.count("dont_restrict_to_prior")) params->m_RestrictToPrior = !us::any_cast(parsedArgs["dont_restrict_to_prior"]); if (parsedArgs.count("no_new_directions_from_prior")) params->m_NewDirectionsFromPrior = !us::any_cast(parsedArgs["no_new_directions_from_prior"]); if (parsedArgs.count("sharpen_odfs")) params->m_SharpenOdfs = us::any_cast(parsedArgs["sharpen_odfs"]); if (parsedArgs.count("no_data_interpolation")) params->m_InterpolateTractographyData = !us::any_cast(parsedArgs["no_data_interpolation"]); params->m_InterpolateRoiImages = true; if (parsedArgs.count("no_mask_interpolation")) params->m_InterpolateRoiImages = !us::any_cast(parsedArgs["no_mask_interpolation"]); bool use_sh_features = false; if (parsedArgs.count("use_sh_features")) use_sh_features = us::any_cast(parsedArgs["use_sh_features"]); if (parsedArgs.count("use_stop_votes")) params->m_StopVotes = us::any_cast(parsedArgs["use_stop_votes"]); if (parsedArgs.count("use_only_forward_samples")) params->m_OnlyForwardSamples = us::any_cast(parsedArgs["use_only_forward_samples"]); if (parsedArgs.count("flip_x")) params->m_FlipX = us::any_cast(parsedArgs["flip_x"]); if (parsedArgs.count("flip_y")) params->m_FlipY = us::any_cast(parsedArgs["flip_y"]); if (parsedArgs.count("flip_z")) params->m_FlipZ = us::any_cast(parsedArgs["flip_z"]); if (parsedArgs.count("prior_flip_x")) params->m_PriorFlipX = us::any_cast(parsedArgs["prior_flip_x"]); if (parsedArgs.count("prior_flip_y")) params->m_PriorFlipY = us::any_cast(parsedArgs["prior_flip_y"]); if (parsedArgs.count("prior_flip_z")) params->m_PriorFlipZ = us::any_cast(parsedArgs["prior_flip_z"]); if (parsedArgs.count("apply_image_rotation")) params->m_ApplyDirectionMatrix = us::any_cast(parsedArgs["apply_image_rotation"]); if (parsedArgs.count("compress")) params->m_CompressFibers = us::any_cast(parsedArgs["compress"]); if (parsedArgs.count("min_tract_length")) params->m_MinTractLengthMm = us::any_cast(parsedArgs["min_tract_length"]); if (parsedArgs.count("loop_check")) params->SetLoopCheckDeg(us::any_cast(parsedArgs["loop_check"])); std::string forestFile; if (parsedArgs.count("forest")) forestFile = us::any_cast(parsedArgs["forest"]); std::string maskFile = ""; if (parsedArgs.count("tracking_mask")) maskFile = us::any_cast(parsedArgs["tracking_mask"]); std::string seedFile = ""; if (parsedArgs.count("seed_image")) seedFile = us::any_cast(parsedArgs["seed_image"]); std::string targetFile = ""; if (parsedArgs.count("target_image")) targetFile = us::any_cast(parsedArgs["target_image"]); std::string exclusionFile = ""; if (parsedArgs.count("exclusion_image")) exclusionFile = us::any_cast(parsedArgs["exclusion_image"]); std::string stopFile = ""; if (parsedArgs.count("stop_image")) stopFile = us::any_cast(parsedArgs["stop_image"]); if (parsedArgs.count("ep_constraint")) { if (us::any_cast(parsedArgs["ep_constraint"]) == "NONE") params->m_EpConstraints = mitk::StreamlineTractographyParameters::EndpointConstraints::NONE; else if (us::any_cast(parsedArgs["ep_constraint"]) == "EPS_IN_TARGET") params->m_EpConstraints = mitk::StreamlineTractographyParameters::EndpointConstraints::EPS_IN_TARGET; else if (us::any_cast(parsedArgs["ep_constraint"]) == "EPS_IN_TARGET_LABELDIFF") params->m_EpConstraints = mitk::StreamlineTractographyParameters::EndpointConstraints::EPS_IN_TARGET_LABELDIFF; else if (us::any_cast(parsedArgs["ep_constraint"]) == "EPS_IN_SEED_AND_TARGET") params->m_EpConstraints = mitk::StreamlineTractographyParameters::EndpointConstraints::EPS_IN_SEED_AND_TARGET; else if (us::any_cast(parsedArgs["ep_constraint"]) == "MIN_ONE_EP_IN_TARGET") params->m_EpConstraints = mitk::StreamlineTractographyParameters::EndpointConstraints::MIN_ONE_EP_IN_TARGET; else if (us::any_cast(parsedArgs["ep_constraint"]) == "ONE_EP_IN_TARGET") params->m_EpConstraints = mitk::StreamlineTractographyParameters::EndpointConstraints::ONE_EP_IN_TARGET; else if (us::any_cast(parsedArgs["ep_constraint"]) == "NO_EP_IN_TARGET") params->m_EpConstraints = mitk::StreamlineTractographyParameters::EndpointConstraints::NO_EP_IN_TARGET; } if (parsedArgs.count("cutoff")) params->m_Cutoff = us::any_cast(parsedArgs["cutoff"]); if (parsedArgs.count("odf_cutoff")) params->m_OdfCutoff = us::any_cast(parsedArgs["odf_cutoff"]); if (parsedArgs.count("peak_jitter")) params->m_PeakJitter = us::any_cast(parsedArgs["peak_jitter"]); if (parsedArgs.count("step_size")) params->SetStepSizeVox(us::any_cast(parsedArgs["step_size"])); if (parsedArgs.count("sampling_distance")) params->SetSamplingDistanceVox(us::any_cast(parsedArgs["sampling_distance"])); if (parsedArgs.count("num_samples")) params->m_NumSamples = static_cast(us::any_cast(parsedArgs["num_samples"])); if (parsedArgs.count("seeds")) params->m_SeedsPerVoxel = us::any_cast(parsedArgs["seeds"]); if (parsedArgs.count("trials_per_seed")) params->m_TrialsPerSeed = static_cast(us::any_cast(parsedArgs["trials_per_seed"])); if (parsedArgs.count("tend_f")) params->m_F = us::any_cast(parsedArgs["tend_f"]); if (parsedArgs.count("tend_g")) params->m_G = us::any_cast(parsedArgs["tend_g"]); if (parsedArgs.count("angular_threshold")) params->SetAngularThresholdDeg(us::any_cast(parsedArgs["angular_threshold"])); if (parsedArgs.count("max_tracts")) params->m_MaxNumFibers = us::any_cast(parsedArgs["max_tracts"]); std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); if (ext != ".fib" && ext != ".trk") { MITK_INFO << "Output file format not supported. Use one of .fib, .trk, .nii, .nii.gz, .nrrd"; return EXIT_FAILURE; } // LOAD DATASETS mitkDiffusionCommandLineParser::StringContainerType addFiles; if (parsedArgs.count("additional_images")) addFiles = us::any_cast(parsedArgs["additional_images"]); typedef itk::Image ItkFloatImgType; ItkFloatImgType::Pointer mask = nullptr; if (!maskFile.empty()) { MITK_INFO << "loading mask image"; mitk::Image::Pointer img = mitk::IOUtil::Load(maskFile); mask = ItkFloatImgType::New(); mitk::CastToItkImage(img, mask); } ItkFloatImgType::Pointer seed = nullptr; if (!seedFile.empty()) { MITK_INFO << "loading seed ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(seedFile); seed = ItkFloatImgType::New(); mitk::CastToItkImage(img, seed); } ItkFloatImgType::Pointer stop = nullptr; if (!stopFile.empty()) { MITK_INFO << "loading stop ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(stopFile); stop = ItkFloatImgType::New(); mitk::CastToItkImage(img, stop); } ItkFloatImgType::Pointer target = nullptr; if (!targetFile.empty()) { MITK_INFO << "loading target ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(targetFile); target = ItkFloatImgType::New(); mitk::CastToItkImage(img, target); } ItkFloatImgType::Pointer exclusion = nullptr; if (!exclusionFile.empty()) { MITK_INFO << "loading exclusion ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(exclusionFile); exclusion = ItkFloatImgType::New(); mitk::CastToItkImage(img, exclusion); } MITK_INFO << "loading additional images"; std::vector< std::vector< ItkFloatImgType::Pointer > > addImages; addImages.push_back(std::vector< ItkFloatImgType::Pointer >()); for (auto file : addFiles) { mitk::Image::Pointer img = mitk::IOUtil::Load(file); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); addImages.at(0).push_back(itkimg); } // ////////////////////////////////////////////////////////////////// // omp_set_num_threads(1); typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); if (!prior_image.empty()) { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, std::vector()); mitk::PeakImage::Pointer priorImage = mitk::IOUtil::Load(prior_image, &functor); if (priorImage.IsNull()) { MITK_INFO << "Only peak images are supported as prior at the moment!"; return EXIT_FAILURE; } mitk::TrackingDataHandler* priorhandler = new mitk::TrackingHandlerPeaks(); typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(priorImage); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); std::shared_ptr< mitk::StreamlineTractographyParameters > prior_params = std::make_shared< mitk::StreamlineTractographyParameters >(*params); prior_params->m_FlipX = params->m_PriorFlipX; prior_params->m_FlipY = params->m_PriorFlipY; prior_params->m_FlipZ = params->m_PriorFlipZ; prior_params->m_Cutoff = 0.0; dynamic_cast(priorhandler)->SetPeakImage(itkImg); priorhandler->SetParameters(prior_params); tracker->SetTrackingPriorHandler(priorhandler); } mitk::TrackingDataHandler* handler; mitk::Image::Pointer reference_image; if (type == "RF") { mitk::TractographyForest::Pointer forest = mitk::IOUtil::Load(forestFile); if (forest.IsNull()) mitkThrow() << "Forest file " << forestFile << " could not be read."; std::vector include = {"Diffusion Weighted Images"}; std::vector exclude = {}; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor(include, exclude); auto input = mitk::IOUtil::Load(input_files.at(0), &functor); reference_image = input; if (use_sh_features) { handler = new mitk::TrackingHandlerRandomForest<6,28>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } else { handler = new mitk::TrackingHandlerRandomForest<6,100>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } } else if (type == "Peaks") { handler = new mitk::TrackingHandlerPeaks(); MITK_INFO << "loading input peak image"; mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(input_files.at(0)); reference_image = mitkImage; mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = mitk::convert::GetItkPeakFromPeakImage(mitkImage); dynamic_cast(handler)->SetPeakImage(itkImg); } else if (type == "Tensor" && params->m_Mode == mitk::StreamlineTractographyParameters::MODE::DETERMINISTIC) { handler = new mitk::TrackingHandlerTensor(); MITK_INFO << "loading input tensor images"; std::vector< mitk::Image::Pointer > input_images; for (unsigned int i=0; i(input_files.at(i)); reference_image = mitkImage; mitk::TensorImage::ItkTensorImageType::Pointer itkImg = mitk::convert::GetItkTensorFromTensorImage(mitkImage); dynamic_cast(handler)->AddTensorImage(itkImg.GetPointer()); } if (addImages.at(0).size()>0) dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); } else if (type == "ODF" || type == "ODF-DIPY/FSL" || (type == "Tensor" && params->m_Mode == mitk::StreamlineTractographyParameters::MODE::PROBABILISTIC)) { handler = new mitk::TrackingHandlerOdf(); mitk::OdfImage::ItkOdfImageType::Pointer itkImg = nullptr; if (type == "Tensor") { MITK_INFO << "Converting Tensor to ODF image"; auto input = mitk::IOUtil::Load(input_files.at(0)); reference_image = input; itkImg = mitk::convert::GetItkOdfFromTensorImage(input); dynamic_cast(handler)->SetIsOdfFromTensor(true); } else { std::vector include = {"SH Image", "ODF Image"}; std::vector exclude = {}; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor(include, exclude); auto input = mitk::IOUtil::Load(input_files.at(0), &functor)[0]; reference_image = dynamic_cast(input.GetPointer()); if (dynamic_cast(input.GetPointer())) { MITK_INFO << "Converting SH to ODF image"; mitk::ShImage::Pointer mitkShImage = dynamic_cast(input.GetPointer()); if (type == "ODF-DIPY/FSL") mitkShImage->SetShConvention(mitk::ShImage::SH_CONVENTION::FSL); mitk::Image::Pointer mitkImg = dynamic_cast(mitkShImage.GetPointer()); itkImg = mitk::convert::GetItkOdfFromShImage(mitkImg); } else if (dynamic_cast(input.GetPointer())) { mitk::Image::Pointer mitkImg = dynamic_cast(input.GetPointer()); itkImg = mitk::convert::GetItkOdfFromOdfImage(mitkImg); } else mitkThrow() << ""; } dynamic_cast(handler)->SetOdfImage(itkImg); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); } else { MITK_INFO << "Unknown tractography algorithm (" + type+"). Known types are Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF."; return EXIT_FAILURE; } -// float max_size = 0; -// handler->m_ -// for (int i=0; i<3; ++i) -// if (dynamic_cast(m_InputImageNodes.at(0)->GetData())->GetGeometry()->GetExtentInMM(i)>max_size) -// max_size = dynamic_cast(m_InputImageNodes.at(0)->GetData())->GetGeometry()->GetExtentInMM(i); -// if (params->m_MinTractLengthMm >= max_size) -// { -// MITK_INFO << "Max. image size: " << max_size << "mm"; -// MITK_INFO << "Min. tract length: " << params->m_MinTractLengthMm << "mm"; -// QMessageBox::information(nullptr, "Error", "Minimum tract length exceeds the maximum image extent! Recommended value is about 1/10 of the image extent."); -// StartStopTrackingGui(false); -// return; -// } -// else if (params->m_MinTractLengthMm > max_size/10) -// { -// MITK_INFO << "Max. image size: " << max_size << "mm"; -// MITK_INFO << "Min. tract length: " << params->m_MinTractLengthMm << "mm"; -// MITK_WARN << "Minimum tract length is larger than 1/10 the maximum image extent! Decrease recommended."; -// } + float max_size = 0; + for (int i=0; i<3; ++i) + if (reference_image->GetGeometry()->GetExtentInMM(i)>max_size) + max_size = reference_image->GetGeometry()->GetExtentInMM(i); + if (params->m_MinTractLengthMm >= max_size) + { + MITK_INFO << "Max. image size: " << max_size << "mm"; + MITK_INFO << "Min. tract length: " << params->m_MinTractLengthMm << "mm"; + MITK_ERROR << "Minimum tract length exceeds the maximum image extent! Recommended value is about 1/10 of the image extent."; + return EXIT_FAILURE; + } + else if (params->m_MinTractLengthMm > max_size/10) + { + MITK_INFO << "Max. image size: " << max_size << "mm"; + MITK_INFO << "Min. tract length: " << params->m_MinTractLengthMm << "mm"; + MITK_WARN << "Minimum tract length is larger than 1/10 the maximum image extent! Decrease recommended."; + } MITK_INFO << "Tractography algorithm: " << type; tracker->SetMaskImage(mask); tracker->SetSeedImage(seed); tracker->SetStoppingRegions(stop); tracker->SetTargetRegions(target); tracker->SetExclusionRegions(exclusion); tracker->SetTrackingHandler(handler); if (ext != ".fib" && ext != ".trk") params->m_OutputProbMap = true; tracker->SetParameters(params); tracker->Update(); if (ext == ".fib" || ext == ".trk") { vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); if (params->m_CompressFibers) { float min_sp = 999; auto spacing = handler->GetSpacing(); if (spacing[0] < min_sp) min_sp = spacing[0]; if (spacing[1] < min_sp) min_sp = spacing[1]; if (spacing[2] < min_sp) min_sp = spacing[2]; params->m_Compression = min_sp/10; outFib->Compress(params->m_Compression); } outFib->SetTrackVisHeader(reference_image->GetGeometry()); mitk::IOUtil::Save(outFib, outFile); } else { TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (ext != ".nii" && ext != ".nii.gz" && ext != ".nrrd") outFile += ".nii.gz"; mitk::IOUtil::Save(img, outFile); } delete handler; return EXIT_SUCCESS; } diff --git a/Modules/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 6bd9546..a29649f 100644 --- a/Modules/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,999 +1,999 @@ /*=================================================================== 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 "itkStreamlineTrackingFilter.h" #include #include #include #include "itkPointShell.h" #include #include #include #include #include #include #include #include namespace itk { StreamlineTrackingFilter ::StreamlineTrackingFilter() : m_PauseTracking(false) , m_AbortTracking(false) , m_BuildFibersFinished(false) , m_BuildFibersReady(0) , m_FiberPolyData(nullptr) , m_Points(nullptr) , m_Cells(nullptr) , m_StoppingRegions(nullptr) , m_TargetRegions(nullptr) , m_SeedImage(nullptr) , m_MaskImage(nullptr) , m_ExclusionRegions(nullptr) , m_OutputProbabilityMap(nullptr) , m_Verbose(true) , m_DemoMode(false) , m_CurrentTracts(0) , m_Progress(0) , m_StopTracking(false) , m_TrackingPriorHandler(nullptr) { this->SetNumberOfRequiredInputs(0); } std::string StreamlineTrackingFilter::GetStatusText() { std::string status = "Seedpoints processed: " + boost::lexical_cast(m_Progress) + "/" + boost::lexical_cast(m_SeedPoints.size()); if (m_SeedPoints.size()>0) status += " (" + boost::lexical_cast(100*m_Progress/m_SeedPoints.size()) + "%)"; if (m_Parameters->m_MaxNumFibers>0) status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts) + "/" + boost::lexical_cast(m_Parameters->m_MaxNumFibers); else status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts); return status; } void StreamlineTrackingFilter::BeforeTracking() { m_StopTracking = false; m_TrackingHandler->SetParameters(m_Parameters); m_TrackingHandler->InitForTracking(); m_FiberPolyData = PolyDataType::New(); m_Points = vtkSmartPointer< vtkPoints >::New(); m_Cells = vtkSmartPointer< vtkCellArray >::New(); if (m_TrackingPriorHandler!=nullptr) { m_TrackingPriorHandler->InitForTracking(); } m_PolyDataContainer.clear(); for (unsigned int i=0; iGetNumberOfThreads(); i++) { PolyDataType poly = PolyDataType::New(); m_PolyDataContainer.push_back(poly); } auto imageSpacing = m_TrackingHandler->GetSpacing(); if (m_Parameters->m_OutputProbMap) { m_OutputProbabilityMap = ItkDoubleImgType::New(); m_OutputProbabilityMap->SetSpacing(imageSpacing); m_OutputProbabilityMap->SetOrigin(m_TrackingHandler->GetOrigin()); m_OutputProbabilityMap->SetDirection(m_TrackingHandler->GetDirection()); m_OutputProbabilityMap->SetRegions(m_TrackingHandler->GetLargestPossibleRegion()); m_OutputProbabilityMap->Allocate(); m_OutputProbabilityMap->FillBuffer(0); } m_MaskInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_StopInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_SeedInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_TargetInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_ExclusionInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); if (m_StoppingRegions.IsNull()) { m_StoppingRegions = ItkFloatImgType::New(); m_StoppingRegions->SetSpacing( imageSpacing ); m_StoppingRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); m_StoppingRegions->SetDirection( m_TrackingHandler->GetDirection() ); m_StoppingRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_StoppingRegions->Allocate(); m_StoppingRegions->FillBuffer(0); } else std::cout << "StreamlineTracking - Using stopping region image" << std::endl; m_StopInterpolator->SetInputImage(m_StoppingRegions); if (m_ExclusionRegions.IsNotNull()) { std::cout << "StreamlineTracking - Using exclusion region image" << std::endl; m_ExclusionInterpolator->SetInputImage(m_ExclusionRegions); } if (m_TargetRegions.IsNull()) { m_TargetImageSet = false; m_TargetRegions = ItkFloatImgType::New(); m_TargetRegions->SetSpacing( imageSpacing ); m_TargetRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); m_TargetRegions->SetDirection( m_TrackingHandler->GetDirection() ); m_TargetRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_TargetRegions->Allocate(); m_TargetRegions->FillBuffer(1); } else { m_TargetImageSet = true; m_TargetInterpolator->SetInputImage(m_TargetRegions); std::cout << "StreamlineTracking - Using target region image" << std::endl; } if (m_SeedImage.IsNull()) { m_SeedImageSet = false; m_SeedImage = ItkFloatImgType::New(); m_SeedImage->SetSpacing( imageSpacing ); m_SeedImage->SetOrigin( m_TrackingHandler->GetOrigin() ); m_SeedImage->SetDirection( m_TrackingHandler->GetDirection() ); m_SeedImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_SeedImage->Allocate(); m_SeedImage->FillBuffer(1); } else { m_SeedImageSet = true; std::cout << "StreamlineTracking - Using seed image" << std::endl; } m_SeedInterpolator->SetInputImage(m_SeedImage); if (m_MaskImage.IsNull()) { // initialize mask image m_MaskImage = ItkFloatImgType::New(); m_MaskImage->SetSpacing( imageSpacing ); m_MaskImage->SetOrigin( m_TrackingHandler->GetOrigin() ); m_MaskImage->SetDirection( m_TrackingHandler->GetDirection() ); m_MaskImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } else std::cout << "StreamlineTracking - Using mask image" << std::endl; m_MaskInterpolator->SetInputImage(m_MaskImage); // Autosettings for endpoint constraints if (m_Parameters->m_EpConstraints==EndpointConstraints::NONE && m_TargetImageSet && m_SeedImageSet) { MITK_INFO << "No endpoint constraint chosen but seed and target image set --> setting constraint to EPS_IN_SEED_AND_TARGET"; m_Parameters->m_EpConstraints = EndpointConstraints::EPS_IN_SEED_AND_TARGET; } else if (m_Parameters->m_EpConstraints==EndpointConstraints::NONE && m_TargetImageSet) { MITK_INFO << "No endpoint constraint chosen but target image set --> setting constraint to EPS_IN_TARGET"; m_Parameters->m_EpConstraints = EndpointConstraints::EPS_IN_TARGET; } // Check if endpoint constraints are valid FiberType test_fib; itk::Point p; p.Fill(0); test_fib.push_back(p); test_fib.push_back(p); IsValidFiber(&test_fib); if (m_SeedPoints.empty()) GetSeedPointsFromSeedImage(); m_BuildFibersReady = 0; m_BuildFibersFinished = false; m_Tractogram.clear(); m_SamplingPointset = mitk::PointSet::New(); m_AlternativePointset = mitk::PointSet::New(); m_StopVotePointset = mitk::PointSet::New(); m_StartTime = std::chrono::system_clock::now(); if (m_DemoMode) omp_set_num_threads(1); if (m_Parameters->m_Mode==mitk::TrackingDataHandler::MODE::DETERMINISTIC) std::cout << "StreamlineTracking - Mode: deterministic" << std::endl; else if(m_Parameters->m_Mode==mitk::TrackingDataHandler::MODE::PROBABILISTIC) { std::cout << "StreamlineTracking - Mode: probabilistic" << std::endl; std::cout << "StreamlineTracking - Trials per seed: " << m_Parameters->m_TrialsPerSeed << std::endl; } else std::cout << "StreamlineTracking - Mode: ???" << std::endl; if (m_Parameters->m_EpConstraints==EndpointConstraints::NONE) std::cout << "StreamlineTracking - Endpoint constraint: NONE" << std::endl; else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_TARGET" << std::endl; else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_TARGET_LABELDIFF) std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_TARGET_LABELDIFF" << std::endl; else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_SEED_AND_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_SEED_AND_TARGET" << std::endl; else if (m_Parameters->m_EpConstraints==EndpointConstraints::MIN_ONE_EP_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: MIN_ONE_EP_IN_TARGET" << std::endl; else if (m_Parameters->m_EpConstraints==EndpointConstraints::ONE_EP_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: ONE_EP_IN_TARGET" << std::endl; else if (m_Parameters->m_EpConstraints==EndpointConstraints::NO_EP_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: NO_EP_IN_TARGET" << std::endl; - std::cout << "StreamlineTracking - Angular threshold: " << m_Parameters->GetAngularThresholdDot() << "°" << std::endl; + std::cout << "StreamlineTracking - Angular threshold: " << m_Parameters->GetAngularThresholdDeg() << "°" << std::endl; std::cout << "StreamlineTracking - Stepsize: " << m_Parameters->GetStepSizeMm() << "mm (" << m_Parameters->GetStepSizeMm()/m_Parameters->GetMinVoxelSizeMm() << "*vox)" << std::endl; std::cout << "StreamlineTracking - Seeds per voxel: " << m_Parameters->m_SeedsPerVoxel << std::endl; std::cout << "StreamlineTracking - Max. tract length: " << m_Parameters->m_MaxTractLengthMm << "mm" << std::endl; std::cout << "StreamlineTracking - Min. tract length: " << m_Parameters->m_MinTractLengthMm << "mm" << std::endl; std::cout << "StreamlineTracking - Max. num. tracts: " << m_Parameters->m_MaxNumFibers << std::endl; std::cout << "StreamlineTracking - Loop check: " << m_Parameters->GetLoopCheckDeg() << "°" << std::endl; std::cout << "StreamlineTracking - Num. neighborhood samples: " << m_Parameters->m_NumSamples << std::endl; std::cout << "StreamlineTracking - Max. sampling distance: " << m_Parameters->GetSamplingDistanceMm() << "mm (" << m_Parameters->GetSamplingDistanceMm()/m_Parameters->GetMinVoxelSizeMm() << "*vox)" << std::endl; std::cout << "StreamlineTracking - Deflection modifier: " << m_Parameters->m_DeflectionMod << std::endl; std::cout << "StreamlineTracking - Use stop votes: " << m_Parameters->m_StopVotes << std::endl; std::cout << "StreamlineTracking - Only frontal samples: " << m_Parameters->m_OnlyForwardSamples << std::endl; if (m_TrackingPriorHandler!=nullptr) std::cout << "StreamlineTracking - Using directional prior for tractography (w=" << m_Parameters->m_Weight << ")" << std::endl; if (m_DemoMode) { std::cout << "StreamlineTracking - Running in demo mode"; std::cout << "StreamlineTracking - Starting streamline tracking using 1 thread" << std::endl; } else std::cout << "StreamlineTracking - Starting streamline tracking using " << omp_get_max_threads() << " threads" << std::endl; } void StreamlineTrackingFilter::CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir) { pos[0] += dir[0]*m_Parameters->GetStepSizeMm(); pos[1] += dir[1]*m_Parameters->GetStepSizeMm(); pos[2] += dir[2]*m_Parameters->GetStepSizeMm(); } std::vector< vnl_vector_fixed > StreamlineTrackingFilter::CreateDirections(unsigned int NPoints) { std::vector< vnl_vector_fixed > pointshell; if (NPoints<2) return pointshell; std::vector< double > theta; theta.resize(NPoints); std::vector< double > phi; phi.resize(NPoints); auto C = sqrt(4*itk::Math::pi); phi[0] = 0.0; phi[NPoints-1] = 0.0; for(unsigned int i=0; i0 && i d; d[0] = static_cast(cos(theta[i]) * cos(phi[i])); d[1] = static_cast(cos(theta[i]) * sin(phi[i])); d[2] = static_cast(sin(theta[i])); pointshell.push_back(d); } return pointshell; } vnl_vector_fixed StreamlineTrackingFilter::GetNewDirection(const itk::Point &pos, std::deque >& olddirs, itk::Index<3> &oldIndex) { if (m_DemoMode) { m_SamplingPointset->Clear(); m_AlternativePointset->Clear(); m_StopVotePointset->Clear(); } vnl_vector_fixed direction; direction.fill(0); if (mitk::imv::IsInsideMask(pos, m_Parameters->m_InterpolateRoiImages, m_MaskInterpolator) && !mitk::imv::IsInsideMask(pos, m_Parameters->m_InterpolateRoiImages, m_StopInterpolator)) direction = m_TrackingHandler->ProposeDirection(pos, olddirs, oldIndex); // get direction proposal at current streamline position else return direction; int stop_votes = 0; int possible_stop_votes = 0; if (!olddirs.empty()) { vnl_vector_fixed olddir = olddirs.back(); std::vector< vnl_vector_fixed > probeVecs = CreateDirections(m_Parameters->m_NumSamples); itk::Point sample_pos; unsigned int alternatives = 1; for (unsigned int i=0; i d; bool is_stop_voter = false; if (!m_Parameters->m_FixRandomSeed && m_Parameters->m_RandomSampling) { d[0] = static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); d[1] = static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); d[2] = static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); d.normalize(); d *= static_cast(m_TrackingHandler->GetRandDouble(0, static_cast(m_Parameters->GetSamplingDistanceMm()))); } else { d = probeVecs.at(i); float dot = dot_product(d, olddir); if (m_Parameters->m_StopVotes && dot>0.7f) { is_stop_voter = true; possible_stop_votes++; } else if (m_Parameters->m_OnlyForwardSamples && dot<0) continue; d *= m_Parameters->GetSamplingDistanceMm(); } sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (mitk::imv::IsInsideMask(sample_pos, m_Parameters->m_InterpolateRoiImages, m_MaskInterpolator)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>static_cast(mitk::eps)) { direction += tempDir; if(m_DemoMode) m_SamplingPointset->InsertPoint(i, sample_pos); } else if (m_Parameters->m_AvoidStop && olddir.magnitude()>0.5f) // out of white matter { if (is_stop_voter) stop_votes++; if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); float dot = dot_product(d, olddir); if (dot >= 0.0f) // in front of plane defined by pos and olddir d = -d + 2*dot*olddir; // reflect else d = -d; // invert // look a bit further into the other direction sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; alternatives++; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (mitk::imv::IsInsideMask(sample_pos, m_Parameters->m_InterpolateRoiImages, m_MaskInterpolator)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>static_cast(mitk::eps)) // are we back in the white matter? { direction += d * m_Parameters->m_DeflectionMod; // go into the direction of the white matter direction += tempDir; // go into the direction of the white matter direction at this location if(m_DemoMode) m_AlternativePointset->InsertPoint(alternatives, sample_pos); } else { if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); } } else { if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); if (is_stop_voter) stop_votes++; } } } bool valid = false; if (direction.magnitude()>0.001f && (possible_stop_votes==0 || static_cast(stop_votes)/possible_stop_votes<0.5f) ) { direction.normalize(); valid = true; } else direction.fill(0); if (m_TrackingPriorHandler!=nullptr && (m_Parameters->m_NewDirectionsFromPrior || valid)) { vnl_vector_fixed prior = m_TrackingPriorHandler->ProposeDirection(pos, olddirs, oldIndex); if (prior.magnitude()>0.001f) { prior.normalize(); if (dot_product(prior,direction)<0) prior *= -1; direction = (1.0f-m_Parameters->m_Weight) * direction + m_Parameters->m_Weight * prior; direction.normalize(); } else if (m_Parameters->m_RestrictToPrior) direction.fill(0.0); } return direction; } float StreamlineTrackingFilter::FollowStreamline(itk::Point pos, vnl_vector_fixed dir, FiberType* fib, DirectionContainer* container, float tractLength, bool front, bool &exclude) { vnl_vector_fixed zero_dir; zero_dir.fill(0.0); std::deque< vnl_vector_fixed > last_dirs; for (unsigned int i=0; im_NumPreviousDirections-1; i++) last_dirs.push_back(zero_dir); for (int step=0; step< 5000; step++) { itk::Index<3> oldIndex; m_TrackingHandler->WorldToIndex(pos, oldIndex); // get new position CalculateNewPosition(pos, dir); if (m_ExclusionRegions.IsNotNull() && mitk::imv::IsInsideMask(pos, m_Parameters->m_InterpolateRoiImages, m_ExclusionInterpolator)) { exclude = true; return tractLength; } if (m_AbortTracking) return tractLength; // if yes, add new point to streamline dir.normalize(); if (front) { fib->push_front(pos); container->push_front(dir); } else { fib->push_back(pos); container->push_back(dir); } tractLength += m_Parameters->GetStepSizeMm(); if (m_Parameters->GetLoopCheckDeg()>=0 && CheckCurvature(container, front)>m_Parameters->GetLoopCheckDeg()) return tractLength; if (tractLength>m_Parameters->m_MaxTractLengthMm) return tractLength; if (m_DemoMode && !m_Parameters->m_OutputProbMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? { #pragma omp critical { m_BuildFibersReady++; m_Tractogram.push_back(*fib); BuildFibers(true); m_Stop = true; while (m_Stop){ } } } last_dirs.push_back(dir); if (last_dirs.size()>m_Parameters->m_NumPreviousDirections) last_dirs.pop_front(); dir = GetNewDirection(pos, last_dirs, oldIndex); while (m_PauseTracking){} if (dir.magnitude()<0.0001f) return tractLength; } return tractLength; } float StreamlineTrackingFilter::CheckCurvature(DirectionContainer* fib, bool front) { if (fib->size()<8) return 0; float m_Distance = std::max(m_Parameters->GetMinVoxelSizeMm()*4, m_Parameters->GetStepSizeMm()*8); float dist = 0; std::vector< vnl_vector_fixed< float, 3 > > vectors; vnl_vector_fixed< float, 3 > meanV; meanV.fill(0); float dev = 0; if (front) { int c = 0; while(dist(fib->size())-1) { dist += m_Parameters->GetStepSizeMm(); vnl_vector_fixed< float, 3 > v = fib->at(static_cast(c)); if (dot_product(v,meanV)<0) v = -v; vectors.push_back(v); meanV += v; c++; } } else { int c = static_cast(fib->size())-1; while(dist=0) { dist += m_Parameters->GetStepSizeMm(); vnl_vector_fixed< float, 3 > v = fib->at(static_cast(c)); if (dot_product(v,meanV)<0) v = -v; vectors.push_back(v); meanV += v; c--; } } meanV.normalize(); for (unsigned int c=0; c1.0f) angle = 1.0; dev += acos(angle)*180.0f/static_cast(itk::Math::pi); } if (vectors.size()>0) dev /= vectors.size(); return dev; } std::shared_ptr StreamlineTrackingFilter::GetParameters() const { return m_Parameters; } void StreamlineTrackingFilter::SetParameters(std::shared_ptr< mitk::StreamlineTractographyParameters > Parameters) { m_Parameters = Parameters; } void StreamlineTrackingFilter::SetTrackingPriorHandler(mitk::TrackingDataHandler *TrackingPriorHandler) { m_TrackingPriorHandler = TrackingPriorHandler; } void StreamlineTrackingFilter::GetSeedPointsFromSeedImage() { MITK_INFO << "StreamlineTracking - Calculating seed points."; m_SeedPoints.clear(); typedef ImageRegionConstIterator< ItkFloatImgType > MaskIteratorType; MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); sit.GoToBegin(); while (!sit.IsAtEnd()) { if (sit.Value()>0) { ItkFloatImgType::IndexType index = sit.GetIndex(); itk::ContinuousIndex start; start[0] = index[0]; start[1] = index[1]; start[2] = index[2]; itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); if ( mitk::imv::IsInsideMask(worldPos, m_Parameters->m_InterpolateRoiImages, m_MaskInterpolator) ) { m_SeedPoints.push_back(worldPos); for (unsigned int s = 1; s < m_Parameters->m_SeedsPerVoxel; s++) { start[0] = index[0] + static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); start[1] = index[1] + static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); start[2] = index[2] + static_cast(m_TrackingHandler->GetRandDouble(-0.5, 0.5)); itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); m_SeedPoints.push_back(worldPos); } } } ++sit; } if (m_SeedPoints.empty()) mitkThrow() << "No valid seed point in seed image! Is your seed image registered with the image you are tracking on?"; } void StreamlineTrackingFilter::GenerateData() { this->BeforeTracking(); if (!m_Parameters->m_FixRandomSeed) std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); m_CurrentTracts = 0; int num_seeds = static_cast(m_SeedPoints.size()); itk::Index<3> zeroIndex; zeroIndex.Fill(0); m_Progress = 0; int i = 0; int print_interval = num_seeds/100; if (print_interval<100) m_Verbose=false; unsigned int trials_per_seed = 1; if(m_Parameters->m_Mode==mitk::TrackingDataHandler::MODE::PROBABILISTIC) trials_per_seed = m_Parameters->m_TrialsPerSeed; #pragma omp parallel while (i=num_seeds || m_StopTracking) continue; else if (m_Verbose && i%print_interval==0) #pragma omp critical { m_Progress += static_cast(print_interval); std::cout << " \r"; if (m_Parameters->m_MaxNumFibers>0) std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << "/" << m_Parameters->m_MaxNumFibers << '\r'; else std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << '\r'; cout.flush(); } const itk::Point worldPos = m_SeedPoints.at(static_cast(temp_i)); for (unsigned int trials=0; trials dir; dir.fill(0.0); std::deque< vnl_vector_fixed > olddirs; dir = GetNewDirection(worldPos, olddirs, zeroIndex) * 0.5f; bool exclude = false; if (m_ExclusionRegions.IsNotNull() && mitk::imv::IsInsideMask(worldPos, m_Parameters->m_InterpolateRoiImages, m_ExclusionInterpolator)) exclude = true; bool success = false; if (dir.magnitude()>0.0001f && !exclude) { // forward tracking tractLength = FollowStreamline(worldPos, dir, &fib, &direction_container, 0, false, exclude); fib.push_front(worldPos); // backward tracking if (!exclude) tractLength = FollowStreamline(worldPos, -dir, &fib, &direction_container, tractLength, true, exclude); counter = fib.size(); if (tractLength>=m_Parameters->m_MinTractLengthMm && counter>=2 && !exclude) { #pragma omp critical if ( IsValidFiber(&fib) ) { if (!m_StopTracking) { if (!m_Parameters->m_OutputProbMap) m_Tractogram.push_back(fib); else FiberToProbmap(&fib); m_CurrentTracts++; success = true; } if (m_Parameters->m_MaxNumFibers > 0 && m_CurrentTracts>=static_cast(m_Parameters->m_MaxNumFibers)) { if (!m_StopTracking) { std::cout << " \r"; MITK_INFO << "Reconstructed maximum number of tracts (" << m_CurrentTracts << "). Stopping tractography."; } m_StopTracking = true; } } } } if (success || m_Parameters->m_Mode!=MODE::PROBABILISTIC) break; // we only try one seed point multiple times if we use a probabilistic tracker and have not found a valid streamline yet }// trials per seed }// seed points this->AfterTracking(); } bool StreamlineTrackingFilter::IsValidFiber(FiberType* fib) { if (m_Parameters->m_EpConstraints==EndpointConstraints::NONE) { return true; } else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) && mitk::imv::IsInsideMask(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint EPS_IN_TARGET chosen!"; } else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_TARGET_LABELDIFF) { if (m_TargetImageSet) { float v1 = mitk::imv::GetImageValue(fib->front(), false, m_TargetInterpolator); float v2 = mitk::imv::GetImageValue(fib->back(), false, m_TargetInterpolator); if ( v1>0.0f && v2>0.0f && v1!=v2 ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint EPS_IN_TARGET_LABELDIFF chosen!"; } else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_SEED_AND_TARGET) { if (m_TargetImageSet && m_SeedImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_Parameters->m_InterpolateRoiImages, m_SeedInterpolator) && mitk::imv::IsInsideMask(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) ) return true; if ( mitk::imv::IsInsideMask(fib->back(), m_Parameters->m_InterpolateRoiImages, m_SeedInterpolator) && mitk::imv::IsInsideMask(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target or seed image set but endpoint constraint EPS_IN_SEED_AND_TARGET chosen!"; } else if (m_Parameters->m_EpConstraints==EndpointConstraints::MIN_ONE_EP_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) || mitk::imv::IsInsideMask(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint MIN_ONE_EP_IN_TARGET chosen!"; } else if (m_Parameters->m_EpConstraints==EndpointConstraints::ONE_EP_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) && !mitk::imv::IsInsideMask(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) ) return true; if ( !mitk::imv::IsInsideMask(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) && mitk::imv::IsInsideMask(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) ) return true; return false; } else mitkThrow() << "No target image set but endpoint constraint ONE_EP_IN_TARGET chosen!"; } else if (m_Parameters->m_EpConstraints==EndpointConstraints::NO_EP_IN_TARGET) { if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) || mitk::imv::IsInsideMask(fib->back(), m_Parameters->m_InterpolateRoiImages, m_TargetInterpolator) ) return false; return true; } else mitkThrow() << "No target image set but endpoint constraint NO_EP_IN_TARGET chosen!"; } return true; } void StreamlineTrackingFilter::FiberToProbmap(FiberType* fib) { ItkDoubleImgType::IndexType last_idx; last_idx.Fill(0); for (auto p : *fib) { ItkDoubleImgType::IndexType idx; m_OutputProbabilityMap->TransformPhysicalPointToIndex(p, idx); if (idx != last_idx) { if (m_OutputProbabilityMap->GetLargestPossibleRegion().IsInside(idx)) m_OutputProbabilityMap->SetPixel(idx, m_OutputProbabilityMap->GetPixel(idx)+1); last_idx = idx; } } } void StreamlineTrackingFilter::BuildFibers(bool check) { if (m_BuildFibersReady::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); for (unsigned int i=0; i container = vtkSmartPointer::New(); FiberType fib = m_Tractogram.at(i); for (FiberType::iterator it = fib.begin(); it!=fib.end(); ++it) { vtkIdType id = vNewPoints->InsertNextPoint((*it).GetDataPointer()); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } if (check) for (int i=0; iSetPoints(vNewPoints); m_FiberPolyData->SetLines(vNewLines); m_BuildFibersFinished = true; } void StreamlineTrackingFilter::AfterTracking() { if (m_Verbose) std::cout << " \r"; if (!m_Parameters->m_OutputProbMap) { MITK_INFO << "Reconstructed " << m_Tractogram.size() << " fibers."; MITK_INFO << "Generating polydata "; BuildFibers(false); } else { itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::Pointer filter = itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); filter->SetInput(m_OutputProbabilityMap); filter->SetOutputMaximum(1.0); filter->SetOutputMinimum(0.0); filter->Update(); m_OutputProbabilityMap = filter->GetOutput(); } MITK_INFO << "done"; 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); std::chrono::seconds ss = std::chrono::duration_cast(m_EndTime - m_StartTime); mm %= 60; ss %= 60; MITK_INFO << "Tracking took " << hh.count() << "h, " << mm.count() << "m and " << ss.count() << "s"; m_SeedPoints.clear(); } void StreamlineTrackingFilter::SetDicomProperties(mitk::FiberBundle::Pointer fib) { std::string model_code_value = "-"; std::string model_code_meaning = "-"; std::string algo_code_value = "-"; std::string algo_code_meaning = "-"; if ( m_Parameters->m_Mode==MODE::DETERMINISTIC && dynamic_cast(m_TrackingHandler)) { algo_code_value = "sup181_ee01"; algo_code_meaning = "Deterministic"; if (m_Parameters->m_F > 0.99 && m_Parameters->m_G < 0.01) { if (m_Parameters->m_InterpolateTractographyData) { algo_code_value = "sup181_ee08"; algo_code_meaning = "Euler"; } else { algo_code_value = "sup181_ee04"; algo_code_meaning = "FACT"; } } else if (m_Parameters->m_G > 0.99 && m_Parameters->m_F < 0.01) { algo_code_value = "sup181_ee06"; algo_code_meaning = "TEND"; } } else if (m_Parameters->m_Mode==MODE::DETERMINISTIC) { algo_code_value = "sup181_ee01"; algo_code_meaning = "Deterministic"; } else if (m_Parameters->m_Mode==MODE::PROBABILISTIC) { algo_code_value = "sup181_ee02"; algo_code_meaning = "Probabilistic"; } if (dynamic_cast(m_TrackingHandler) || (dynamic_cast(m_TrackingHandler) && dynamic_cast(m_TrackingHandler)->GetIsOdfFromTensor() ) ) { if ( dynamic_cast(m_TrackingHandler) && dynamic_cast(m_TrackingHandler)->GetNumTensorImages()>1 ) { model_code_value = "sup181_bb02"; model_code_meaning = "Multi Tensor"; } else { model_code_value = "sup181_bb01"; model_code_meaning = "Single Tensor"; } } else if (dynamic_cast*>(m_TrackingHandler) || dynamic_cast*>(m_TrackingHandler)) { model_code_value = "sup181_bb03"; model_code_meaning = "Model Free"; } else if (dynamic_cast(m_TrackingHandler)) { model_code_value = "-"; model_code_meaning = "ODF"; } else if (dynamic_cast(m_TrackingHandler)) { model_code_value = "-"; model_code_meaning = "Peaks"; } fib->SetProperty("DICOM.anatomy.value", mitk::StringProperty::New("T-A0095")); fib->SetProperty("DICOM.anatomy.meaning", mitk::StringProperty::New("White matter of brain and spinal cord")); fib->SetProperty("DICOM.algo_family_code.value", mitk::StringProperty::New(algo_code_value)); fib->SetProperty("DICOM.algo_family_code.meaning", mitk::StringProperty::New(algo_code_meaning)); fib->SetProperty("DICOM.model_code.value", mitk::StringProperty::New(model_code_value)); fib->SetProperty("DICOM.model_code.meaning", mitk::StringProperty::New(model_code_meaning)); } }