diff --git a/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/StreamlineTractography.cpp b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/StreamlineTractography.cpp index 6b24985008..e486ae680f 100755 --- a/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionImaging/DiffusionCmdApps/Tractography/StreamlineTractography.cpp @@ -1,592 +1,563 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #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[]) { mitkCommandLineParser 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", mitkCommandLineParser::StringList, "Input:", "input image (multiple possible for 'DetTensor' algorithm)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument("", "o", mitkCommandLineParser::String, "Output:", "output fiberbundle/probability map", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument("algorithm", "", mitkCommandLineParser::String, "Algorithm:", "which algorithm to use (DetPeaks; ProbPeaks; DetTensor; ProbTensor; DetODF; ProbODF; DetRF; ProbRF)", us::Any(), false); parser.endGroup(); parser.beginGroup("2. Seeding:"); parser.addArgument("seeds", "", mitkCommandLineParser::Int, "Seeds per voxel:", "number of seed points per voxel", 1); parser.addArgument("seed_image", "", mitkCommandLineParser::String, "Seed image:", "mask image defining seed voxels", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("trials_per_seed", "", mitkCommandLineParser::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", "", mitkCommandLineParser::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", "", mitkCommandLineParser::String, "Mask image:", "streamlines leaving the mask will stop immediately", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("stop_image", "", mitkCommandLineParser::String, "Stop ROI image:", "streamlines entering the mask will stop immediately", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("exclusion_image", "", mitkCommandLineParser::String, "Exclusion ROI image:", "streamlines entering the mask will be discarded", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("ep_constraint", "", mitkCommandLineParser::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", "", mitkCommandLineParser::String, "Target ROI image:", "effact depends on the chosen endpoint constraint (option ep_constraint)", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.endGroup(); parser.beginGroup("4. Streamline integration parameters:"); parser.addArgument("sharpen_odfs", "", mitkCommandLineParser::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", "", mitkCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 0.1); parser.addArgument("odf_cutoff", "", mitkCommandLineParser::Float, "ODF Cutoff:", "threshold on the ODF magnitude. this is useful in case of CSD fODF tractography.", 0.0); parser.addArgument("step_size", "", mitkCommandLineParser::Float, "Step size:", "step size (in voxels)", 0.5); parser.addArgument("min_tract_length", "", mitkCommandLineParser::Float, "Min. tract length:", "minimum fiber length (in mm)", 20); parser.addArgument("angular_threshold", "", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold between two successive steps, (default: 90° * step_size, minimum 15°)"); parser.addArgument("loop_check", "", mitkCommandLineParser::Float, "Check for loops:", "threshold on angular stdev over the last 4 voxel lengths"); parser.endGroup(); parser.beginGroup("5. Tractography prior:"); parser.addArgument("prior_image", "", mitkCommandLineParser::String, "Peak prior:", "tractography prior in thr for of a peak image", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("prior_weight", "", mitkCommandLineParser::Float, "Prior weight", "weighting factor between prior and data.", 0.5); - parser.addArgument("restrict_to_prior", "", mitkCommandLineParser::Bool, "Restrict to prior:", "restrict tractography to regions where the prior is valid."); - parser.addArgument("new_directions_from_prior", "", mitkCommandLineParser::Bool, "New directios from prior:", "the prior can create directions where there are none in the data."); + parser.addArgument("dont_restrict_to_prior", "", mitkCommandLineParser::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", "", mitkCommandLineParser::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", "", mitkCommandLineParser::Bool, "Prior Flip X:", "multiply x-coordinate of prior direction by -1"); parser.addArgument("prior_flip_y", "", mitkCommandLineParser::Bool, "Prior Flip Y:", "multiply y-coordinate of prior direction by -1"); parser.addArgument("prior_flip_z", "", mitkCommandLineParser::Bool, "Prior Flip Z:", "multiply z-coordinate of prior direction by -1"); parser.endGroup(); parser.beginGroup("6. Neighborhood sampling:"); parser.addArgument("num_samples", "", mitkCommandLineParser::Int, "Num. neighborhood samples:", "number of neighborhood samples that are use to determine the next progression direction", 0); parser.addArgument("sampling_distance", "", mitkCommandLineParser::Float, "Sampling distance:", "distance of neighborhood sampling points (in voxels)", 0.25); parser.addArgument("use_stop_votes", "", mitkCommandLineParser::Bool, "Use stop votes:", "use stop votes"); parser.addArgument("use_only_forward_samples", "", mitkCommandLineParser::Bool, "Use only forward samples:", "use only forward samples"); parser.endGroup(); parser.beginGroup("7. Tensor tractography specific:"); parser.addArgument("tend_f", "", mitkCommandLineParser::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", "", mitkCommandLineParser::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", "", mitkCommandLineParser::String, "Forest:", "input random forest (HDF5 file)", us::Any(), true, false, false, mitkCommandLineParser::Input); parser.addArgument("use_sh_features", "", mitkCommandLineParser::Bool, "Use SH features:", "use SH features"); parser.endGroup(); parser.beginGroup("9. Additional input:"); parser.addArgument("additional_images", "", mitkCommandLineParser::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, mitkCommandLineParser::Input); parser.endGroup(); parser.beginGroup("10. Misc:"); parser.addArgument("flip_x", "", mitkCommandLineParser::Bool, "Flip X:", "multiply x-coordinate of direction proposal by -1"); parser.addArgument("flip_y", "", mitkCommandLineParser::Bool, "Flip Y:", "multiply y-coordinate of direction proposal by -1"); parser.addArgument("flip_z", "", mitkCommandLineParser::Bool, "Flip Z:", "multiply z-coordinate of direction proposal by -1"); parser.addArgument("no_data_interpolation", "", mitkCommandLineParser::Bool, "Don't interpolate input data:", "don't interpolate input image values"); parser.addArgument("no_mask_interpolation", "", mitkCommandLineParser::Bool, "Don't interpolate masks:", "don't interpolate mask image values"); parser.addArgument("compress", "", mitkCommandLineParser::Float, "Compress:", "compress output fibers using the given error threshold (in mm)"); parser.addArgument("fix_seed", "", mitkCommandLineParser::Bool, "Fix Random Seed:", "always use the same random numbers"); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["i"]); std::string outFile = us::any_cast(parsedArgs["o"]); std::string algorithm = us::any_cast(parsedArgs["algorithm"]); + std::shared_ptr< mitk::StreamlineTractographyParameters > params = std::make_shared(); + std::string prior_image = ""; if (parsedArgs.count("prior_image")) prior_image = us::any_cast(parsedArgs["prior_image"]); - float prior_weight = 0.5; if (parsedArgs.count("prior_weight")) - prior_weight = us::any_cast(parsedArgs["prior_weight"]); + params->m_Weight = us::any_cast(parsedArgs["prior_weight"]); - bool fix_seed = false; if (parsedArgs.count("fix_seed")) - fix_seed = us::any_cast(parsedArgs["fix_seed"]); + params->m_FixRandomSeed = us::any_cast(parsedArgs["fix_seed"]); - bool restrict_to_prior = false; - if (parsedArgs.count("restrict_to_prior")) - restrict_to_prior = us::any_cast(parsedArgs["restrict_to_prior"]); + params->m_RestrictToPrior = true; + if (parsedArgs.count("dont_restrict_to_prior")) + params->m_RestrictToPrior = !us::any_cast(parsedArgs["dont_restrict_to_prior"]); - bool new_directions_from_prior = false; - if (parsedArgs.count("new_directions_from_prior")) - new_directions_from_prior = us::any_cast(parsedArgs["new_directions_from_prior"]); + params->m_NewDirectionsFromPrior = true; + if (parsedArgs.count("no_new_directions_from_prior")) + params->m_NewDirectionsFromPrior = !us::any_cast(parsedArgs["no_new_directions_from_prior"]); - bool sharpen_odfs = false; + params->m_SharpenOdfs = false; if (parsedArgs.count("sharpen_odfs")) - sharpen_odfs = us::any_cast(parsedArgs["sharpen_odfs"]); + params->m_SharpenOdfs = us::any_cast(parsedArgs["sharpen_odfs"]); - bool interpolate = true; + params->m_InterpolateTractographyData = true; if (parsedArgs.count("no_data_interpolation")) - interpolate = !us::any_cast(parsedArgs["no_data_interpolation"]); + params->m_InterpolateTractographyData = !us::any_cast(parsedArgs["no_data_interpolation"]); - bool mask_interpolation = true; + params->m_InterpolateRoiImages = true; if (parsedArgs.count("no_mask_interpolation")) - interpolate = !us::any_cast(parsedArgs["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"]); - bool use_stop_votes = false; + params->m_StopVotes = false; if (parsedArgs.count("use_stop_votes")) - use_stop_votes = us::any_cast(parsedArgs["use_stop_votes"]); + params->m_StopVotes = us::any_cast(parsedArgs["use_stop_votes"]); - bool use_only_forward_samples = false; + params->m_OnlyForwardSamples = false; if (parsedArgs.count("use_only_forward_samples")) - use_only_forward_samples = us::any_cast(parsedArgs["use_only_forward_samples"]); + params->m_OnlyForwardSamples = us::any_cast(parsedArgs["use_only_forward_samples"]); - bool flip_x = false; + params->m_FlipX = false; if (parsedArgs.count("flip_x")) - flip_x = us::any_cast(parsedArgs["flip_x"]); - bool flip_y = false; + params->m_FlipX = us::any_cast(parsedArgs["flip_x"]); + params->m_FlipY = false; if (parsedArgs.count("flip_y")) - flip_y = us::any_cast(parsedArgs["flip_y"]); - bool flip_z = false; + params->m_FlipY = us::any_cast(parsedArgs["flip_y"]); + params->m_FlipZ = false; if (parsedArgs.count("flip_z")) - flip_z = us::any_cast(parsedArgs["flip_z"]); + params->m_FlipZ = us::any_cast(parsedArgs["flip_z"]); bool prior_flip_x = false; if (parsedArgs.count("prior_flip_x")) prior_flip_x = us::any_cast(parsedArgs["prior_flip_x"]); bool prior_flip_y = false; if (parsedArgs.count("prior_flip_y")) prior_flip_y = us::any_cast(parsedArgs["prior_flip_y"]); bool prior_flip_z = false; if (parsedArgs.count("prior_flip_z")) prior_flip_z = us::any_cast(parsedArgs["prior_flip_z"]); - bool apply_image_rotation = false; + params->m_ApplyDirectionMatrix = false; if (parsedArgs.count("apply_image_rotation")) - apply_image_rotation = us::any_cast(parsedArgs["apply_image_rotation"]); + params->m_ApplyDirectionMatrix = us::any_cast(parsedArgs["apply_image_rotation"]); float compress = -1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); - float min_tract_length = 20; + params->m_MinTractLength = 20; if (parsedArgs.count("min_tract_length")) - min_tract_length = us::any_cast(parsedArgs["min_tract_length"]); + params->m_MinTractLength = us::any_cast(parsedArgs["min_tract_length"]); - float loop_check = -1; + params->SetLoopCheckDeg(-1); if (parsedArgs.count("loop_check")) - loop_check = us::any_cast(parsedArgs["loop_check"]); + params->SetLoopCheckDeg(us::any_cast(parsedArgs["loop_check"])); std::string forestFile; if (parsedArgs.count("forest")) forestFile = us::any_cast(parsedArgs["forest"]); std::string maskFile = ""; if (parsedArgs.count("tracking_mask")) maskFile = us::any_cast(parsedArgs["tracking_mask"]); std::string seedFile = ""; if (parsedArgs.count("seed_image")) seedFile = us::any_cast(parsedArgs["seed_image"]); std::string targetFile = ""; if (parsedArgs.count("target_image")) targetFile = us::any_cast(parsedArgs["target_image"]); std::string exclusionFile = ""; if (parsedArgs.count("exclusion_image")) exclusionFile = us::any_cast(parsedArgs["exclusion_image"]); std::string stopFile = ""; if (parsedArgs.count("stop_image")) stopFile = us::any_cast(parsedArgs["stop_image"]); std::string ep_constraint = "NONE"; if (parsedArgs.count("ep_constraint")) ep_constraint = us::any_cast(parsedArgs["ep_constraint"]); - float cutoff = 0.1f; + params->m_Cutoff = 0.1f; if (parsedArgs.count("cutoff")) - cutoff = us::any_cast(parsedArgs["cutoff"]); + params->m_Cutoff = us::any_cast(parsedArgs["cutoff"]); - float odf_cutoff = 0.0; + params->m_OdfCutoff = 0.0; if (parsedArgs.count("odf_cutoff")) - odf_cutoff = us::any_cast(parsedArgs["odf_cutoff"]); + params->m_OdfCutoff = us::any_cast(parsedArgs["odf_cutoff"]); - float stepsize = -1; + params->SetStepSizeVox(-1); if (parsedArgs.count("step_size")) - stepsize = us::any_cast(parsedArgs["step_size"]); + params->SetStepSizeVox(us::any_cast(parsedArgs["step_size"])); - float sampling_distance = -1; + params->SetSamplingDistance(-1); if (parsedArgs.count("sampling_distance")) - sampling_distance = us::any_cast(parsedArgs["sampling_distance"]); + params->SetSamplingDistance(us::any_cast(parsedArgs["sampling_distance"])); - unsigned int num_samples = 0; + params->m_NumSamples = 0; if (parsedArgs.count("num_samples")) - num_samples = static_cast(us::any_cast(parsedArgs["num_samples"])); + params->m_NumSamples = static_cast(us::any_cast(parsedArgs["num_samples"])); - int num_seeds = 1; + params->m_SeedsPerVoxel = 1; if (parsedArgs.count("seeds")) - num_seeds = us::any_cast(parsedArgs["seeds"]); + params->m_SeedsPerVoxel = us::any_cast(parsedArgs["seeds"]); - unsigned int trials_per_seed = 10; + params->m_TrialsPerSeed = 10; if (parsedArgs.count("trials_per_seed")) - trials_per_seed = static_cast(us::any_cast(parsedArgs["trials_per_seed"])); + params->m_TrialsPerSeed = static_cast(us::any_cast(parsedArgs["trials_per_seed"])); - float tend_f = 1; + params->m_F = 1; if (parsedArgs.count("tend_f")) - tend_f = us::any_cast(parsedArgs["tend_f"]); + params->m_F = us::any_cast(parsedArgs["tend_f"]); - float tend_g = 0; + params->m_G = 0; if (parsedArgs.count("tend_g")) - tend_g = us::any_cast(parsedArgs["tend_g"]); + params->m_G = us::any_cast(parsedArgs["tend_g"]); - float angular_threshold = -1; + params->SetAngularThresholdDeg(-1); if (parsedArgs.count("angular_threshold")) - angular_threshold = us::any_cast(parsedArgs["angular_threshold"]); + params->SetAngularThresholdDeg(us::any_cast(parsedArgs["angular_threshold"])); - int max_tracts = -1; + params->m_MaxNumFibers = -1; if (parsedArgs.count("max_tracts")) - max_tracts = us::any_cast(parsedArgs["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 mitkCommandLineParser::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); + // omp_set_num_threads(1); typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); if (!prior_image.empty()) { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, {}); mitk::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_Cutoff = 0.0; + prior_params->m_Mode = mitk::StreamlineTractographyParameters::MODE::DETERMINISTIC; dynamic_cast(priorhandler)->SetPeakImage(itkImg); dynamic_cast(priorhandler)->SetPeakThreshold(0.0); dynamic_cast(priorhandler)->SetInterpolate(interpolate); + dynamic_cast(priorhandler)->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); - priorhandler->SetFlipX(prior_flip_x); - priorhandler->SetFlipY(prior_flip_y); - priorhandler->SetFlipZ(prior_flip_z); - if (algorithm == "ProbPeaks") - priorhandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); - else - priorhandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); - + priorhandler->SetParameters(prior_params); tracker->SetTrackingPriorHandler(priorhandler); tracker->SetTrackingPriorWeight(prior_weight); tracker->SetTrackingPriorAsMask(restrict_to_prior); tracker->SetIntroduceDirectionsFromPrior(new_directions_from_prior); } mitk::TrackingDataHandler* handler; if (algorithm == "DetRF" || algorithm == "ProbRF") { mitk::TractographyForest::Pointer forest = mitk::IOUtil::Load(forestFile); if (forest.IsNull()) mitkThrow() << "Forest file " << forestFile << " could not be read."; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); auto input = mitk::IOUtil::Load(input_files.at(0), &functor); 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); } if (algorithm == "ProbRF") - handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); + params->m_Mode = mitk::TrackingDataHandler::MODE::PROBABILISTIC; } else if (algorithm == "DetPeaks" or algorithm == "ProbPeaks") { handler = new mitk::TrackingHandlerPeaks(); MITK_INFO << "loading input peak image"; mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(input_files.at(0)); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = mitk::convert::GetItkPeakFromPeakImage(mitkImage); if (algorithm == "ProbPeaks") handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); else handler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); dynamic_cast(handler)->SetPeakImage(itkImg); - dynamic_cast(handler)->SetApplyDirectionMatrix(apply_image_rotation); - dynamic_cast(handler)->SetPeakThreshold(cutoff); } else if (algorithm == "DetTensor") { 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)); mitk::TensorImage::ItkTensorImageType::Pointer itkImg = mitk::convert::GetItkTensorFromTensorImage(mitkImage); dynamic_cast(handler)->AddTensorImage(itkImg.GetPointer()); } - dynamic_cast(handler)->SetFaThreshold(cutoff); - dynamic_cast(handler)->SetF(tend_f); - dynamic_cast(handler)->SetG(tend_g); - if (addImages.at(0).size()>0) dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); } else if (algorithm == "DetODF" || algorithm == "ProbODF" || algorithm == "ProbTensor") { handler = new mitk::TrackingHandlerOdf(); mitk::OdfImage::ItkOdfImageType::Pointer itkImg = nullptr; if (algorithm == "ProbTensor") { MITK_INFO << "Converting Tensor to ODF image"; auto input = mitk::IOUtil::Load(input_files.at(0)); itkImg = mitk::convert::GetItkOdfFromTensorImage(input); - sharpen_odfs = true; - odf_cutoff = 0; dynamic_cast(handler)->SetIsOdfFromTensor(true); } else { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image", "ODF Image"}, {}); auto input = mitk::IOUtil::Load(input_files.at(0), &functor)[0]; if (dynamic_cast(input.GetPointer())) { MITK_INFO << "Converting SH to ODF image"; mitk::Image::Pointer mitkImg = dynamic_cast(input.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); - dynamic_cast(handler)->SetGfaThreshold(cutoff); - dynamic_cast(handler)->SetOdfThreshold(odf_cutoff); - dynamic_cast(handler)->SetSharpenOdfs(sharpen_odfs); if (algorithm == "ProbODF" || algorithm == "ProbTensor") - dynamic_cast(handler)->SetMode(mitk::TrackingHandlerOdf::MODE::PROBABILISTIC); + params->m_Mode = mitk::TrackingDataHandler::MODE::PROBABILISTIC; if (addImages.at(0).size()>0) dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); } else { - MITK_INFO << "Unknown tractography algorithm (" + algorithm+"). Known types are DetPeaks, ProbPeaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF."; + MITK_INFO << "Unknown tractography algorithm (" + algorithm+"). Known types are Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF."; return EXIT_FAILURE; } - handler->SetInterpolate(interpolate); - handler->SetFlipX(flip_x); - handler->SetFlipY(flip_y); - handler->SetFlipZ(flip_z); if (ep_constraint=="NONE") - tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NONE); + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::NONE; else if (ep_constraint=="EPS_IN_TARGET") - tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET); + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET; else if (ep_constraint=="EPS_IN_TARGET_LABELDIFF") - tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF); + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF; else if (ep_constraint=="EPS_IN_SEED_AND_TARGET") - tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET); + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET; else if (ep_constraint=="MIN_ONE_EP_IN_TARGET") - tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET); + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET; else if (ep_constraint=="ONE_EP_IN_TARGET") - tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET); + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET; else if (ep_constraint=="NO_EP_IN_TARGET") - tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET); + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET; MITK_INFO << "Tractography algorithm: " << algorithm; - tracker->SetInterpolateMasks(mask_interpolation); - tracker->SetNumberOfSamples(num_samples); - tracker->SetAngularThreshold(angular_threshold); tracker->SetMaskImage(mask); tracker->SetSeedImage(seed); tracker->SetStoppingRegions(stop); tracker->SetTargetRegions(target); tracker->SetExclusionRegions(exclusion); - tracker->SetSeedsPerVoxel(num_seeds); - tracker->SetStepSize(stepsize); - tracker->SetSamplingDistance(sampling_distance); - tracker->SetUseStopVotes(use_stop_votes); - tracker->SetOnlyForwardSamples(use_only_forward_samples); - tracker->SetLoopCheck(loop_check); - tracker->SetMaxNumTracts(max_tracts); - tracker->SetTrialsPerSeed(trials_per_seed); tracker->SetTrackingHandler(handler); if (ext != ".fib" && ext != ".trk") - tracker->SetUseOutputProbabilityMap(true); - tracker->SetMinTractLength(min_tract_length); - tracker->SetRandom(!fix_seed); + params->m_OutputProbMap = true; + tracker->SetParameters(params); tracker->Update(); if (ext == ".fib" || ext == ".trk") { vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); if (compress > 0) outFib->Compress(compress); 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/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp index 94f25ae860..bae0ea1b8d 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp @@ -1,34 +1,27 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTrackingDataHandler.h" namespace mitk { TrackingDataHandler::TrackingDataHandler() - : m_AngularThreshold(0.7) - , m_Interpolate(true) - , m_FlipX(false) - , m_FlipY(false) - , m_FlipZ(false) - , m_Mode(MODE::DETERMINISTIC) - , m_RngItk(ItkRngType::New()) + : m_RngItk(ItkRngType::New()) , m_NeedsDataInit(true) - , m_Random(true) { } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h index 886a8792b3..b2ff010f48 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingDataHandler.h @@ -1,123 +1,105 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _TrackingDataHandler #define _TrackingDataHandler #include #include #include #include #include #include #include #include #include #include +#include namespace mitk { /** * \brief Abstract class for tracking handler. A tracking handler deals with determining the next progression direction of a streamline fiber. There are different handlers for tensor images, peak images, ... */ class MITKFIBERTRACKING_EXPORT TrackingDataHandler { public: - enum MODE { - DETERMINISTIC, - PROBABILISTIC - }; - TrackingDataHandler(); virtual ~TrackingDataHandler(){} typedef itk::Statistics::MersenneTwisterRandomVariateGenerator ItkRngType; typedef boost::mt19937 BoostRngType; typedef itk::Image ItkUcharImgType; typedef itk::Image ItkShortImgType; typedef itk::Image ItkFloatImgType; typedef itk::Image ItkDoubleImgType; typedef vnl_vector_fixed< float, 3 > TrackingDirectionType; + typedef mitk::StreamlineTractographyParameters::MODE MODE; virtual TrackingDirectionType ProposeDirection(const itk::Point& pos, std::deque< TrackingDirectionType >& olddirs, itk::Index<3>& oldIndex) = 0; ///< predicts next progression direction at the given position virtual void InitForTracking() = 0; virtual itk::Vector GetSpacing() = 0; virtual itk::Point GetOrigin() = 0; virtual itk::Matrix GetDirection() = 0; virtual itk::ImageRegion<3> GetLargestPossibleRegion() = 0; virtual bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) = 0; - virtual void SetMode(MODE m) = 0; - MODE GetMode() const { return m_Mode; } - - void SetInterpolate( bool interpolate ){ m_Interpolate = interpolate; } - bool GetInterpolate() const { return m_Interpolate; } - void SetAngularThreshold( float a ){ m_AngularThreshold = a; } - void SetFlipX( bool f ){ m_FlipX = f; } - void SetFlipY( bool f ){ m_FlipY = f; } - void SetFlipZ( bool f ){ m_FlipZ = f; } - void SetRandom( bool random ) + void SetParameters(std::shared_ptr< mitk::StreamlineTractographyParameters > parameters) { - m_Random = random; - if (!random) + m_Parameters = parameters; + + if (m_Parameters->m_FixRandomSeed) { m_Rng.seed(0); std::srand(0); m_RngItk->SetSeed(0); } else { m_Rng.seed(); m_RngItk->SetSeed(); std::srand(std::time(nullptr)); } } double GetRandDouble(const double & a, const double & b) { return m_RngItk->GetUniformVariate(a, b); } - bool GetFlipX() const { return m_FlipX; } - bool GetFlipY() const { return m_FlipY; } - bool GetFlipZ() const { return m_FlipZ; } protected: - float m_AngularThreshold; - bool m_Interpolate; - bool m_FlipX; - bool m_FlipY; - bool m_FlipZ; - MODE m_Mode; BoostRngType m_Rng; ItkRngType::Pointer m_RngItk; bool m_NeedsDataInit; - bool m_Random; + std::shared_ptr< mitk::StreamlineTractographyParameters > m_Parameters; void DataModified() { m_NeedsDataInit = true; } + vnl_matrix_fixed m_FloatImageRotation; + }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp index 2fed13346f..d89a5144b4 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp @@ -1,294 +1,306 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTrackingHandlerOdf.h" #include #include #include #include #include namespace mitk { TrackingHandlerOdf::TrackingHandlerOdf() - : m_GfaThreshold(0.2) - , m_OdfThreshold(0.1) - , m_SharpenOdfs(false) - , m_NumProbSamples(1) + : m_NumProbSamples(1) , m_OdfFromTensor(false) { m_GfaInterpolator = itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::New(); m_OdfInterpolator = itk::LinearInterpolateImageFunction< itk::Image< ItkOdfImageType::PixelType, 3 >, float >::New(); } TrackingHandlerOdf::~TrackingHandlerOdf() { } bool TrackingHandlerOdf::WorldToIndex(itk::Point& pos, itk::Index<3>& index) { m_OdfImage->TransformPhysicalPointToIndex(pos, index); return m_OdfImage->GetLargestPossibleRegion().IsInside(index); } void TrackingHandlerOdf::InitForTracking() { MITK_INFO << "Initializing ODF tracker."; if (m_NeedsDataInit) { m_OdfHemisphereIndices.clear(); itk::OrientationDistributionFunction< float, ODF_SAMPLING_SIZE > odf; vnl_vector_fixed ref; ref.fill(0); ref[0]=1; for (int i=0; i0) m_OdfHemisphereIndices.push_back(i); m_OdfFloatDirs.set_size(m_OdfHemisphereIndices.size(), 3); + auto double_dir = m_OdfImage->GetDirection().GetVnlMatrix(); + for (int r=0; r<3; r++) + for (int c=0; c<3; c++) + { + m_FloatImageRotation[r][c] = double_dir[r][c]; + } + for (unsigned int i=0; i GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(m_OdfImage); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); m_GfaImage = gfaFilter->GetOutput(); } m_NeedsDataInit = false; } + if (m_OdfFromTensor) + { + m_Parameters->m_OdfCutoff = 0; + m_Parameters->m_SharpenOdfs = false; + } + m_GfaInterpolator->SetInputImage(m_GfaImage); m_OdfInterpolator->SetInputImage(m_OdfImage); - std::cout << "TrackingHandlerOdf - GFA threshold: " << m_GfaThreshold << std::endl; - std::cout << "TrackingHandlerOdf - ODF threshold: " << m_OdfThreshold << std::endl; - if (m_SharpenOdfs) + std::cout << "TrackingHandlerOdf - GFA threshold: " << m_Parameters->m_Cutoff << std::endl; + std::cout << "TrackingHandlerOdf - ODF threshold: " << m_Parameters->m_OdfCutoff << std::endl; + if (m_Parameters->m_SharpenOdfs) std::cout << "TrackingHandlerOdf - Sharpening ODfs" << std::endl; } int TrackingHandlerOdf::SampleOdf(vnl_vector< float >& probs, vnl_vector< float >& angles) { boost::random::discrete_distribution dist(probs.begin(), probs.end()); int sampled_idx = 0; int max_sample_idx = -1; float max_prob = 0; int trials = 0; for (int i=0; i> sampler(m_Rng, dist); sampled_idx = sampler(); } - if (probs[sampled_idx]>max_prob && probs[sampled_idx]>m_OdfThreshold && fabs(angles[sampled_idx])>=m_AngularThreshold) + if (probs[sampled_idx]>max_prob && probs[sampled_idx]>m_Parameters->m_OdfCutoff && fabs(angles[sampled_idx])>=m_Parameters->GetAngularThreshold()) { max_prob = probs[sampled_idx]; max_sample_idx = sampled_idx; } - else if ( (probs[sampled_idx]<=m_OdfThreshold || fabs(angles[sampled_idx])m_OdfCutoff || fabs(angles[sampled_idx])GetAngularThreshold()) && trials<50) // we allow 50 trials to exceed the ODF threshold i--; } return max_sample_idx; } void TrackingHandlerOdf::SetIsOdfFromTensor(bool OdfFromTensor) { m_OdfFromTensor = OdfFromTensor; } bool TrackingHandlerOdf::GetIsOdfFromTensor() const { return m_OdfFromTensor; } vnl_vector_fixed TrackingHandlerOdf::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_OdfImage->TransformPhysicalPointToIndex(pos, idx); if ( !m_OdfImage->GetLargestPossibleRegion().IsInside(idx) ) return output_direction; // check GFA threshold for termination - float gfa = mitk::imv::GetImageValue(pos, m_Interpolate, m_GfaInterpolator); - if (gfa(pos, m_Parameters->m_InterpolateTractographyData, m_GfaInterpolator); + if (gfam_Cutoff) return output_direction; vnl_vector_fixed last_dir; if (!olddirs.empty()) last_dir = olddirs.back(); - if (!m_Interpolate && oldIndex==idx) + if (!m_Parameters->m_InterpolateTractographyData && oldIndex==idx) return last_dir; - ItkOdfImageType::PixelType odf_values = mitk::imv::GetImageValue(pos, m_Interpolate, m_OdfInterpolator); + ItkOdfImageType::PixelType odf_values = mitk::imv::GetImageValue(pos, m_Parameters->m_InterpolateTractographyData, m_OdfInterpolator); vnl_vector< float > probs; probs.set_size(m_OdfHemisphereIndices.size()); vnl_vector< float > angles; angles.set_size(m_OdfHemisphereIndices.size()); angles.fill(1.0); // Find ODF maximum and remove <0 values float max_odf_val = 0; float min_odf_val = 999; int max_idx_d = -1; int c = 0; for (int i : m_OdfHemisphereIndices) { if (odf_values[i]<0) odf_values[i] = 0; if (odf_values[i]>max_odf_val) { max_odf_val = odf_values[i]; max_idx_d = c; } if (odf_values[i]m_SharpenOdfs) { // sharpen ODF probs -= min_odf_val; probs /= (max_odf_val-min_odf_val); for (unsigned int i=0; i0) { probs /= odf_sum; max_odf_val /= odf_sum; } } // no previous direction - if (max_odf_val>m_OdfThreshold && (olddirs.empty() || last_dir.magnitude()<=0.5)) + if (max_odf_val>m_Parameters->m_OdfCutoff && (olddirs.empty() || last_dir.magnitude()<=0.5)) { - if (m_Mode==MODE::DETERMINISTIC) // return maximum peak + if (m_Parameters->m_Mode==MODE::DETERMINISTIC) // return maximum peak { output_direction = m_OdfFloatDirs.get_row(max_idx_d); return output_direction * max_odf_val; } - else if (m_Mode==MODE::PROBABILISTIC) // sample from complete ODF + else if (m_Parameters->m_Mode==MODE::PROBABILISTIC) // sample from complete ODF { int max_sample_idx = SampleOdf(probs, angles); if (max_sample_idx>=0) output_direction = m_OdfFloatDirs.get_row(max_sample_idx) * probs[max_sample_idx]; return output_direction; } } - else if (max_odf_val<=m_OdfThreshold) // return (0,0,0) + else if (max_odf_val<=m_Parameters->m_OdfCutoff) // return (0,0,0) { return output_direction; } // correct previous direction - if (m_FlipX) + if (m_Parameters->m_FlipX) last_dir[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) last_dir[1] *= -1; - if (m_FlipZ) + if (m_Parameters->m_FlipZ) last_dir[2] *= -1; // calculate angles between previous direction and ODF directions angles = m_OdfFloatDirs*last_dir; float probs_sum = 0; float max_prob = 0; for (unsigned int i=0; imax_prob && odf_val>m_OdfThreshold) + if (m_Parameters->m_Mode==MODE::DETERMINISTIC && odf_val>max_prob && odf_val>m_Parameters->m_OdfCutoff) { // use maximum peak of the ODF weighted with the directional prior max_prob = odf_val; vnl_vector_fixed d = m_OdfFloatDirs.get_row(i); if (angle<0) d *= -1; output_direction = odf_val*d; } - else if (m_Mode==MODE::PROBABILISTIC) + else if (m_Parameters->m_Mode==MODE::PROBABILISTIC) { // update ODF probabilties with the ODF values pow(abs_angle, m_DirPriorPower) probs[i] = odf_val; probs_sum += probs[i]; } } // do probabilistic sampling - if (m_Mode==MODE::PROBABILISTIC && probs_sum>0.0001) + if (m_Parameters->m_Mode==MODE::PROBABILISTIC && probs_sum>0.0001) { int max_sample_idx = SampleOdf(probs, angles); if (max_sample_idx>=0) { output_direction = m_OdfFloatDirs.get_row(max_sample_idx); if (angles[max_sample_idx]<0) output_direction *= -1; output_direction *= probs[max_sample_idx]; } } // check hard angular threshold float mag = output_direction.magnitude(); if (mag>=0.0001) { output_direction.normalize(); float a = dot_product(output_direction, last_dir); - if (aGetAngularThreshold()) output_direction.fill(0); } else output_direction.fill(0); - if (m_FlipX) + if (m_Parameters->m_FlipX) output_direction[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) output_direction[1] *= -1; - if (m_FlipZ) + if (m_Parameters->m_FlipZ) output_direction[2] *= -1; + if (m_Parameters->m_ApplyDirectionMatrix) + output_direction = m_FloatImageRotation*output_direction; return output_direction; } void TrackingHandlerOdf::SetNumProbSamples(int NumProbSamples) { m_NumProbSamples = NumProbSamples; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h index 593c02520e..4a70fd4f0c 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h @@ -1,87 +1,80 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _TrackingHandlerOdf #define _TrackingHandlerOdf #include "mitkTrackingDataHandler.h" #include #include #include #include namespace mitk { /** * \brief Enables streamline tracking on sampled ODF images. */ class MITKFIBERTRACKING_EXPORT TrackingHandlerOdf : public TrackingDataHandler { public: TrackingHandlerOdf(); ~TrackingHandlerOdf() override; typedef OdfImage::ItkOdfImageType ItkOdfImageType; typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; void InitForTracking() override; ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex) override; ///< predicts next progression direction at the given position bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) override; - void SetSharpenOdfs(bool doSharpen) { m_SharpenOdfs=doSharpen; } - void SetOdfThreshold(float odfThreshold){ m_OdfThreshold = odfThreshold; } - void SetGfaThreshold(float gfaThreshold){ m_GfaThreshold = gfaThreshold; } void SetOdfImage( ItkOdfImageType::Pointer img ){ m_OdfImage = img; DataModified(); } void SetGfaImage( ItkFloatImgType::Pointer img ){ m_GfaImage = img; DataModified(); } - void SetMode( MODE m ) override{ m_Mode = m; } ItkUcharImgType::SpacingType GetSpacing() override{ return m_OdfImage->GetSpacing(); } itk::Point GetOrigin() override{ return m_OdfImage->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection() override{ return m_OdfImage->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion() override{ return m_OdfImage->GetLargestPossibleRegion(); } int OdfPower() const; void SetNumProbSamples(int NumProbSamples); bool GetIsOdfFromTensor() const; void SetIsOdfFromTensor(bool OdfFromTensor); protected: int SampleOdf(vnl_vector< float >& probs, vnl_vector< float >& angles); - float m_GfaThreshold; - float m_OdfThreshold; - bool m_SharpenOdfs; ItkFloatImgType::Pointer m_GfaImage; ///< GFA image used to determine streamline termination. ItkOdfImageType::Pointer m_OdfImage; ///< Input odf image. ItkOdfImageType::Pointer m_WorkingOdfImage; ///< Modified odf image. std::vector< int > m_OdfHemisphereIndices; vnl_matrix< float > m_OdfFloatDirs; int m_NumProbSamples; bool m_OdfFromTensor; itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::Pointer m_GfaInterpolator; itk::LinearInterpolateImageFunction< itk::Image< ItkOdfImageType::PixelType, 3 >, float >::Pointer m_OdfInterpolator; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp index 09a1db5760..623aca610d 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp @@ -1,304 +1,302 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTrackingHandlerPeaks.h" namespace mitk { TrackingHandlerPeaks::TrackingHandlerPeaks() - : m_PeakThreshold(0.1) - , m_ApplyDirectionMatrix(false) { } TrackingHandlerPeaks::~TrackingHandlerPeaks() { } bool TrackingHandlerPeaks::WorldToIndex(itk::Point& pos, itk::Index<3>& index) { m_DummyImage->TransformPhysicalPointToIndex(pos, index); return m_DummyImage->GetLargestPossibleRegion().IsInside(index); } void TrackingHandlerPeaks::InitForTracking() { MITK_INFO << "Initializing peak tracker."; if (m_NeedsDataInit) { itk::Vector spacing4 = m_PeakImage->GetSpacing(); itk::Point origin4 = m_PeakImage->GetOrigin(); itk::Matrix direction4 = m_PeakImage->GetDirection(); itk::ImageRegion<4> imageRegion4 = m_PeakImage->GetLargestPossibleRegion(); spacing3[0] = spacing4[0]; spacing3[1] = spacing4[1]; spacing3[2] = spacing4[2]; origin3[0] = origin4[0]; origin3[1] = origin4[1]; origin3[2] = origin4[2]; for (int r=0; r<3; r++) for (int c=0; c<3; c++) { direction3[r][c] = direction4[r][c]; m_FloatImageRotation[r][c] = direction4[r][c]; } imageRegion3.SetSize(0, imageRegion4.GetSize()[0]); imageRegion3.SetSize(1, imageRegion4.GetSize()[1]); imageRegion3.SetSize(2, imageRegion4.GetSize()[2]); m_DummyImage = ItkUcharImgType::New(); m_DummyImage->SetSpacing( spacing3 ); m_DummyImage->SetOrigin( origin3 ); m_DummyImage->SetDirection( direction3 ); m_DummyImage->SetRegions( imageRegion3 ); m_DummyImage->Allocate(); m_DummyImage->FillBuffer(0.0); m_NumDirs = imageRegion4.GetSize(3)/3; m_NeedsDataInit = false; } - std::cout << "TrackingHandlerPeaks - Peak threshold: " << m_PeakThreshold << std::endl; + std::cout << "TrackingHandlerPeaks - Peak threshold: " << m_Parameters->m_Cutoff << std::endl; } vnl_vector_fixed TrackingHandlerPeaks::GetMatchingDirection(itk::Index<3> idx3, vnl_vector_fixed& oldDir) { vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0; float mag = oldDir.magnitude(); if (magm_FixRandomSeed) { // try m_NumDirs times to get a non-zero random direction for (int j=0; jGetIntegerVariate(m_NumDirs-1); out_dir = GetDirection(idx3, i); if (out_dir.magnitude()>mitk::eps) { oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; found = true; break; } } } if (!found) { // if you didn't find a non-zero random direction, take first non-zero direction you find for (int i=0; imitk::eps) { oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; break; } } } } else { for (int i=0; i dir = GetDirection(idx3, i); mag = dir.magnitude(); if (mag>mitk::eps) dir.normalize(); float a = dot_product(dir, oldDir); if (fabs(a)>angle) { angle = fabs(a); if (a<0) out_dir = -dir; else out_dir = dir; out_dir *= mag; out_dir *= angle; // shrink contribution of direction if is less parallel to previous direction } } } return out_dir; } vnl_vector_fixed TrackingHandlerPeaks::GetDirection(itk::Index<3> idx3, int dirIdx) { vnl_vector_fixed dir; dir.fill(0.0); if ( !m_DummyImage->GetLargestPossibleRegion().IsInside(idx3) ) return dir; PeakImgType::IndexType idx4; idx4.SetElement(0,idx3[0]); idx4.SetElement(1,idx3[1]); idx4.SetElement(2,idx3[2]); for (int k=0; k<3; k++) { idx4.SetElement(3, dirIdx*3 + k); dir[k] = m_PeakImage->GetPixel(idx4); } - if (m_FlipX) + if (m_Parameters->m_FlipX) dir[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) dir[1] *= -1; - if (m_FlipZ) + if (m_Parameters->m_FlipZ) dir[2] *= -1; - if (m_ApplyDirectionMatrix) + if (m_Parameters->m_ApplyDirectionMatrix) dir = m_FloatImageRotation*dir; return dir; } vnl_vector_fixed TrackingHandlerPeaks::GetDirection(itk::Point itkP, bool interpolate, vnl_vector_fixed oldDir){ // transform physical point to index coordinates itk::Index<3> idx3; itk::ContinuousIndex< float, 3> cIdx; m_DummyImage->TransformPhysicalPointToIndex(itkP, idx3); m_DummyImage->TransformPhysicalPointToContinuousIndex(itkP, cIdx); vnl_vector_fixed dir; dir.fill(0.0); if ( !m_DummyImage->GetLargestPossibleRegion().IsInside(idx3) ) return dir; if (interpolate) { float frac_x = cIdx[0] - idx3[0]; float frac_y = cIdx[1] - idx3[1]; float frac_z = cIdx[2] - idx3[2]; if (frac_x<0) { idx3[0] -= 1; frac_x += 1; } if (frac_y<0) { idx3[1] -= 1; frac_y += 1; } if (frac_z<0) { idx3[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx3[0] >= 0 && idx3[0] < static_cast(m_DummyImage->GetLargestPossibleRegion().GetSize(0) - 1) && idx3[1] >= 0 && idx3[1] < static_cast(m_DummyImage->GetLargestPossibleRegion().GetSize(1) - 1) && idx3[2] >= 0 && idx3[2] < static_cast(m_DummyImage->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); dir = GetMatchingDirection(idx3, oldDir) * interpWeights[0]; itk::Index<3> tmpIdx = idx3; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[1]; tmpIdx = idx3; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[2]; tmpIdx = idx3; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[3]; tmpIdx = idx3; tmpIdx[0]++; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[4]; tmpIdx = idx3; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[5]; tmpIdx = idx3; tmpIdx[2]++; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[6]; tmpIdx = idx3; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[7]; } } else dir = GetMatchingDirection(idx3, oldDir); return dir; } vnl_vector_fixed TrackingHandlerPeaks::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { // CHECK: wann wird wo normalisiert vnl_vector_fixed output_direction; output_direction.fill(0); itk::Index<3> index; m_DummyImage->TransformPhysicalPointToIndex(pos, index); vnl_vector_fixed oldDir; oldDir.fill(0.0); if (!olddirs.empty()) oldDir = olddirs.back(); float old_mag = oldDir.magnitude(); - if (!m_Interpolate && oldIndex==index) + if (!m_Parameters->m_InterpolateTractographyData && oldIndex==index) return oldDir; - output_direction = GetDirection(pos, m_Interpolate, oldDir); + output_direction = GetDirection(pos, m_Parameters->m_InterpolateTractographyData, oldDir); float mag = output_direction.magnitude(); - if (mag>=m_PeakThreshold) + if (mag>=m_Parameters->m_Cutoff) { if (m_Mode == MODE::PROBABILISTIC) { output_direction[0] += this->m_RngItk->GetNormalVariate(0, fabs(output_direction[0])*0.01); output_direction[1] += this->m_RngItk->GetNormalVariate(0, fabs(output_direction[1])*0.01); output_direction[2] += this->m_RngItk->GetNormalVariate(0, fabs(output_direction[2])*0.01); mag = output_direction.magnitude(); } output_direction.normalize(); float a = 1; if (old_mag>0.5) a = dot_product(output_direction, oldDir); - if (a>=m_AngularThreshold) + if (a>=m_Parameters->GetAngularThreshold()) output_direction *= mag; else output_direction.fill(0); } else output_direction.fill(0); return output_direction; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h index c4e92f3b85..6143ac66f2 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h @@ -1,81 +1,71 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _TrackingHandlerPeaks #define _TrackingHandlerPeaks #include "mitkTrackingDataHandler.h" #include #include namespace mitk { /** * \brief Enables deterministic streamline tracking on MRtrix style peak images (4D float images) */ class MITKFIBERTRACKING_EXPORT TrackingHandlerPeaks : public TrackingDataHandler { public: TrackingHandlerPeaks(); ~TrackingHandlerPeaks() override; typedef itk::Image< float, 4 > PeakImgType; ///< MRtrix peak image type void InitForTracking() override; ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex) override; ///< predicts next progression direction at the given position bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) override; - void SetPeakThreshold(float thr){ m_PeakThreshold = thr; } void SetPeakImage( PeakImgType::Pointer image ){ m_PeakImage = image; DataModified(); } - void SetApplyDirectionMatrix( bool applyDirectionMatrix ){ m_ApplyDirectionMatrix = applyDirectionMatrix; } itk::Vector GetSpacing() override{ return spacing3; } itk::Point GetOrigin() override{ return origin3; } itk::Matrix GetDirection() override{ return direction3; } itk::ImageRegion<3> GetLargestPossibleRegion() override{ return imageRegion3; } - void SetMode( MODE m ) override - { - m_Mode = m; - } protected: vnl_vector_fixed GetDirection(itk::Point itkP, bool interpolate, vnl_vector_fixed oldDir); vnl_vector_fixed GetMatchingDirection(itk::Index<3> idx3, vnl_vector_fixed& oldDir); vnl_vector_fixed GetDirection(itk::Index<3> idx3, int dirIdx); PeakImgType::ConstPointer m_PeakImage; - float m_PeakThreshold; int m_NumDirs; itk::Vector spacing3; itk::Point origin3; itk::Matrix direction3; itk::ImageRegion<3> imageRegion3; - vnl_matrix_fixed m_FloatImageRotation; ItkUcharImgType::Pointer m_DummyImage; - - bool m_ApplyDirectionMatrix; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp index dfeee974d3..b6be79a25c 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.cpp @@ -1,900 +1,908 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #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_NumPreviousDirections(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; } } 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_Interpolate && oldIndex==idx) + 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_Interpolate, m_DwiFeatureImageInterpolator); + 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_FlipX) + if (m_Parameters->m_FlipX) tempD[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) tempD[1] *= -1; - if (m_FlipZ) + 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_NumPreviousDirections*3+c) = v; + 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_Mode==MODE::DETERMINISTIC) + 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_AngularThreshold) // is angle between the directions smaller than our hard threshold? + if (abs_angle>=m_Parameters->GetAngularThreshold()) // 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_Mode==MODE::PROBABILISTIC && pNonFib<0.5) + 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_AngularThreshold)) ) + if ( probs2[sampled_idx]>0.1 && (!check_last_dir || (check_last_dir && fabs(angles[sampled_idx])>=m_Parameters->GetAngularThreshold())) ) 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_FlipX) + if (m_Parameters->m_FlipX) output_direction[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) output_direction[1] *= -1; - if (m_FlipZ) + 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_NumPreviousDirections+additional_features) ) - MITK_ERROR << "Wrong number of features in forest: got " << m_Forest->GetNumFeatures() << ", expected " << (NumberOfSignalFeatures+3*m_NumPreviousDirections+additional_features); + 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_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->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_NumPreviousDirections+1); + 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_NumPreviousDirections*3+m_AdditionalFeatureImages.at(0).size()) ); + 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_NumPreviousDirections); // how many directions should be zero? - 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_NumPreviousDirections*3+add_feat_c) = v; + 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_NumPreviousDirections); ++n) + for (int n = 0; n <= static_cast(m_Parameters->m_NumPreviousDirections); ++n) { if (!m_ZeroDirWmFeatures) - n = m_NumPreviousDirections; + n = m_Parameters->m_NumPreviousDirections; for (bool reverse : {false, true}) { for (int j=1; j itkP1, itkP2; - int num_nonzero_dirs = m_NumPreviousDirections; + 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_NumPreviousDirections-num_nonzero_dirs); fm_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_Interpolate, dwi_interp); + 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/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h index 5b911770d4..136c6be6f0 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h @@ -1,166 +1,163 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _TrackingForestHandler #define _TrackingForestHandler #include #include #include #include #include #include #include #include #undef DIFFERENCE #define VIGRA_STATIC_LIB #include #include #include #include #include #include #define _USE_MATH_DEFINES #include namespace mitk { /** * \brief Manages random forests for fiber tractography. The preparation of the features from the inputa data and the training process are handled here. The data preprocessing and actual prediction for the tracking process is also performed here. The tracking itself is performed in MLBSTrackingFilter. */ template< int ShOrder, int NumberOfSignalFeatures > class TrackingHandlerRandomForest : public TrackingDataHandler { public: TrackingHandlerRandomForest(); ~TrackingHandlerRandomForest() override; typedef itk::Image< itk::Vector< float, NumberOfSignalFeatures > , 3 > DwiFeatureImageType; typedef itk::LinearInterpolateImageFunction< DwiFeatureImageType, float > DwiFeatureImageInterpolatorType; typedef itk::LinearInterpolateImageFunction< ItkFloatImgType, float > FloatImageInterpolatorType; typedef mitk::ThresholdSplit >,int,vigra::ClassificationTag> DefaultSplitType; void SetDwis( std::vector< Image::Pointer > images ){ m_InputDwis = images; DataModified(); } void AddDwi( Image::Pointer img ){ m_InputDwis.push_back(img); DataModified(); } void SetTractograms( std::vector< FiberBundle::Pointer > tractograms ) { m_Tractograms.clear(); for (auto fib : tractograms) m_Tractograms.push_back(fib); DataModified(); } void SetMaskImages( std::vector< ItkUcharImgType::Pointer > images ){ m_MaskImages = images; DataModified(); } void SetWhiteMatterImages( std::vector< ItkUcharImgType::Pointer > images ){ m_WhiteMatterImages = images; DataModified(); } void SetFiberVolumeModImages( std::vector< ItkFloatImgType::Pointer > images ){ m_FiberVolumeModImages = images; DataModified(); } void SetAdditionalFeatureImages( std::vector< std::vector< ItkFloatImgType::Pointer > > images ){ m_AdditionalFeatureImages = images; DataModified(); } void StartTraining(); - void SetMode( MODE m ) override{ m_Mode = m; } void SetForest(mitk::TractographyForest::Pointer forest){ m_Forest = forest; } void SetMaxNumWmSamples(int num){ m_MaxNumWmSamples=num; } - void SetNumPreviousDirections( unsigned int num ){ m_NumPreviousDirections=num; } void SetNumTrees(int num){ m_NumTrees = num; } void SetMaxTreeDepth(int depth){ m_MaxTreeDepth = depth; } void SetFiberSamplingStep(float step){ m_WmSampleDistance = step; } void SetGrayMatterSamplesPerVoxel(int samples){ m_GmSamplesPerVoxel = samples; } void SetSampleFraction(float fraction){ m_SampleFraction = fraction; } void SetBidirectionalFiberSampling(bool val) { m_BidirectionalFiberSampling = val; } void SetZeroDirWmFeatures(bool val) { m_ZeroDirWmFeatures = val; } void InitForTracking() override; ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex) override; ///< predicts next progression direction at the given position bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) override; bool IsForestValid(); ///< true is forest is not null, has more than 0 trees and the correct number of features (NumberOfSignalFeatures + 3) mitk::TractographyForest::Pointer GetForest(){ return m_Forest; } ItkUcharImgType::SpacingType GetSpacing() override{ return m_DwiFeatureImages.at(0)->GetSpacing(); } itk::Point GetOrigin() override{ return m_DwiFeatureImages.at(0)->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection() override{ return m_DwiFeatureImages.at(0)->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion() override{ return m_DwiFeatureImages.at(0)->GetLargestPossibleRegion(); } protected: void InputDataValidForTracking(); ///< check if raw data is set and tracking forest is valid template typename std::enable_if::type InitDwiImageFeatures(mitk::Image::Pointer mitk_dwi); template typename std::enable_if= 100, T>::type InitDwiImageFeatures(mitk::Image::Pointer mitk_dwi); void InputDataValidForTraining(); ///< Check if everything is tehere for training (raw datasets, fiber tracts) void InitForTraining(); ///< Generate masks if necessary, resample fibers, spherically interpolate raw DWIs void CalculateTrainingSamples(); ///< Calculate GM and WM features using the interpolated raw data, the WM masks and the fibers std::vector< Image::Pointer > m_InputDwis; ///< original input DWI data mitk::TractographyForest::Pointer m_Forest; ///< random forest classifier std::chrono::time_point m_StartTime; std::chrono::time_point m_EndTime; std::vector< typename DwiFeatureImageType::Pointer > m_DwiFeatureImages; std::vector< std::vector< ItkFloatImgType::Pointer > > m_AdditionalFeatureImages; std::vector< ItkFloatImgType::Pointer > m_FiberVolumeModImages; ///< used to correct the fiber density std::vector< FiberBundle::Pointer > m_Tractograms; ///< training tractograms std::vector< ItkUcharImgType::Pointer > m_MaskImages; ///< binary mask images to constrain training to a certain area (e.g. brain mask) std::vector< ItkUcharImgType::Pointer > m_WhiteMatterImages; ///< defines white matter voxels. if not set, theses mask images are automatically generated from the input tractograms float m_WmSampleDistance; ///< deterines the number of white matter samples (distance of sampling points on each fiber). int m_NumTrees; ///< number of trees in random forest int m_MaxTreeDepth; ///< limits the tree depth float m_SampleFraction; ///< fraction of samples used to train each tree unsigned int m_NumberOfSamples; ///< stores overall number of samples used for training std::vector< unsigned int > m_GmSamples; ///< number of gray matter samples int m_GmSamplesPerVoxel; ///< number of gray matter samplees per voxel. if -1, then the number is automatically chosen to gain an overall number of GM samples close to the number of WM samples. vigra::MultiArray<2, float> m_FeatureData; ///< vigra container for training features - unsigned int m_NumPreviousDirections; ///< How many "old" directions should be used as classification features? // only for tracking vigra::MultiArray<2, float> m_LabelData; ///< vigra container for training labels vigra::MultiArray<2, double> m_Weights; ///< vigra container for training sample weights std::vector< vnl_vector_fixed > m_DirectionContainer; vnl_matrix< float > m_OdfFloatDirs; bool m_BidirectionalFiberSampling; bool m_ZeroDirWmFeatures; int m_MaxNumWmSamples; std::vector< std::vector< bool > > m_SampleUsage; vnl_matrix_fixed m_ImageDirection; vnl_matrix_fixed m_ImageDirectionInverse; typename DwiFeatureImageInterpolatorType::Pointer m_DwiFeatureImageInterpolator; std::vector< std::vector< FloatImageInterpolatorType::Pointer > > m_AdditionalFeatureImageInterpolators; }; } #include "mitkTrackingHandlerRandomForest.cpp" #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp index 0581e41ec7..e21d5651ac 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp @@ -1,371 +1,377 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTrackingHandlerTensor.h" namespace mitk { TrackingHandlerTensor::TrackingHandlerTensor() - : m_FaThreshold(0.1) - , m_F(1.0) - , m_G(0.0) - , m_InterpolateTensors(true) + : m_InterpolateTensors(true) , m_NumberOfInputs(0) { m_FaInterpolator = itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::New(); } TrackingHandlerTensor::~TrackingHandlerTensor() { } void TrackingHandlerTensor::InitForTracking() { MITK_INFO << "Initializing tensor tracker."; if (m_NeedsDataInit) { m_NumberOfInputs = m_TensorImages.size(); for (int i=0; iSetSpacing( m_TensorImages.at(0)->GetSpacing() ); pdImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); pdImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); pdImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); pdImage->Allocate(); m_PdImage.push_back(pdImage); ItkDoubleImgType::Pointer emaxImage = ItkDoubleImgType::New(); emaxImage->SetSpacing( m_TensorImages.at(0)->GetSpacing() ); emaxImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); emaxImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); emaxImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); emaxImage->Allocate(); emaxImage->FillBuffer(1.0); m_EmaxImage.push_back(emaxImage); } bool useUserFaImage = true; if (m_FaImage.IsNull()) { m_FaImage = ItkFloatImgType::New(); m_FaImage->SetSpacing( m_TensorImages.at(0)->GetSpacing() ); m_FaImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); m_FaImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); m_FaImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); m_FaImage->Allocate(); m_FaImage->FillBuffer(0.0); useUserFaImage = false; } typedef itk::DiffusionTensor3D TensorType; for (int x=0; x<(int)m_TensorImages.at(0)->GetLargestPossibleRegion().GetSize()[0]; x++) for (int y=0; y<(int)m_TensorImages.at(0)->GetLargestPossibleRegion().GetSize()[1]; y++) for (int z=0; z<(int)m_TensorImages.at(0)->GetLargestPossibleRegion().GetSize()[2]; z++) { ItkTensorImageType::IndexType index; index[0] = x; index[1] = y; index[2] = z; for (int i=0; i dir; tensor = m_TensorImages.at(i)->GetPixel(index); tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); if (dir.magnitude()>mitk::eps) dir.normalize(); else dir.fill(0.0); m_PdImage.at(i)->SetPixel(index, dir); if (!useUserFaImage) m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)+tensor.GetFractionalAnisotropy()); m_EmaxImage.at(i)->SetPixel(index, 2/eigenvalues[2]); } if (!useUserFaImage) m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)/m_NumberOfInputs); } + auto double_dir = m_TensorImages.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; } - if (m_F+m_G>1.0) + if (m_Parameters->m_F+m_Parameters->m_G>1.0) { - float temp = m_F+m_G; - m_F /= temp; - m_G /= temp; + float temp = m_Parameters->m_F+m_Parameters->m_G; + m_Parameters->m_F /= temp; + m_Parameters->m_G /= temp; } m_FaInterpolator->SetInputImage(m_FaImage); - std::cout << "TrackingHandlerTensor - FA threshold: " << m_FaThreshold << std::endl; - std::cout << "TrackingHandlerTensor - f: " << m_F << std::endl; - std::cout << "TrackingHandlerTensor - g: " << m_G << std::endl; + std::cout << "TrackingHandlerTensor - FA threshold: " << m_Parameters->m_Cutoff << std::endl; + std::cout << "TrackingHandlerTensor - f: " << m_Parameters->m_F << std::endl; + std::cout << "TrackingHandlerTensor - g: " << m_Parameters->m_G << std::endl; } vnl_vector_fixed TrackingHandlerTensor::GetMatchingDirection(itk::Index<3> idx, vnl_vector_fixed& oldDir, int& image_num) { vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0; float mag = oldDir.magnitude(); if (magGetPixel(idx); if (out_dir.magnitude()>0.5) { image_num = i; oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; break; } } } else { for (unsigned int i=0; i dir = m_PdImage.at(i)->GetPixel(idx); float a = dot_product(dir, oldDir); if (fabs(a)>angle) { image_num = i; angle = fabs(a); if (a<0) out_dir = -dir; else out_dir = dir; out_dir *= angle; // shrink contribution of direction if is less parallel to previous direction } } } return out_dir; } bool TrackingHandlerTensor::WorldToIndex(itk::Point& pos, itk::Index<3>& index) { m_TensorImages.at(0)->TransformPhysicalPointToIndex(pos, index); return m_TensorImages.at(0)->GetLargestPossibleRegion().IsInside(index); } vnl_vector_fixed TrackingHandlerTensor::GetDirection(itk::Point itkP, vnl_vector_fixed oldDir, TensorType& tensor) { // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; m_FaImage->TransformPhysicalPointToIndex(itkP, idx); m_FaImage->TransformPhysicalPointToContinuousIndex(itkP, cIdx); vnl_vector_fixed dir; dir.fill(0.0); if ( !m_FaImage->GetLargestPossibleRegion().IsInside(idx) ) return dir; int image_num = -1; - if (!m_Interpolate) + if (!m_Parameters->m_InterpolateTractographyData) { dir = GetMatchingDirection(idx, oldDir, image_num); if (image_num>=0) tensor = m_TensorImages[image_num]->GetPixel(idx) * m_EmaxImage[image_num]->GetPixel(idx); } else { float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx[0] >= 0 && idx[0] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(0) - 1) && idx[1] >= 0 && idx[1] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(1) - 1) && idx[2] >= 0 && idx[2] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); dir = GetMatchingDirection(idx, oldDir, image_num) * interpWeights[0]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(idx) * m_EmaxImage[image_num]->GetPixel(idx) * interpWeights[0]; itk::Index<3> tmpIdx = idx; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[1]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[2]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[3]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[4]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[5]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[6]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[7]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[7]; } } return dir; } vnl_vector_fixed TrackingHandlerTensor::GetLargestEigenvector(TensorType& tensor) { vnl_vector_fixed dir; TensorType::EigenValuesArrayType eigenvalues; TensorType::EigenVectorsMatrixType eigenvectors; tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); if (dir.magnitude() TrackingHandlerTensor::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { vnl_vector_fixed output_direction; output_direction.fill(0); TensorType tensor; tensor.Fill(0); try { itk::Index<3> index; m_TensorImages.at(0)->TransformPhysicalPointToIndex(pos, index); - float fa = mitk::imv::GetImageValue(pos, m_Interpolate, m_FaInterpolator); - if (fa(pos, m_Parameters->m_InterpolateTractographyData, m_FaInterpolator); + if (fam_Cutoff) return output_direction; vnl_vector_fixed oldDir; oldDir.fill(0.0); if (!olddirs.empty()) oldDir = olddirs.back(); - if (m_FlipX) + if (m_Parameters->m_FlipX) oldDir[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) oldDir[1] *= -1; - if (m_FlipZ) + if (m_Parameters->m_FlipZ) oldDir[2] *= -1; float old_mag = oldDir.magnitude(); - if (!m_Interpolate && oldIndex==index) + if (!m_Parameters->m_InterpolateTractographyData && oldIndex==index) return oldDir; output_direction = GetDirection(pos, oldDir, tensor); float mag = output_direction.magnitude(); if (mag>=mitk::eps) { output_direction.normalize(); - if (old_mag>0.5 && m_G>mitk::eps) // TEND tracking + if (old_mag>0.5 && m_Parameters->m_G>mitk::eps) // TEND tracking { - output_direction[0] = m_F*output_direction[0] + (1-m_F)*( (1-m_G)*oldDir[0] + m_G*(tensor[0]*oldDir[0] + tensor[1]*oldDir[1] + tensor[2]*oldDir[2])); - output_direction[1] = m_F*output_direction[1] + (1-m_F)*( (1-m_G)*oldDir[1] + m_G*(tensor[1]*oldDir[0] + tensor[3]*oldDir[1] + tensor[4]*oldDir[2])); - output_direction[2] = m_F*output_direction[2] + (1-m_F)*( (1-m_G)*oldDir[2] + m_G*(tensor[2]*oldDir[0] + tensor[4]*oldDir[1] + tensor[5]*oldDir[2])); + output_direction[0] = m_Parameters->m_F*output_direction[0] + (1-m_Parameters->m_F)*( (1-m_Parameters->m_G)*oldDir[0] + m_Parameters->m_G*(tensor[0]*oldDir[0] + tensor[1]*oldDir[1] + tensor[2]*oldDir[2])); + output_direction[1] = m_Parameters->m_F*output_direction[1] + (1-m_Parameters->m_F)*( (1-m_Parameters->m_G)*oldDir[1] + m_Parameters->m_G*(tensor[1]*oldDir[0] + tensor[3]*oldDir[1] + tensor[4]*oldDir[2])); + output_direction[2] = m_Parameters->m_F*output_direction[2] + (1-m_Parameters->m_F)*( (1-m_Parameters->m_G)*oldDir[2] + m_Parameters->m_G*(tensor[2]*oldDir[0] + tensor[4]*oldDir[1] + tensor[5]*oldDir[2])); output_direction.normalize(); } float a = 1; if (old_mag>0.5) a = dot_product(output_direction, oldDir); - if (a>=m_AngularThreshold) + if (a>=m_Parameters->GetAngularThreshold()) output_direction *= mag; else output_direction.fill(0); } else output_direction.fill(0); } catch(...) { } - if (m_FlipX) + if (m_Parameters->m_FlipX) output_direction[0] *= -1; - if (m_FlipY) + if (m_Parameters->m_FlipY) output_direction[1] *= -1; - if (m_FlipZ) + if (m_Parameters->m_FlipZ) output_direction[2] *= -1; + if (m_Parameters->m_ApplyDirectionMatrix) + output_direction = m_FloatImageRotation*output_direction; return output_direction; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h index 5884ee8be1..6175bec856 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h @@ -1,94 +1,80 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _TrackingHandlerTensor #define _TrackingHandlerTensor #include "mitkTrackingDataHandler.h" #include #include #include #include namespace mitk { /** * \brief Enables streamline tracking on tensor images. Supports multi tensor tracking by adding multiple tensor images. */ class MITKFIBERTRACKING_EXPORT TrackingHandlerTensor : public TrackingDataHandler { public: TrackingHandlerTensor(); ~TrackingHandlerTensor() override; typedef TensorImage::PixelType TensorType; typedef TensorImage::ItkTensorImageType ItkTensorImageType; typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; void InitForTracking() override; ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex) override; ///< predicts next progression direction at the given position bool WorldToIndex(itk::Point& pos, itk::Index<3>& index) override; - void SetF(float f){ m_F = f; } - void SetG(float g){ m_G = g; } - void SetFaThreshold(float FaThreshold){ m_FaThreshold = FaThreshold; } void AddTensorImage( ItkTensorImageType::ConstPointer img ){ m_TensorImages.push_back(img); DataModified(); } void SetTensorImage( ItkTensorImageType::ConstPointer img ){ m_TensorImages.clear(); m_TensorImages.push_back(img); DataModified(); } void ClearTensorImages(){ m_TensorImages.clear(); DataModified(); } void SetFaImage( ItkFloatImgType::Pointer img ){ m_FaImage = img; DataModified(); } void SetInterpolateTensors( bool interpolateTensors ){ m_InterpolateTensors = interpolateTensors; } - void SetMode( MODE m ) override - { - if (m==MODE::DETERMINISTIC) - m_Mode = m; - else - mitkThrow() << "Tensor tracker is only implemented for deterministic mode."; - } int GetNumTensorImages() const { return m_TensorImages.size(); } ItkUcharImgType::SpacingType GetSpacing() override{ return m_FaImage->GetSpacing(); } itk::Point GetOrigin() override{ return m_FaImage->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection() override{ return m_FaImage->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion() override{ return m_FaImage->GetLargestPossibleRegion(); } protected: vnl_vector_fixed GetMatchingDirection(itk::Index<3> idx, vnl_vector_fixed& oldDir, int& image_num); vnl_vector_fixed GetDirection(itk::Point itkP, vnl_vector_fixed oldDir, TensorType& tensor); vnl_vector_fixed GetLargestEigenvector(TensorType& tensor); - float m_FaThreshold; - float m_F; - float m_G; - std::vector< ItkDoubleImgType::Pointer > m_EmaxImage; ///< Stores largest eigenvalues per voxel (one for each tensor) ItkFloatImgType::Pointer m_FaImage; ///< FA image used to determine streamline termination. std::vector< ItkPDImgType::Pointer > m_PdImage; ///< Stores principal direction of each tensor in each voxel. std::vector< ItkTensorImageType::ConstPointer > m_TensorImages; ///< Input tensor images. For multi tensor tracking provide multiple tensor images. bool m_InterpolateTensors; ///< If false, then the peaks are interpolated. Otherwiese, The tensors are interpolated. int m_NumberOfInputs; itk::LinearInterpolateImageFunction< itk::Image< float, 3 >, float >::Pointer m_FaInterpolator; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 7d76d2624c..92fe1c7996 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,1025 +1,986 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #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_MinVoxelSize(-1) - , m_AngularThresholdDeg(-1) - , m_StepSizeVox(-1) - , m_SamplingDistanceVox(-1) - , m_AngularThreshold(-1) - , m_StepSize(0) - , m_MaxLength(10000) - , m_MinTractLength(20.0) - , m_MaxTractLength(400.0) - , m_SeedsPerVoxel(1) - , m_AvoidStop(true) - , m_RandomSampling(false) - , m_SamplingDistance(-1) - , m_DeflectionMod(1.0) - , m_OnlyForwardSamples(true) - , m_UseStopVotes(true) - , m_NumberOfSamples(30) - , m_NumPreviousDirections(1) - , m_MaxNumTracts(-1) , m_Verbose(true) - , m_LoopCheck(-1) , m_DemoMode(false) - , m_Random(true) - , m_UseOutputProbabilityMap(false) , m_CurrentTracts(0) , m_Progress(0) , m_StopTracking(false) - , m_InterpolateMasks(true) - , m_TrialsPerSeed(10) - , m_EndpointConstraint(EndpointConstraints::NONE) - , m_IntroduceDirectionsFromPrior(true) - , m_TrackingPriorAsMask(true) - , m_TrackingPriorWeight(1.0) , 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_MaxNumTracts>0) - status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts) + "/" + boost::lexical_cast(m_MaxNumTracts); + 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->SetRandom(m_Random); - m_TrackingHandler->InitForTracking(); - m_FiberPolyData = PolyDataType::New(); - m_Points = vtkSmartPointer< vtkPoints >::New(); - m_Cells = vtkSmartPointer< vtkCellArray >::New(); - itk::Vector< double, 3 > imageSpacing = m_TrackingHandler->GetSpacing(); - if(imageSpacing[0](imageSpacing[0]); else if (imageSpacing[1] < imageSpacing[2]) m_MinVoxelSize = static_cast(imageSpacing[1]); else m_MinVoxelSize = static_cast(imageSpacing[2]); - if (m_StepSizeVox(mitk::eps)) - m_StepSize = 0.5f*m_MinVoxelSize; - else - m_StepSize = m_StepSizeVox*m_MinVoxelSize; + m_Parameters->SetMinVoxelSize(m_MinVoxelSize); - if (m_AngularThresholdDeg<0) - { - if (m_StepSize/m_MinVoxelSize<=0.966f) // minimum 15° for automatic estimation - m_AngularThreshold = static_cast(std::cos( 0.5 * itk::Math::pi * static_cast(m_StepSize/m_MinVoxelSize) )); - else - m_AngularThreshold = static_cast(std::cos( 0.5 * itk::Math::pi * 0.966 )); - } - else - m_AngularThreshold = static_cast(std::cos( static_cast(m_AngularThresholdDeg)*itk::Math::pi/180.0 )); - m_TrackingHandler->SetAngularThreshold(m_AngularThreshold); + 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->SetRandom(m_Random); m_TrackingPriorHandler->InitForTracking(); - m_TrackingPriorHandler->SetAngularThreshold(m_AngularThreshold); } - if (m_SamplingDistanceVox(mitk::eps)) - m_SamplingDistance = m_MinVoxelSize*0.25f; - else - m_SamplingDistance = m_SamplingDistanceVox * m_MinVoxelSize; - m_PolyDataContainer.clear(); for (unsigned int i=0; iGetNumberOfThreads(); i++) { PolyDataType poly = PolyDataType::New(); m_PolyDataContainer.push_back(poly); } - if (m_UseOutputProbabilityMap) + 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_EndpointConstraint==EndpointConstraints::NONE && m_TargetImageSet && m_SeedImageSet) + 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_EndpointConstraint = EndpointConstraints::EPS_IN_SEED_AND_TARGET; + m_Parameters->m_EpConstraints = EndpointConstraints::EPS_IN_SEED_AND_TARGET; } - else if (m_EndpointConstraint==EndpointConstraints::NONE && m_TargetImageSet) + 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_EndpointConstraint = EndpointConstraints::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_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::DETERMINISTIC) + if (m_Parameters->m_Mode==mitk::TrackingDataHandler::MODE::DETERMINISTIC) std::cout << "StreamlineTracking - Mode: deterministic" << std::endl; - else if(m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::PROBABILISTIC) + else if(m_Parameters->m_Mode==mitk::TrackingDataHandler::MODE::PROBABILISTIC) { std::cout << "StreamlineTracking - Mode: probabilistic" << std::endl; - std::cout << "StreamlineTracking - Trials per seed: " << m_TrialsPerSeed << std::endl; + std::cout << "StreamlineTracking - Trials per seed: " << m_Parameters->m_TrialsPerSeed << std::endl; } else std::cout << "StreamlineTracking - Mode: ???" << std::endl; - if (m_EndpointConstraint==EndpointConstraints::NONE) + if (m_Parameters->m_EpConstraints==EndpointConstraints::NONE) std::cout << "StreamlineTracking - Endpoint constraint: NONE" << std::endl; - else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET) + else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_TARGET) std::cout << "StreamlineTracking - Endpoint constraint: EPS_IN_TARGET" << std::endl; - else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET_LABELDIFF) + 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_EndpointConstraint==EndpointConstraints::EPS_IN_SEED_AND_TARGET) + 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_EndpointConstraint==EndpointConstraints::MIN_ONE_EP_IN_TARGET) + 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_EndpointConstraint==EndpointConstraints::ONE_EP_IN_TARGET) + 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_EndpointConstraint==EndpointConstraints::NO_EP_IN_TARGET) + 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_AngularThreshold << " (" << 180*std::acos( static_cast(m_AngularThreshold) )/itk::Math::pi << "°)" << std::endl; - std::cout << "StreamlineTracking - Stepsize: " << m_StepSize << "mm (" << m_StepSize/m_MinVoxelSize << "*vox)" << std::endl; - std::cout << "StreamlineTracking - Seeds per voxel: " << m_SeedsPerVoxel << std::endl; - std::cout << "StreamlineTracking - Max. tract length: " << m_MaxTractLength << "mm" << std::endl; - std::cout << "StreamlineTracking - Min. tract length: " << m_MinTractLength << "mm" << std::endl; - std::cout << "StreamlineTracking - Max. num. tracts: " << m_MaxNumTracts << std::endl; - std::cout << "StreamlineTracking - Loop check: " << m_LoopCheck << "°" << std::endl; + std::cout << "StreamlineTracking - Angular threshold: " << m_Parameters->GetAngularThreshold() << "°" << std::endl; + std::cout << "StreamlineTracking - Stepsize: " << m_Parameters->GetStepSize() << "mm (" << m_Parameters->GetStepSize()/m_MinVoxelSize << "*vox)" << std::endl; + std::cout << "StreamlineTracking - Seeds per voxel: " << m_Parameters->m_SeedsPerVoxel << std::endl; + std::cout << "StreamlineTracking - Max. tract length: " << m_Parameters->m_MaxTractLength << "mm" << std::endl; + std::cout << "StreamlineTracking - Min. tract length: " << m_Parameters->m_MinTractLength << "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_NumberOfSamples << std::endl; - std::cout << "StreamlineTracking - Max. sampling distance: " << m_SamplingDistance << "mm (" << m_SamplingDistance/m_MinVoxelSize << "*vox)" << std::endl; - std::cout << "StreamlineTracking - Deflection modifier: " << m_DeflectionMod << std::endl; + std::cout << "StreamlineTracking - Num. neighborhood samples: " << m_Parameters->m_NumSamples << std::endl; + std::cout << "StreamlineTracking - Max. sampling distance: " << m_Parameters->GetSamplingDistance() << "mm (" << m_Parameters->GetSamplingDistance()/m_MinVoxelSize << "*vox)" << std::endl; + std::cout << "StreamlineTracking - Deflection modifier: " << m_Parameters->m_DeflectionMod << std::endl; - std::cout << "StreamlineTracking - Use stop votes: " << m_UseStopVotes << std::endl; - std::cout << "StreamlineTracking - Only frontal samples: " << m_OnlyForwardSamples << 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_TrackingPriorWeight << ")" << std::endl; + 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_StepSize; - pos[1] += dir[1]*m_StepSize; - pos[2] += dir[2]*m_StepSize; + pos[0] += dir[0]*m_Parameters->GetStepSize(); + pos[1] += dir[1]*m_Parameters->GetStepSize(); + pos[2] += dir[2]*m_Parameters->GetStepSize(); } 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_InterpolateMasks, m_MaskInterpolator) && !mitk::imv::IsInsideMask(pos, m_InterpolateMasks, m_StopInterpolator)) + 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_NumberOfSamples); + 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_Random && m_RandomSampling) + 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_SamplingDistance))); + d *= static_cast(m_TrackingHandler->GetRandDouble(0, static_cast(m_Parameters->GetSamplingDistance()))); } else { d = probeVecs.at(i); float dot = dot_product(d, olddir); - if (m_UseStopVotes && dot>0.7f) + if (m_Parameters->m_StopVotes && dot>0.7f) { is_stop_voter = true; possible_stop_votes++; } - else if (m_OnlyForwardSamples && dot<0) + else if (m_Parameters->m_OnlyForwardSamples && dot<0) continue; - d *= m_SamplingDistance; + d *= m_Parameters->GetSamplingDistance(); } 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_InterpolateMasks, m_MaskInterpolator)) + 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_AvoidStop && olddir.magnitude()>0.5f) // out of white matter + 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_InterpolateMasks, m_MaskInterpolator)) + 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_DeflectionMod; // go into the direction of 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_IntroduceDirectionsFromPrior || valid)) + 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_TrackingPriorWeight) * direction + m_TrackingPriorWeight * prior; + direction = (1.0f-m_Parameters->m_Weight) * direction + m_Parameters->m_Weight * prior; direction.normalize(); } - else if (m_TrackingPriorAsMask) + 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< m_MaxLength/2; step++) + 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_InterpolateMasks, m_ExclusionInterpolator)) + 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_StepSize; + tractLength += m_Parameters->GetStepSize(); - if (m_LoopCheck>=0 && CheckCurvature(container, front)>m_LoopCheck) + if (m_Parameters->GetLoopCheckDeg()>=0 && CheckCurvature(container, front)>m_Parameters->GetLoopCheckDeg()) return tractLength; - if (tractLength>m_MaxTractLength) + if (tractLength>m_Parameters->m_MaxTractLength) return tractLength; - if (m_DemoMode && !m_UseOutputProbabilityMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? + 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_NumPreviousDirections) + 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_MinVoxelSize*4, m_StepSize*8); + float m_Distance = std::max(m_MinVoxelSize*4, m_Parameters->GetStepSize()*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_StepSize; + dist += m_Parameters->GetStepSize(); 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_StepSize; + dist += m_Parameters->GetStepSize(); 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_InterpolateMasks, m_MaskInterpolator) ) + if ( mitk::imv::IsInsideMask(worldPos, m_Parameters->m_InterpolateRoiImages, m_MaskInterpolator) ) { m_SeedPoints.push_back(worldPos); - for (int s = 1; s < m_SeedsPerVoxel; s++) + 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_Random) + 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; #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_MaxNumTracts>0) - std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << "/" << m_MaxNumTracts << '\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; trialsm_TrialsPerSeed; ++trials) { FiberType fib; DirectionContainer direction_container; float tractLength = 0; unsigned long counter = 0; // get starting direction vnl_vector_fixed 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_InterpolateMasks, m_ExclusionInterpolator)) + 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_MinTractLength && counter>=2 && !exclude) + if (tractLength>=m_Parameters->m_MinTractLength && counter>=2 && !exclude) { #pragma omp critical if ( IsValidFiber(&fib) ) { if (!m_StopTracking) { - if (!m_UseOutputProbabilityMap) + if (!m_Parameters->m_OutputProbMap) m_Tractogram.push_back(fib); else FiberToProbmap(&fib); m_CurrentTracts++; success = true; } - if (m_MaxNumTracts > 0 && m_CurrentTracts>=static_cast(m_MaxNumTracts)) + 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_TrackingHandler->GetMode()!=mitk::TrackingDataHandler::PROBABILISTIC) + 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_EndpointConstraint==EndpointConstraints::NONE) + if (m_Parameters->m_EpConstraints==EndpointConstraints::NONE) { return true; } - else if (m_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET) + else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_TARGET) { if (m_TargetImageSet) { - if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) - && mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) + 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_EndpointConstraint==EndpointConstraints::EPS_IN_TARGET_LABELDIFF) + 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_EndpointConstraint==EndpointConstraints::EPS_IN_SEED_AND_TARGET) + else if (m_Parameters->m_EpConstraints==EndpointConstraints::EPS_IN_SEED_AND_TARGET) { if (m_TargetImageSet && m_SeedImageSet) { - if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_SeedInterpolator) - && mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) + 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_InterpolateMasks, m_SeedInterpolator) - && mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) ) + 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_EndpointConstraint==EndpointConstraints::MIN_ONE_EP_IN_TARGET) + else if (m_Parameters->m_EpConstraints==EndpointConstraints::MIN_ONE_EP_IN_TARGET) { if (m_TargetImageSet) { - if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) - || mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) + 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_EndpointConstraint==EndpointConstraints::ONE_EP_IN_TARGET) + else if (m_Parameters->m_EpConstraints==EndpointConstraints::ONE_EP_IN_TARGET) { if (m_TargetImageSet) { - if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) - && !mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) + 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_InterpolateMasks, m_TargetInterpolator) - && mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) ) + 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_EndpointConstraint==EndpointConstraints::NO_EP_IN_TARGET) + else if (m_Parameters->m_EpConstraints==EndpointConstraints::NO_EP_IN_TARGET) { if (m_TargetImageSet) { - if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMasks, m_TargetInterpolator) - || mitk::imv::IsInsideMask(fib->back(), m_InterpolateMasks, m_TargetInterpolator) ) + 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_UseOutputProbabilityMap) + 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_TrackingHandler->GetMode()==mitk::TrackingDataHandler::DETERMINISTIC && dynamic_cast(m_TrackingHandler) && !m_TrackingHandler->GetInterpolate()) + if ( m_Parameters->m_Mode==MODE::DETERMINISTIC && dynamic_cast(m_TrackingHandler) && !m_Parameters->m_InterpolateTractographyData) { algo_code_value = "sup181_ee04"; algo_code_meaning = "FACT"; } - else if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::DETERMINISTIC) + else if (m_Parameters->m_Mode==MODE::DETERMINISTIC) { algo_code_value = "sup181_ee01"; algo_code_meaning = "Deterministic"; } - else if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::PROBABILISTIC) + 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_code.value", mitk::StringProperty::New(algo_code_value)); fib->SetProperty("DICOM.algo_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)); } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h index 3817102215..3875db4794 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h @@ -1,258 +1,198 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkMLBSTrackingFilter_h_ #define __itkMLBSTrackingFilter_h_ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include namespace itk{ /** * \brief Performs streamline tracking on the input image. Depending on the tracking handler this can be a tensor, peak or machine learning based tracking. */ class MITKFIBERTRACKING_EXPORT StreamlineTrackingFilter : public ProcessObject { public: - enum EndpointConstraints { - NONE, ///< No constraints on endpoint locations - EPS_IN_TARGET, ///< Both EPs are required to be located in the target image - EPS_IN_TARGET_LABELDIFF, ///< Both EPs are required to be located in the target image and the image values at the respective position needs to be distinct - EPS_IN_SEED_AND_TARGET, ///< One EP is required to be located in the seed image and one in the target image - MIN_ONE_EP_IN_TARGET, ///< At least one EP is required to be located in the target image - ONE_EP_IN_TARGET, ///< Exactly one EP is required to be located in the target image - NO_EP_IN_TARGET ///< No EP is allowed to be located in the target image - }; - typedef StreamlineTrackingFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ProcessObject Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(MLBSTrackingFilter, ImageToImageFilter) typedef itk::Image ItkUcharImgType; typedef itk::Image ItkUintImgType; typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkFloatImgType; typedef vtkSmartPointer< vtkPolyData > PolyDataType; + typedef mitk::StreamlineTractographyParameters::EndpointConstraints EndpointConstraints; + typedef mitk::StreamlineTractographyParameters::MODE MODE; typedef std::deque< vnl_vector_fixed > DirectionContainer; typedef std::deque< itk::Point > FiberType; typedef std::vector< FiberType > BundleType; volatile bool m_PauseTracking; bool m_AbortTracking; bool m_BuildFibersFinished; int m_BuildFibersReady; volatile bool m_Stop; mitk::PointSet::Pointer m_SamplingPointset; mitk::PointSet::Pointer m_StopVotePointset; mitk::PointSet::Pointer m_AlternativePointset; - void SetStepSize(float v) ///< Integration step size in voxels, default is 0.5 * voxel - { m_StepSizeVox = v; } - void SetAngularThreshold(float v) ///< Angular threshold per step (in degree), default is 90deg x stepsize - { m_AngularThresholdDeg = v; } - void SetSamplingDistance(float v) ///< Maximum distance of sampling points in voxels, default is 0.25 * voxel - { m_SamplingDistanceVox = v; } - void SetDicomProperties(mitk::FiberBundle::Pointer fib); itkGetMacro( OutputProbabilityMap, ItkDoubleImgType::Pointer) ///< Output probability map itkGetMacro( FiberPolyData, PolyDataType ) ///< Output fibers - itkGetMacro( UseOutputProbabilityMap, bool) itkGetMacro( MinVoxelSize, float) - itkGetMacro( EndpointConstraint, EndpointConstraints) itkSetMacro( SeedImage, ItkFloatImgType::Pointer) ///< Seeds are only placed inside of this mask. itkSetMacro( MaskImage, ItkFloatImgType::Pointer) ///< Tracking is only performed inside of this mask image. itkSetMacro( ExclusionRegions, ItkFloatImgType::Pointer)///< Fibers passing any of the ROIs in this image are discarded. - itkSetMacro( SeedsPerVoxel, int) ///< One seed placed in the center of each voxel or multiple seeds randomly placed inside each voxel. - itkSetMacro( MinTractLength, float ) ///< Shorter tracts are discarded. - itkSetMacro( MaxTractLength, float ) ///< Streamline progression stops if tract is longer than specified. - itkSetMacro( EndpointConstraint, EndpointConstraints) ///< Determines what fibers are accepted based on their endpoint location - - itkSetMacro( UseStopVotes, bool ) ///< Frontal sampling points can vote for stopping the streamline even if the remaining sampling points keep pushing - itkSetMacro( OnlyForwardSamples, bool ) ///< Don't use sampling points behind the current position in progression direction - itkSetMacro( DeflectionMod, float ) ///< Deflection distance modifier + itkSetMacro( StoppingRegions, ItkFloatImgType::Pointer) ///< Streamlines entering a stopping region will stop immediately itkSetMacro( TargetRegions, ItkFloatImgType::Pointer) ///< Only streamline starting and ending in this mask are retained itkSetMacro( DemoMode, bool ) - itkSetMacro( NumberOfSamples, unsigned int ) ///< Number of neighborhood sampling points - itkSetMacro( LoopCheck, float ) ///< Checks fiber curvature (angular deviation across 5mm) is larger than 30°. If yes, the streamline progression is stopped. - itkSetMacro( AvoidStop, bool ) ///< Use additional sampling points to avoid premature streamline termination - itkSetMacro( RandomSampling, bool ) ///< If true, the sampling points are distributed randomly around the current position, not sphericall in the specified sampling distance. - itkSetMacro( NumPreviousDirections, unsigned int ) ///< How many "old" steps do we want to consider in our decision where to go next? - itkSetMacro( MaxNumTracts, int ) ///< Tracking is stopped if the maximum number of tracts is exceeded - itkSetMacro( Random, bool ) ///< If true, seedpoints are shuffled randomly before tracking itkSetMacro( Verbose, bool ) ///< If true, output tracking progress (might be slower) - itkSetMacro( UseOutputProbabilityMap, bool) ///< If true, no tractogram but a probability map is created as output. itkSetMacro( StopTracking, bool ) - itkSetMacro( InterpolateMasks, bool ) - itkSetMacro( TrialsPerSeed, unsigned int ) ///< When using probabilistic tractography, each seed point is used N times until a valid streamline that is compliant with all thresholds etc. is found - itkSetMacro( TrackingPriorWeight, float) ///< Weight between prior and data [0-1]. One mean tracking only on the prior peaks, zero only on the data. - itkSetMacro( TrackingPriorAsMask, bool) ///< If true, data directions in voxels where prior directions are invalid are set to zero - itkSetMacro( IntroduceDirectionsFromPrior, bool) ///< If false, prior voxels with invalid data voxel are ignored ///< Use manually defined points in physical space as seed points instead of seed image void SetSeedPoints( const std::vector< itk::Point >& sP) { m_SeedPoints = sP; } void SetTrackingHandler( mitk::TrackingDataHandler* h ) ///< { m_TrackingHandler = h; } void Update() override{ this->GenerateData(); } std::string GetStatusText(); void SetTrackingPriorHandler(mitk::TrackingDataHandler *TrackingPriorHandler); + std::shared_ptr< mitk::StreamlineTractographyParameters > GetParameters() const; + void SetParameters(std::shared_ptr< mitk::StreamlineTractographyParameters > Parameters); + protected: void GenerateData() override; StreamlineTrackingFilter(); ~StreamlineTrackingFilter() override {} bool IsValidFiber(FiberType* fib); ///< Check endpoints void FiberToProbmap(FiberType* fib); void GetSeedPointsFromSeedImage(); void CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir); ///< Calculate next integration step. float FollowStreamline(itk::Point start_pos, vnl_vector_fixed dir, FiberType* fib, DirectionContainer* container, float tractLength, bool front, bool& exclude); ///< Start streamline in one direction. vnl_vector_fixed GetNewDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< Determine new direction by sample voting at the current position taking the last progression direction into account. std::vector< vnl_vector_fixed > CreateDirections(unsigned int NPoints); void BeforeTracking(); void AfterTracking(); PolyDataType m_FiberPolyData; vtkSmartPointer m_Points; vtkSmartPointer m_Cells; BundleType m_Tractogram; BundleType m_GmStubs; ItkFloatImgType::Pointer m_StoppingRegions; ItkFloatImgType::Pointer m_TargetRegions; ItkFloatImgType::Pointer m_SeedImage; ItkFloatImgType::Pointer m_MaskImage; ItkFloatImgType::Pointer m_ExclusionRegions; ItkDoubleImgType::Pointer m_OutputProbabilityMap; float m_MinVoxelSize; - float m_AngularThresholdDeg; - float m_StepSizeVox; - float m_SamplingDistanceVox; - float m_AngularThreshold; - float m_StepSize; - int m_MaxLength; - float m_MinTractLength; - float m_MaxTractLength; - int m_SeedsPerVoxel; - - bool m_AvoidStop; - bool m_RandomSampling; - float m_SamplingDistance; - float m_DeflectionMod; - bool m_OnlyForwardSamples; - bool m_UseStopVotes; - unsigned int m_NumberOfSamples; - - unsigned int m_NumPreviousDirections; - int m_MaxNumTracts; bool m_Verbose; - float m_LoopCheck; bool m_DemoMode; - bool m_Random; - bool m_UseOutputProbabilityMap; std::vector< itk::Point > m_SeedPoints; unsigned int m_CurrentTracts; unsigned int m_Progress; bool m_StopTracking; - bool m_InterpolateMasks; - unsigned int m_TrialsPerSeed; - EndpointConstraints m_EndpointConstraint; void BuildFibers(bool check); float CheckCurvature(DirectionContainer *fib, bool front); // decision forest mitk::TrackingDataHandler* m_TrackingHandler; std::vector< PolyDataType > m_PolyDataContainer; std::chrono::time_point m_StartTime; std::chrono::time_point m_EndTime; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_MaskInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_StopInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_TargetInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_SeedInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_ExclusionInterpolator; bool m_SeedImageSet; bool m_TargetImageSet; + std::shared_ptr< mitk::StreamlineTractographyParameters > m_Parameters; + // Directional priors - bool m_IntroduceDirectionsFromPrior; - bool m_TrackingPriorAsMask; - float m_TrackingPriorWeight; mitk::TrackingDataHandler* m_TrackingPriorHandler; private: }; } //#ifndef ITK_MANUAL_INSTANTIATION //#include "itkMLBSTrackingFilter.cpp" //#endif #endif //__itkMLBSTrackingFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkStreamlineTractographyParameters.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkStreamlineTractographyParameters.cpp new file mode 100644 index 0000000000..ca0ffe25e1 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkStreamlineTractographyParameters.cpp @@ -0,0 +1,348 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +mitk::StreamlineTractographyParameters::StreamlineTractographyParameters() +{ + SetStepSizeVox(0.5); + SetAngularThresholdDeg(-1); + SetLoopCheckDeg(-1); +} + +mitk::StreamlineTractographyParameters::~StreamlineTractographyParameters() +{ + +} + +void mitk::StreamlineTractographyParameters::SaveParameters(std::string filename) +{ + if(filename.empty()) + return; + if(".stp"!=filename.substr(filename.size()-4, 4)) + filename += ".stp"; + + const std::string& locale = "C"; + const std::string& currLocale = setlocale( LC_ALL, nullptr ); + + if ( locale.compare(currLocale)!=0 ) + { + try + { + setlocale(LC_ALL, locale.c_str()); + } + catch(...) + { + MITK_INFO << "Could not set locale " << locale; + } + } + + boost::property_tree::ptree parameters; + + parameters.put("seeding.seeds_per_voxel", m_SeedsPerVoxel); + parameters.put("seeding.trials_per_seed", m_TrialsPerSeed); + parameters.put("seeding.max_num_fibers", m_MaxNumFibers); + parameters.put("seeding.interactive_radius", m_InteractiveRadius); + parameters.put("seeding.num_interactive_seeds", m_NumInteractiveSeeds); + parameters.put("seeding.enable_interactive", m_EnableInteractive); + + parameters.put("roi_constraints.ep_constraints", m_EpConstraints); + + parameters.put("tractography.mode", m_Mode); + parameters.put("tractography.sharpen_odfs", m_SharpenOdfs); + parameters.put("tractography.cutoff", m_Cutoff); + parameters.put("tractography.odf_cutoff", m_OdfCutoff); + parameters.put("tractography.step_size", m_StepSizeVox); + parameters.put("tractography.min_tract_length", m_MinTractLength); + parameters.put("tractography.angluar_threshold", m_AngularThresholdDeg); + parameters.put("tractography.loop_check", m_LoopCheckDeg); + parameters.put("tractography.f", m_F); + parameters.put("tractography.g", m_G); + parameters.put("tractography.fix_seed", m_FixRandomSeed); + + parameters.put("prior.weight", m_Weight); + parameters.put("prior.restrict_to_prior", m_RestrictToPrior); + parameters.put("prior.new_directions", m_NewDirectionsFromPrior); + + parameters.put("nsampling.num_samples", m_NumSamples); + parameters.put("nsampling.sampling_distance", m_SamplingDistance); + parameters.put("nsampling.only_frontal", m_OnlyForwardSamples); + parameters.put("nsampling.stop_votes", m_StopVotes); + + parameters.put("data_handling.flip_x", m_FlipX); + parameters.put("data_handling.flip_y", m_FlipY); + parameters.put("data_handling.flip_z", m_FlipZ); + parameters.put("data_handling.interpolate_tracto_data", m_InterpolateTractographyData); + parameters.put("data_handling.interpolate_roi_images", m_InterpolateRoiImages); + + parameters.put("output.compress", m_CompressFibers); + parameters.put("output.compression", m_Compression); + parameters.put("output.prob_map", m_OutputProbMap); + + boost::property_tree::json_parser::write_json(filename, parameters, std::locale(), true); + +// try{ +// itk::ImageFileWriter::Pointer writer = itk::ImageFileWriter::New(); +// writer->SetFileName(filename+"_FMAP.nii.gz"); +// writer->SetInput(m_SignalGen.m_FrequencyMap); +// writer->Update(); +// } +// catch(...) +// { +// MITK_INFO << "No frequency map saved."; +// } + + setlocale(LC_ALL, currLocale.c_str()); +} + + + +template< class ParameterType > +ParameterType mitk::StreamlineTractographyParameters::ReadVal(boost::property_tree::ptree::value_type const& v, std::string tag, ParameterType defaultValue, bool essential) +{ + try + { + return v.second.get(tag); + } + catch (...) + { + if (essential) + { + mitkThrow() << "Parameter file corrupted. Essential tag is missing: '" << tag << "'"; + } + MITK_INFO << "Tag '" << tag << "' not found. Using default value '" << defaultValue << "'."; + return defaultValue; + } +} + +void mitk::StreamlineTractographyParameters::LoadParameters(std::string filename) +{ + if(filename.empty()) { return; } + + const std::string& locale = "C"; + const std::string& currLocale = setlocale( LC_ALL, nullptr ); + + if ( locale.compare(currLocale)!=0 ) + { + try + { + setlocale(LC_ALL, locale.c_str()); + } + catch(...) + { + MITK_INFO << "Could not set locale " << locale; + } + } + + + boost::property_tree::ptree parameterTree; + boost::property_tree::json_parser::read_json( filename, parameterTree ); + + BOOST_FOREACH( boost::property_tree::ptree::value_type const& v1, parameterTree.get_child("fiberfox") ) + { + if( v1.first == "seeding" ) + { + m_SeedsPerVoxel = ReadVal(v1,"seeds_per_voxel", m_SeedsPerVoxel); + m_TrialsPerSeed = ReadVal(v1,"trials_per_seed", m_TrialsPerSeed); + m_MaxNumFibers = ReadVal(v1,"max_num_fibers", m_MaxNumFibers); + m_InteractiveRadius = ReadVal(v1,"interactive_radius", m_InteractiveRadius); + m_NumInteractiveSeeds = ReadVal(v1,"num_interactive_seeds", m_NumInteractiveSeeds); + m_EnableInteractive = ReadVal(v1,"enable_interactive", m_EnableInteractive); + } + else if( v1.first == "roi_constraints" ) + { + switch( ReadVal(v1,"ep_constraints", 0) ) + { + default: + m_EpConstraints = EndpointConstraints::NONE; + break; + case 1: + m_EpConstraints = EndpointConstraints::EPS_IN_TARGET; + break; + case 2: + m_EpConstraints = EndpointConstraints::EPS_IN_TARGET_LABELDIFF; + break; + case 3: + m_EpConstraints = EndpointConstraints::EPS_IN_SEED_AND_TARGET; + break; + case 4: + m_EpConstraints = EndpointConstraints::MIN_ONE_EP_IN_TARGET; + break; + case 5: + m_EpConstraints = EndpointConstraints::ONE_EP_IN_TARGET; + break; + case 6: + m_EpConstraints = EndpointConstraints::NO_EP_IN_TARGET; + break; + } + } + else if( v1.first == "tractography" ) + { + if(ReadVal(v1,"mode", 0) == 0) + m_Mode = MODE::DETERMINISTIC; + else + m_Mode = MODE::PROBABILISTIC; + + m_SharpenOdfs = ReadVal(v1,"sharpen_odfs", m_SharpenOdfs); + m_Cutoff = ReadVal(v1,"cutoff", m_Cutoff); + m_OdfCutoff = ReadVal(v1,"odf_cutoff", m_OdfCutoff); + SetStepSizeVox(ReadVal(v1,"step_size", m_StepSizeVox)); + m_MinTractLength = ReadVal(v1,"min_tract_length", m_MinTractLength); + SetAngularThresholdDeg(ReadVal(v1,"angluar_threshold", m_AngularThresholdDeg)); + SetLoopCheckDeg(ReadVal(v1,"loop_check", m_LoopCheckDeg)); + m_F = ReadVal(v1,"f", m_F); + m_G = ReadVal(v1,"g", m_G); + m_FixRandomSeed = ReadVal(v1,"fix_seed", m_FixRandomSeed); + } + else if( v1.first == "prior" ) + { + m_Weight = ReadVal(v1,"weight", m_Weight); + m_RestrictToPrior = ReadVal(v1,"restrict_to_prior", m_RestrictToPrior); + m_NewDirectionsFromPrior = ReadVal(v1,"new_directions", m_NewDirectionsFromPrior); + } + else if( v1.first == "nsampling" ) + { + m_NumSamples = ReadVal(v1,"num_samples", m_NumSamples); + m_SamplingDistance = ReadVal(v1,"sampling_distance", m_SamplingDistance); + m_OnlyForwardSamples = ReadVal(v1,"only_frontal", m_OnlyForwardSamples); + m_StopVotes = ReadVal(v1,"stop_votes", m_StopVotes); + } + else if( v1.first == "data_handling" ) + { + m_FlipX = ReadVal(v1,"flip_x", m_FlipX); + m_FlipY = ReadVal(v1,"flip_y", m_FlipY); + m_FlipZ = ReadVal(v1,"flip_z", m_FlipZ); + m_InterpolateTractographyData = ReadVal(v1,"interpolate_tracto_data", m_InterpolateTractographyData); + m_InterpolateRoiImages = ReadVal(v1,"interpolate_roi_images", m_InterpolateRoiImages); + } + else if( v1.first == "output" ) + { + m_CompressFibers = ReadVal(v1,"compress", m_CompressFibers); + m_Compression = ReadVal(v1,"compression", m_Compression); + m_OutputProbMap = ReadVal(v1,"prob_map", m_OutputProbMap); + } + } + +// try +// { +// itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); +// reader->SetFileName(filename+"_FMAP.nrrd"); +// if ( itksys::SystemTools::FileExists(filename+"_FMAP.nii.gz") ) +// reader->SetFileName(filename+"_FMAP.nii.gz"); +// else if ( itksys::SystemTools::FileExists(filename+"_FMAP.nii") ) +// reader->SetFileName(filename+"_FMAP.nii"); +// else +// reader->SetFileName(filename+"_FMAP.nrrd"); +// reader->Update(); +// m_SignalGen.m_FrequencyMap = reader->GetOutput(); +// MITK_INFO << "Frequency map loaded."; +// } +// catch(...) +// { +// MITK_INFO << "No frequency map found."; +// } + + setlocale(LC_ALL, currLocale.c_str()); +} + +float mitk::StreamlineTractographyParameters::GetSamplingDistance() const +{ + return m_SamplingDistance; +} + +void mitk::StreamlineTractographyParameters::SetSamplingDistance(float SamplingDistance) +{ + m_SamplingDistance = SamplingDistance; + AutoAdjust(); +} + +void mitk::StreamlineTractographyParameters::AutoAdjust() +{ + if (m_StepSizeVox(mitk::eps)) + m_StepSizeVox = 0.5; + m_StepSize = m_StepSizeVox*m_MinVoxelSize; + + if (m_AngularThresholdDeg<0) + { + if (m_StepSize/m_MinVoxelSize<=0.966f) // minimum 15° for automatic estimation + m_AngularThreshold = static_cast(std::cos( 0.5 * itk::Math::pi * static_cast(m_StepSize/m_MinVoxelSize) )); + else + m_AngularThreshold = static_cast(std::cos( 0.5 * itk::Math::pi * 0.966 )); + + m_AngularThresholdDeg = std::acos(m_AngularThreshold)*180.0/itk::Math::pi; + } + else + m_AngularThreshold = static_cast(std::cos( static_cast(m_AngularThresholdDeg)*itk::Math::pi/180.0 )); + + if (m_SamplingDistance(mitk::eps)) + m_SamplingDistance = m_MinVoxelSize*0.25f; +} + +float mitk::StreamlineTractographyParameters::GetStepSize() const +{ + return m_StepSize; +} + +void mitk::StreamlineTractographyParameters::SetMinVoxelSize(float MinVoxelSize) +{ + m_MinVoxelSize = MinVoxelSize; + AutoAdjust(); +} + +float mitk::StreamlineTractographyParameters::GetAngularThreshold() const +{ + return m_AngularThreshold; +} + +float mitk::StreamlineTractographyParameters::GetStepSizeVox() const +{ + return m_StepSizeVox; +} + +void mitk::StreamlineTractographyParameters::SetStepSizeVox(float StepSizeVox) +{ + m_StepSizeVox = StepSizeVox; + AutoAdjust(); +} + +float mitk::StreamlineTractographyParameters::GetLoopCheckDeg() const +{ + return m_LoopCheckDeg; +} + +void mitk::StreamlineTractographyParameters::SetLoopCheckDeg(float LoopCheckDeg) +{ + m_LoopCheckDeg = LoopCheckDeg; +} + +float mitk::StreamlineTractographyParameters::GetAngularThresholdDeg() const +{ + return m_AngularThresholdDeg; +} + +void mitk::StreamlineTractographyParameters::SetAngularThresholdDeg(float AngularThresholdDeg) +{ + m_AngularThresholdDeg = AngularThresholdDeg; + AutoAdjust(); +} diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkStreamlineTractographyParameters.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkStreamlineTractographyParameters.h new file mode 100644 index 0000000000..8ff9fdd9b5 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkStreamlineTractographyParameters.h @@ -0,0 +1,163 @@ +#pragma once + +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include +#include +#include +#include +#include + +namespace mitk +{ + /** + * \brief Datastructure to manage streamline tractography parameters. + * + */ + class MITKFIBERTRACKING_EXPORT StreamlineTractographyParameters + { + + public: + + enum EndpointConstraints { + NONE, ///< No constraints on endpoint locations + EPS_IN_TARGET, ///< Both EPs are required to be located in the target image + EPS_IN_TARGET_LABELDIFF, ///< Both EPs are required to be located in the target image and the image values at the respective position needs to be distinct + EPS_IN_SEED_AND_TARGET, ///< One EP is required to be located in the seed image and one in the target image + MIN_ONE_EP_IN_TARGET, ///< At least one EP is required to be located in the target image + ONE_EP_IN_TARGET, ///< Exactly one EP is required to be located in the target image + NO_EP_IN_TARGET ///< No EP is allowed to be located in the target image + }; + + enum MODE { + DETERMINISTIC, + PROBABILISTIC + }; + + typedef itk::Image ItkFloatImgType; + typedef itk::Image ItkDoubleImgType; + typedef itk::Image ItkUcharImgType; + + StreamlineTractographyParameters(); + StreamlineTractographyParameters(const StreamlineTractographyParameters ¶ms) = default; + ~StreamlineTractographyParameters(); + + void SaveParameters(std::string filename); ///< Save image generation parameters to .stp file. + void LoadParameters(std::string filename); ///< Load image generation parameters from .stp file. + + template< class ParameterType > + ParameterType ReadVal(boost::property_tree::ptree::value_type const& v, std::string tag, ParameterType defaultValue, bool essential=false); + + // seeding + unsigned int m_SeedsPerVoxel = 1; + unsigned int m_TrialsPerSeed = 10; + int m_MaxNumFibers = -1; + // - seed image + + // interactive + float m_InteractiveRadius = 2; + unsigned int m_NumInteractiveSeeds = 50; + bool m_EnableInteractive = false; + + // ROI constraints + EndpointConstraints m_EpConstraints; + // - mask image + // - stop image + // - exclusion image + // - target image + + // tractography + MODE m_Mode; + bool m_SharpenOdfs = false; + float m_Cutoff = 0.1; + // - fa/gfa image + float m_OdfCutoff = 0.00025; + float m_MinTractLength = 20; + float m_MaxTractLength = 400; + float m_F = 1; + float m_G = 0; + bool m_FixRandomSeed = false; + unsigned int m_NumPreviousDirections = 1; + + // prior + // - peak image + float m_Weight = 0.5; + bool m_RestrictToPrior = true; + bool m_NewDirectionsFromPrior = true; + + // neighborhood sampling + unsigned int m_NumSamples = 0; + bool m_OnlyForwardSamples = false; + bool m_StopVotes = false; + bool m_AvoidStop = true; + bool m_RandomSampling = false; + float m_DeflectionMod = 1.0; + + // data handling + bool m_FlipX = false; + bool m_FlipY = false; + bool m_FlipZ = false; + bool m_InterpolateTractographyData = true; + bool m_InterpolateRoiImages; + bool m_ApplyDirectionMatrix = false; + + // output and postprocessing + bool m_CompressFibers = true; + float m_Compression = 0.1; + bool m_OutputProbMap = false; + + float GetAngularThresholdDeg() const; + void SetAngularThresholdDeg(float AngularThresholdDeg); + + float GetLoopCheckDeg() const; + void SetLoopCheckDeg(float LoopCheckDeg); + + float GetStepSizeVox() const; + void SetStepSizeVox(float StepSizeVox); + + void SetMinVoxelSize(float MinVoxelSize); + + float GetAngularThreshold() const; + + float GetSamplingDistance() const; + void SetSamplingDistance(float SamplingDistance); + + float GetStepSize() const; + + private: + + void AutoAdjust(); + + float m_SamplingDistance = -1; + + float m_AngularThresholdDeg; + float m_LoopCheckDeg; + float m_StepSizeVox; + + float m_AngularThreshold; + float m_LoopCheck; + float m_StepSize; // mm + + float m_MinVoxelSize = 1.0; + + +// float m_AngularThreshold = 0; // in deg +// float m_LoopCheck = 0; // in deg +// float m_StepSize = 0; + }; +} + diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp index 2e94c0b96b..a1c6b0b1b3 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp @@ -1,109 +1,110 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include #include #include #include #include #include #include #include #include #include #include +#include #include "mitkTestFixture.h" class mitkMachineLearningTrackingTestSuite : public mitk::TestFixture { - CPPUNIT_TEST_SUITE(mitkMachineLearningTrackingTestSuite); - MITK_TEST(Track1); - CPPUNIT_TEST_SUITE_END(); + CPPUNIT_TEST_SUITE(mitkMachineLearningTrackingTestSuite); + MITK_TEST(Track1); + CPPUNIT_TEST_SUITE_END(); - typedef itk::Image ItkFloatImgType; + typedef itk::Image ItkFloatImgType; private: - /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ - mitk::FiberBundle::Pointer ref; - mitk::TrackingHandlerRandomForest<6, 100>* tfh; - mitk::Image::Pointer dwi; - ItkFloatImgType::Pointer seed; + /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ + mitk::FiberBundle::Pointer ref; + mitk::TrackingHandlerRandomForest<6, 100>* tfh; + mitk::Image::Pointer dwi; + ItkFloatImgType::Pointer seed; public: - void setUp() override - { - ref = nullptr; - tfh = new mitk::TrackingHandlerRandomForest<6,100>(); - - std::vector fibInfile = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/ReferenceTracts.fib")); - mitk::BaseData::Pointer baseData = fibInfile.at(0); - ref = dynamic_cast(baseData.GetPointer()); - - dwi = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/DiffusionImage.dwi")); - - - mitk::Image::Pointer img = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/seed.nrrd")); - seed = ItkFloatImgType::New(); - mitk::CastToItkImage(img, seed); - - mitk::TractographyForest::Pointer forest = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/forest.rf")); - - tfh->SetForest(forest); - tfh->AddDwi(dwi); - } - - void tearDown() override - { - delete tfh; - ref = nullptr; - } - - void Track1() - { - omp_set_num_threads(1); - typedef itk::StreamlineTrackingFilter TrackerType; - TrackerType::Pointer tracker = TrackerType::New(); - tracker->SetDemoMode(false); - tracker->SetInterpolateMasks(false); - tracker->SetSeedImage(seed); - tracker->SetSeedsPerVoxel(1); - tracker->SetStepSize(-1); - tracker->SetAngularThreshold(45); - tracker->SetMinTractLength(20); - tracker->SetMaxTractLength(400); - tracker->SetTrackingHandler(tfh); - tracker->SetAvoidStop(true); - tracker->SetSamplingDistance(0.5); - tracker->SetRandomSampling(false); - tracker->Update(); - vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); - mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); - - //MITK_INFO << mitk::IOUtil::GetTempPath() << "ReferenceTracts.fib"; - if (!ref->Equals(outFib)) - mitk::IOUtil::Save(outFib, mitk::IOUtil::GetTempPath()+"ML_Track1.fib"); - - CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(outFib)); - } + void setUp() override + { + ref = nullptr; + tfh = new mitk::TrackingHandlerRandomForest<6,100>(); + + std::vector fibInfile = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/ReferenceTracts.fib")); + mitk::BaseData::Pointer baseData = fibInfile.at(0); + ref = dynamic_cast(baseData.GetPointer()); + + dwi = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/DiffusionImage.dwi")); + + + mitk::Image::Pointer img = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/seed.nrrd")); + seed = ItkFloatImgType::New(); + mitk::CastToItkImage(img, seed); + + mitk::TractographyForest::Pointer forest = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/forest.rf")); + + tfh->SetForest(forest); + tfh->AddDwi(dwi); + } + + void tearDown() override + { + delete tfh; + ref = nullptr; + } + + void Track1() + { + std::shared_ptr params = std::make_shared(); + params->m_RandomSampling = false; + params->SetSamplingDistance(0.5); + params->m_AvoidStop = true; + params->SetAngularThresholdDeg(45); + params->SetStepSizeVox(-1); + params->m_SeedsPerVoxel = 1; + params->m_InterpolateRoiImages = false; + omp_set_num_threads(1); + typedef itk::StreamlineTrackingFilter TrackerType; + TrackerType::Pointer tracker = TrackerType::New(); + tracker->SetDemoMode(false); + tracker->SetSeedImage(seed); + tracker->SetTrackingHandler(tfh); + tracker->SetParameters(params); + tracker->Update(); + vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); + mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); + + //MITK_INFO << mitk::IOUtil::GetTempPath() << "ReferenceTracts.fib"; + if (!ref->Equals(outFib)) + mitk::IOUtil::Save(outFib, mitk::IOUtil::GetTempPath()+"ML_Track1.fib"); + + CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(outFib)); + } }; MITK_TEST_SUITE_REGISTRATION(mitkMachineLearningTracking) diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp index 508a0df61d..edb22aadd7 100755 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp @@ -1,426 +1,438 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include class mitkStreamlineTractographyTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkStreamlineTractographyTestSuite); MITK_TEST(Test_Peak1); MITK_TEST(Test_Peak2); MITK_TEST(Test_Tensor1); MITK_TEST(Test_Tensor2); MITK_TEST(Test_Tensor3); MITK_TEST(Test_Odf1); MITK_TEST(Test_Odf2); MITK_TEST(Test_Odf3); MITK_TEST(Test_Odf4); MITK_TEST(Test_Odf5); MITK_TEST(Test_Odf6); CPPUNIT_TEST_SUITE_END(); typedef itk::VectorImage< short, 3> ItkDwiType; private: public: /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ typedef itk::Image ItkFloatImgType; mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itk_odf_image; mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itk_tensor_image; mitk::TrackingHandlerPeaks::PeakImgType::Pointer itk_peak_image; ItkFloatImgType::Pointer itk_seed_image; ItkFloatImgType::Pointer itk_mask_image; ItkFloatImgType::Pointer itk_gfa_image; float gfa_threshold; float odf_threshold; float peak_threshold; + std::shared_ptr params; itk::StreamlineTrackingFilter::Pointer tracker; void setUp() override { omp_set_num_threads(1); gfa_threshold = 0.2f; odf_threshold = 0.1f; peak_threshold = 0.1f; mitk::Image::Pointer odf_image = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/qball_image.qbi")); mitk::Image::Pointer tensor_image = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/tensor_image.dti")); mitk::Image::Pointer peak_image = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/qball_peak_image.nii.gz")); mitk::Image::Pointer seed_image = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/seed_image.nii.gz")); mitk::Image::Pointer mask_image = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/mask_image.nii.gz")); mitk::Image::Pointer gfa_image = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/gfa_image.nii.gz")); + params = std::make_shared(); + params->m_FixRandomSeed = true; + params->m_InterpolateRoiImages = false; + { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(peak_image); caster->Update(); itk_peak_image = caster->GetOutput(); } { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(tensor_image); caster->Update(); itk_tensor_image = caster->GetOutput(); } { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(odf_image); caster->Update(); itk_odf_image = caster->GetOutput(); } itk_gfa_image = ItkFloatImgType::New(); mitk::CastToItkImage(gfa_image, itk_gfa_image); itk_seed_image = ItkFloatImgType::New(); mitk::CastToItkImage(seed_image, itk_seed_image); itk_mask_image = ItkFloatImgType::New(); mitk::CastToItkImage(mask_image, itk_mask_image); } mitk::FiberBundle::Pointer LoadReferenceFib(std::string filename) { mitk::FiberBundle::Pointer fib = nullptr; if (itksys::SystemTools::FileExists(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename))) { mitk::BaseData::Pointer baseData = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename)).at(0); fib = dynamic_cast(baseData.GetPointer()); } return fib; } mitk::Image::Pointer LoadReferenceImage(std::string filename) { mitk::Image::Pointer img = nullptr; if (itksys::SystemTools::FileExists(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename))) { img = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename)); } return img; } void SetupTracker(mitk::TrackingDataHandler* handler) { + tracker = itk::StreamlineTrackingFilter::New(); - tracker->SetRandom(false); - tracker->SetInterpolateMasks(false); - tracker->SetNumberOfSamples(0); - tracker->SetAngularThreshold(-1); + +// tracker->SetInterpolateMasks(false); +// tracker->SetNumberOfSamples(0); +// tracker->SetAngularThreshold(-1); tracker->SetMaskImage(itk_mask_image); tracker->SetSeedImage(itk_seed_image); tracker->SetStoppingRegions(nullptr); - tracker->SetSeedsPerVoxel(1); - tracker->SetStepSize(0.5); - tracker->SetSamplingDistance(0.25); - tracker->SetUseStopVotes(true); - tracker->SetOnlyForwardSamples(true); - tracker->SetMinTractLength(20); - tracker->SetMaxNumTracts(-1); +// tracker->SetSeedsPerVoxel(1); +// tracker->SetStepSize(0.5); +// tracker->SetSamplingDistance(0.25); +// tracker->SetUseStopVotes(true); +// tracker->SetOnlyForwardSamples(true); +// tracker->SetMinTractLength(20); +// tracker->SetMaxNumTracts(-1); tracker->SetTrackingHandler(handler); - tracker->SetUseOutputProbabilityMap(false); +// tracker->SetUseOutputProbabilityMap(false); + tracker->SetParameters(params); } void tearDown() override { } void CheckFibResult(std::string ref_file, mitk::FiberBundle::Pointer test_fib) { mitk::FiberBundle::Pointer ref = LoadReferenceFib(ref_file); if (ref.IsNull()) { mitk::IOUtil::Save(test_fib, mitk::IOUtil::GetTempPath()+ref_file); CPPUNIT_FAIL("Reference file not found. Saving test file to " + mitk::IOUtil::GetTempPath() + ref_file); } else { bool is_equal = ref->Equals(test_fib); if (!is_equal) { mitk::IOUtil::Save(test_fib, mitk::IOUtil::GetTempPath()+ref_file); CPPUNIT_FAIL("Tractograms are not equal! Saving test file to " + mitk::IOUtil::GetTempPath() + ref_file); } } } void CheckImageResult(std::string ref_file, mitk::Image::Pointer test_img) { mitk::Image::Pointer ref = LoadReferenceImage(ref_file); if (ref.IsNull()) { mitk::IOUtil::Save(test_img, mitk::IOUtil::GetTempPath()+ref_file); CPPUNIT_FAIL("Reference file not found. Saving test file to " + mitk::IOUtil::GetTempPath() + ref_file); } else { MITK_ASSERT_EQUAL(test_img, ref, "Images should be equal"); } } void Test_Peak1() { mitk::TrackingHandlerPeaks* handler = new mitk::TrackingHandlerPeaks(); handler->SetPeakImage(itk_peak_image); - handler->SetPeakThreshold(peak_threshold); - + params->m_Cutoff = peak_threshold; SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Peak1.fib", outFib); delete handler; } void Test_Peak2() { mitk::TrackingHandlerPeaks* handler = new mitk::TrackingHandlerPeaks(); handler->SetPeakImage(itk_peak_image); - handler->SetPeakThreshold(peak_threshold); - handler->SetInterpolate(false); + params->m_Cutoff = peak_threshold; + params->m_InterpolateTractographyData = false; SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Peak2.fib", outFib); delete handler; } void Test_Tensor1() { mitk::TrackingHandlerTensor* handler = new mitk::TrackingHandlerTensor(); handler->SetTensorImage(itk_tensor_image); - handler->SetFaThreshold(gfa_threshold); - + params->m_Cutoff = gfa_threshold; SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Tensor1.fib", outFib); delete handler; } void Test_Tensor2() { mitk::TrackingHandlerTensor* handler = new mitk::TrackingHandlerTensor(); handler->SetTensorImage(itk_tensor_image); - handler->SetFaThreshold(gfa_threshold); - handler->SetInterpolate(false); + params->m_Cutoff = gfa_threshold; + params->m_InterpolateTractographyData = false; SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Tensor2.fib", outFib); delete handler; } void Test_Tensor3() { mitk::TrackingHandlerTensor* handler = new mitk::TrackingHandlerTensor(); handler->SetTensorImage(itk_tensor_image); - handler->SetFaThreshold(gfa_threshold); - handler->SetInterpolate(false); - handler->SetF(0); - handler->SetG(1); + params->m_Cutoff = gfa_threshold; + params->m_InterpolateTractographyData = false; + params->m_F = 0; + params->m_G = 1; SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Tensor3.fib", outFib); delete handler; } void Test_Odf1() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); - handler->SetGfaThreshold(gfa_threshold); - handler->SetOdfThreshold(0); - handler->SetSharpenOdfs(true); + + params->m_Cutoff = gfa_threshold; + params->m_OdfCutoff = 0; + params->m_SharpenOdfs = true; SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf1.fib", outFib); delete handler; } void Test_Odf2() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); - handler->SetGfaThreshold(gfa_threshold); - handler->SetOdfThreshold(0); - handler->SetSharpenOdfs(false); + + params->m_Cutoff = gfa_threshold; + params->m_OdfCutoff = 0; + params->m_SharpenOdfs = false; SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf2.fib", outFib); delete handler; } void Test_Odf3() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); - handler->SetGfaThreshold(gfa_threshold); - handler->SetOdfThreshold(0); - handler->SetSharpenOdfs(true); - handler->SetInterpolate(false); + + params->m_Cutoff = gfa_threshold; + params->m_OdfCutoff = 0; + params->m_SharpenOdfs = true; + params->m_InterpolateTractographyData = false; SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf3.fib", outFib); delete handler; } void Test_Odf4() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); - handler->SetGfaThreshold(gfa_threshold); - handler->SetOdfThreshold(0); - handler->SetSharpenOdfs(true); + + params->m_Cutoff = gfa_threshold; + params->m_OdfCutoff = 0; + params->m_SharpenOdfs = true; + params->m_SeedsPerVoxel = 3; SetupTracker(handler); - tracker->SetSeedsPerVoxel(3); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf4.fib", outFib); delete handler; } void Test_Odf5() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); - handler->SetGfaThreshold(gfa_threshold); - handler->SetOdfThreshold(0); - handler->SetSharpenOdfs(true); - handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); + + params->m_Cutoff = gfa_threshold; + params->m_OdfCutoff = 0; + params->m_SharpenOdfs = true; + params->m_SeedsPerVoxel = 3; + params->m_Mode = mitk::TrackingDataHandler::MODE::PROBABILISTIC; SetupTracker(handler); - tracker->SetSeedsPerVoxel(3); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf5.fib", outFib); delete handler; } void Test_Odf6() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); - handler->SetGfaThreshold(gfa_threshold); - handler->SetOdfThreshold(0); - handler->SetSharpenOdfs(true); - handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); + + params->m_Cutoff = gfa_threshold; + params->m_OdfCutoff = 0; + params->m_SharpenOdfs = true; + params->m_SeedsPerVoxel = 10; + params->m_Mode = mitk::TrackingDataHandler::MODE::PROBABILISTIC; + params->m_OutputProbMap = true; SetupTracker(handler); - tracker->SetSeedsPerVoxel(10); - tracker->SetUseOutputProbabilityMap(true); tracker->Update(); itk::StreamlineTrackingFilter::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::IOUtil::Save(img, mitk::IOUtil::GetTempPath()+"Test_Odf6.nrrd"); CheckImageResult("Test_Odf6.nrrd", img); delete handler; } }; MITK_TEST_SUITE_REGISTRATION(mitkStreamlineTractography) diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index 19ac812b48..39210be9cf 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,109 +1,111 @@ set(CPP_FILES mitkFiberTrackingModuleActivator.cpp ## IO datastructures IODataStructures/FiberBundle/mitkFiberBundle.cpp IODataStructures/FiberBundle/mitkTrackvis.cpp IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp IODataStructures/mitkTractographyForest.cpp IODataStructures/mitkFiberfoxParameters.cpp + IODataStructures/mitkStreamlineTractographyParameters.cpp # Interactions # Tractography Algorithms/GibbsTracking/mitkParticleGrid.cpp Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp Algorithms/GibbsTracking/mitkEnergyComputer.cpp Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp Algorithms/GibbsTracking/mitkFiberBuilder.cpp Algorithms/GibbsTracking/mitkSphereInterpolator.cpp Algorithms/itkStreamlineTrackingFilter.cpp Algorithms/TrackingHandlers/mitkTrackingDataHandler.cpp Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp ) set(H_FILES # DataStructures -> FiberBundle IODataStructures/FiberBundle/mitkFiberBundle.h IODataStructures/FiberBundle/mitkTrackvis.h IODataStructures/mitkFiberfoxParameters.h IODataStructures/mitkTractographyForest.h + IODataStructures/mitkStreamlineTractographyParameters.h # Algorithms Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkTractsToVectorImageFilter.h Algorithms/itkEvaluateDirectionImagesFilter.h Algorithms/itkEvaluateTractogramDirectionsFilter.h Algorithms/itkFiberCurvatureFilter.h Algorithms/itkFitFibersToImageFilter.h Algorithms/itkTractClusteringFilter.h Algorithms/itkTractDistanceFilter.h Algorithms/itkFiberExtractionFilter.h Algorithms/itkTdiToVolumeFractionFilter.h # Tractography Algorithms/TrackingHandlers/mitkTrackingDataHandler.h Algorithms/TrackingHandlers/mitkTrackingHandlerRandomForest.h Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.h Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.h Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.h Algorithms/itkGibbsTrackingFilter.h Algorithms/itkStochasticTractographyFilter.h Algorithms/GibbsTracking/mitkParticle.h Algorithms/GibbsTracking/mitkParticleGrid.h Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.h Algorithms/GibbsTracking/mitkSimpSamp.h Algorithms/GibbsTracking/mitkEnergyComputer.h Algorithms/GibbsTracking/mitkGibbsEnergyComputer.h Algorithms/GibbsTracking/mitkSphereInterpolator.h Algorithms/GibbsTracking/mitkFiberBuilder.h Algorithms/itkStreamlineTrackingFilter.h # Clustering Algorithms/ClusteringMetrics/mitkClusteringMetric.h Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMean.h Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanMax.h Algorithms/ClusteringMetrics/mitkClusteringMetricEuclideanStd.h Algorithms/ClusteringMetrics/mitkClusteringMetricAnatomic.h Algorithms/ClusteringMetrics/mitkClusteringMetricScalarMap.h Algorithms/ClusteringMetrics/mitkClusteringMetricInnerAngles.h Algorithms/ClusteringMetrics/mitkClusteringMetricLength.h # Fiberfox Fiberfox/itkFibersFromPlanarFiguresFilter.h Fiberfox/itkTractsToDWIImageFilter.h Fiberfox/itkKspaceImageFilter.h Fiberfox/itkDftImageFilter.h Fiberfox/itkFieldmapGeneratorFilter.h Fiberfox/itkRandomPhantomFilter.h Fiberfox/SignalModels/mitkDiffusionSignalModel.h Fiberfox/SignalModels/mitkTensorModel.h Fiberfox/SignalModels/mitkBallModel.h Fiberfox/SignalModels/mitkDotModel.h Fiberfox/SignalModels/mitkAstroStickModel.h Fiberfox/SignalModels/mitkStickModel.h Fiberfox/SignalModels/mitkRawShModel.h Fiberfox/SignalModels/mitkDiffusionNoiseModel.h Fiberfox/SignalModels/mitkRicianNoiseModel.h Fiberfox/SignalModels/mitkChiSquareNoiseModel.h Fiberfox/Sequences/mitkAcquisitionType.h Fiberfox/Sequences/mitkSingleShotEpi.h Fiberfox/Sequences/mitkConventionalSpinEcho.h Fiberfox/Sequences/mitkFastSpinEcho.h ) set(RESOURCE_FILES # Binary directory resources FiberTrackingLUTBaryCoords.bin FiberTrackingLUTIndices.bin ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/files.cmake b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/files.cmake index 42137f0e9c..aa364f9f67 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/files.cmake +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/files.cmake @@ -1,58 +1,47 @@ -set(SRC_CPP_FILES - QmitkMlbstTrainingDataWidget.cpp -) - set(INTERNAL_CPP_FILES QmitkGibbsTrackingView.cpp QmitkStreamlineTrackingView.cpp - QmitkMLBTView.cpp Perspectives/QmitkGibbsTractographyPerspective.cpp Perspectives/QmitkStreamlineTractographyPerspective.cpp mitkPluginActivator.cpp ) set(UI_FILES src/internal/QmitkGibbsTrackingViewControls.ui src/internal/QmitkStreamlineTrackingViewControls.ui - src/internal/QmitkMLBTViewControls.ui - - src/QmitkMlbstTrainingDataWidgetControls.ui ) set(MOC_H_FILES src/internal/mitkPluginActivator.h src/internal/QmitkGibbsTrackingView.h src/internal/QmitkStreamlineTrackingView.h - src/internal/QmitkMLBTView.h - - src/QmitkMlbstTrainingDataWidget.h src/internal/Perspectives/QmitkGibbsTractographyPerspective.h src/internal/Perspectives/QmitkStreamlineTractographyPerspective.h ) set(CACHED_RESOURCE_FILES plugin.xml resources/tract.png resources/tractogram.png resources/ml_tractography.png ) set(QRC_FILES ) set(CPP_FILES ) foreach(file ${SRC_CPP_FILES}) set(CPP_FILES ${CPP_FILES} src/${file}) endforeach(file ${SRC_CPP_FILES}) foreach(file ${INTERNAL_CPP_FILES}) set(CPP_FILES ${CPP_FILES} src/internal/${file}) endforeach(file ${INTERNAL_CPP_FILES}) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/QmitkMlbstTrainingDataWidget.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/QmitkMlbstTrainingDataWidget.cpp deleted file mode 100644 index 94f525fd33..0000000000 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/QmitkMlbstTrainingDataWidget.cpp +++ /dev/null @@ -1,65 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center, -Division of Medical and Biological Informatics. -All rights reserved. - -This software is distributed WITHOUT ANY WARRANTY; without -even the implied warranty of MERCHANTABILITY or FITNESS FOR -A PARTICULAR PURPOSE. - -See LICENSE.txt or http://www.mitk.org for details. - -===================================================================*/ - -//Qmitk headers -#include "QmitkMlbstTrainingDataWidget.h" - -#include -#include -#include -#include -#include -#include -#include - -const std::string QmitkMlbstTrainingDataWidget::VIEW_ID = "org.mitk.views.DotModelParametersWidget"; - -QmitkMlbstTrainingDataWidget::QmitkMlbstTrainingDataWidget( QWidget * parent, Qt::WindowFlags ) - : QWidget(parent) -{ - m_Controls = nullptr; - this->CreateQtPartControl(this); -} - -QmitkMlbstTrainingDataWidget::~QmitkMlbstTrainingDataWidget() -{ -} - -void QmitkMlbstTrainingDataWidget::CreateQtPartControl(QWidget *parent) -{ - if (!m_Controls) - { - // create GUI widgets - m_Controls = new Ui::QmitkMlbstTrainingDataWidgetControls; - m_Controls->setupUi(parent); - - mitk::NodePredicateIsDWI::Pointer isDiffusionImage = mitk::NodePredicateIsDWI::New(); - - mitk::TNodePredicateDataType::Pointer isFib = mitk::TNodePredicateDataType::New(); - - m_Controls->image->SetPredicate(isDiffusionImage); - m_Controls->fibers->SetPredicate(isFib); - - mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); - mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); - mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); - mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); - finalPredicate = mitk::NodePredicateAnd::New(finalPredicate, isBinaryPredicate); - - m_Controls->mask->SetPredicate(finalPredicate); - m_Controls->whiteMatter->SetPredicate(finalPredicate); - } -} diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/QmitkMlbstTrainingDataWidget.h b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/QmitkMlbstTrainingDataWidget.h deleted file mode 100644 index e728c4810a..0000000000 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/QmitkMlbstTrainingDataWidget.h +++ /dev/null @@ -1,69 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center, -Division of Medical and Biological Informatics. -All rights reserved. - -This software is distributed WITHOUT ANY WARRANTY; without -even the implied warranty of MERCHANTABILITY or FITNESS FOR -A PARTICULAR PURPOSE. - -See LICENSE.txt or http://www.mitk.org for details. - -===================================================================*/ -#ifndef _QmitkMlbstTrainingDataWidget_H_INCLUDED -#define _QmitkMlbstTrainingDataWidget_H_INCLUDED - -//QT headers -#include -#include -#include "ui_QmitkMlbstTrainingDataWidgetControls.h" -#include -#include -#include - -/** @brief - */ -class DIFFUSIONIMAGING_TRACTOGRAPHY_EXPORT QmitkMlbstTrainingDataWidget : public QWidget -{ - //this is needed for all Qt objects that should have a MOC object (everything that derives from QObject) - Q_OBJECT - -public: - - static const std::string VIEW_ID; - - QmitkMlbstTrainingDataWidget (QWidget* parent = nullptr, Qt::WindowFlags f = nullptr); - - virtual ~QmitkMlbstTrainingDataWidget(); - - virtual void CreateQtPartControl(QWidget *parent); - - mitk::DataNode::Pointer GetImage(){ return m_Controls->image->GetSelectedNode(); } - mitk::DataNode::Pointer GetFibers(){ return m_Controls->fibers->GetSelectedNode(); } - mitk::DataNode::Pointer GetMask(){ return m_Controls->mask->GetSelectedNode(); } - mitk::DataNode::Pointer GetWhiteMatter(){ return m_Controls->whiteMatter->GetSelectedNode(); } - void SetDataStorage(mitk::DataStorage::Pointer ds) - { - m_Controls->image->SetDataStorage(ds); - m_Controls->fibers->SetDataStorage(ds); - m_Controls->mask->SetDataStorage(ds); - m_Controls->whiteMatter->SetDataStorage(ds); - m_Controls->mask->SetZeroEntryText("--"); - m_Controls->whiteMatter->SetZeroEntryText("--"); - } - -public slots: - -protected: - // member variables - Ui::QmitkMlbstTrainingDataWidgetControls* m_Controls; - -private: - -}; - -#endif // _QmitkMlbstTrainingDataWidget_H_INCLUDED - diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/QmitkMlbstTrainingDataWidgetControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/QmitkMlbstTrainingDataWidgetControls.ui deleted file mode 100644 index b8b4cfe501..0000000000 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/QmitkMlbstTrainingDataWidgetControls.ui +++ /dev/null @@ -1,79 +0,0 @@ - - - QmitkMlbstTrainingDataWidgetControls - - - - 0 - 0 - 349 - 27 - - - - - 0 - 0 - - - - Form - - - - 0 - - - 0 - - - 0 - - - 0 - - - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - QmitkDataStorageComboBox - QComboBox -
QmitkDataStorageComboBox.h
- - - QmitkDataStorageComboBoxWithSelectNone - QComboBox -
QmitkDataStorageComboBoxWithSelectNone.h
- - - - - diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkMLBTView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkMLBTView.cpp deleted file mode 100644 index 3a017d6b22..0000000000 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkMLBTView.cpp +++ /dev/null @@ -1,246 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center, -Division of Medical and Biological Informatics. -All rights reserved. - -This software is distributed WITHOUT ANY WARRANTY; without -even the implied warranty of MERCHANTABILITY or FITNESS FOR -A PARTICULAR PURPOSE. - -See LICENSE.txt or http://www.mitk.org for details. - -===================================================================*/ - -#ifdef _MSC_VER -# pragma warning(push) -# pragma warning(disable: 4244) -#endif - -// Blueberry -#include -#include - -// Qmitk -#include "QmitkMLBTView.h" - -// Qt -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#define _USE_MATH_DEFINES -#include - -const std::string QmitkMLBTView::VIEW_ID = "org.mitk.views.mlbtview"; -using namespace berry; - -QmitkMLBTView::QmitkMLBTView() - : QmitkAbstractView() - , m_Controls( 0 ) -{ -} - -// Destructor -QmitkMLBTView::~QmitkMLBTView() -{ - delete m_ForestHandler; -} - -void QmitkMLBTView::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::QmitkMLBTViewControls; - m_Controls->setupUi( parent ); - - connect( m_Controls->m_StartTrainingButton, SIGNAL ( clicked() ), this, SLOT( StartTrainingThread() ) ); - connect( &m_TrainingWatcher, SIGNAL ( finished() ), this, SLOT( OnTrainingThreadStop() ) ); - connect( m_Controls->m_AddTwButton, SIGNAL ( clicked() ), this, SLOT( AddTrainingWidget() ) ); - connect( m_Controls->m_RemoveTwButton, SIGNAL ( clicked() ), this, SLOT( RemoveTrainingWidget() ) ); - - mitk::NodePredicateIsDWI::Pointer isDiffusionImage = mitk::NodePredicateIsDWI::New(); - - mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); - mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); - mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); - mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); - finalPredicate = mitk::NodePredicateAnd::New(finalPredicate, isBinaryPredicate); - AddTrainingWidget(); - } - -} - -void QmitkMLBTView::SetFocus() -{ -} - -void QmitkMLBTView::AddTrainingWidget() -{ - std::shared_ptr tw = std::make_shared(); - tw->SetDataStorage(this->GetDataStorage()); - m_Controls->m_TwFrame->layout()->addWidget(tw.get()); - m_TrainingWidgets.push_back(tw); -} - -void QmitkMLBTView::RemoveTrainingWidget() -{ - if(m_TrainingWidgets.size()>1) - { - m_TrainingWidgets.back().reset(); - m_TrainingWidgets.pop_back(); - } -} - -void QmitkMLBTView::StartTrainingThread() -{ - if (!this->IsTrainingInputValid()) - { - QMessageBox::warning(nullptr, "Training aborted", "Training could not be started. Not all necessary datasets were selected."); - return; - } - - if (m_Controls->m_FeatureTypeBox->currentIndex()==0) - m_ForestHandler = new mitk::TrackingHandlerRandomForest<6,28>(); - else - m_ForestHandler = new mitk::TrackingHandlerRandomForest<6,100>(); - - QFuture future = QtConcurrent::run( this, &QmitkMLBTView::StartTraining ); - m_TrainingWatcher.setFuture(future); - m_Controls->m_StartTrainingButton->setEnabled(false); - m_Controls->m_StartTrainingButton->setText("Training in progress ..."); - m_Controls->m_StartTrainingButton->setToolTip("Training in progress. This can take up to a couple of hours."); - m_Controls->m_StartTrainingButton->setCursor(Qt::WaitCursor); - QApplication::processEvents(); -} - -void QmitkMLBTView::OnTrainingThreadStop() -{ - m_Controls->m_StartTrainingButton->setEnabled(true); - m_Controls->m_StartTrainingButton->setCursor(Qt::ArrowCursor); - m_Controls->m_StartTrainingButton->setText("Start Training"); - m_Controls->m_StartTrainingButton->setToolTip("Start Training. This can take up to a couple of hours."); - - mitk::DataNode::Pointer node = mitk::DataNode::New(); - - switch (m_Controls->m_FeatureTypeBox->currentIndex()) - { - case 0: - if (dynamic_cast*>(m_ForestHandler)->IsForestValid()) - node->SetData(dynamic_cast*>(m_ForestHandler)->GetForest()); - break; - case 1: - if (dynamic_cast*>(m_ForestHandler)->IsForestValid()) - node->SetData(dynamic_cast*>(m_ForestHandler)->GetForest()); - break; - default: - mitkThrow() << "Unknown feature type!"; - } - - QString name("Forest"); - node->SetName(name.toStdString()); - GetDataStorage()->Add(node); - - delete m_ForestHandler; - - QApplication::processEvents(); -} - -void QmitkMLBTView::StartTraining() -{ - std::vector< mitk::Image::Pointer > m_SelectedDiffImages; - std::vector< mitk::FiberBundle::Pointer > m_SelectedFB; - std::vector< ItkUcharImgType::Pointer > m_MaskImages; - std::vector< ItkUcharImgType::Pointer > m_WhiteMatterImages; - - for (auto w : m_TrainingWidgets) - { - m_SelectedDiffImages.push_back(dynamic_cast(w->GetImage()->GetData())); - m_SelectedFB.push_back(dynamic_cast(w->GetFibers()->GetData())); - - if (w->GetMask().IsNotNull()) - { - mitk::Image::Pointer img = dynamic_cast(w->GetMask()->GetData()); - ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); - mitk::CastToItkImage(img, itkMask); - m_MaskImages.push_back(itkMask); - } - else - m_MaskImages.push_back(nullptr); - - if (w->GetWhiteMatter().IsNotNull()) - { - mitk::Image::Pointer img = dynamic_cast(w->GetWhiteMatter()->GetData()); - ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); - mitk::CastToItkImage(img, itkMask); - m_WhiteMatterImages.push_back(itkMask); - } - else - m_WhiteMatterImages.push_back(nullptr); - } - - switch (m_Controls->m_FeatureTypeBox->currentIndex()) - { - case 0: - dynamic_cast*>(m_ForestHandler)->SetDwis(m_SelectedDiffImages); - dynamic_cast*>(m_ForestHandler)->SetTractograms(m_SelectedFB); - dynamic_cast*>(m_ForestHandler)->SetMaskImages(m_MaskImages); - dynamic_cast*>(m_ForestHandler)->SetWhiteMatterImages(m_WhiteMatterImages); - dynamic_cast*>(m_ForestHandler)->SetNumTrees(m_Controls->m_NumTreesBox->value()); - dynamic_cast*>(m_ForestHandler)->SetMaxTreeDepth(m_Controls->m_MaxDepthBox->value()); - dynamic_cast*>(m_ForestHandler)->SetGrayMatterSamplesPerVoxel(m_Controls->m_GmSamplingBox->value()); - dynamic_cast*>(m_ForestHandler)->SetSampleFraction(m_Controls->m_SampleFractionBox->value()); - dynamic_cast*>(m_ForestHandler)->SetFiberSamplingStep(m_Controls->m_TrainingStepSizeBox->value()); - dynamic_cast*>(m_ForestHandler)->SetNumPreviousDirections(m_Controls->m_NumPrevDirs->value()); - dynamic_cast*>(m_ForestHandler)->StartTraining(); - break; - case 1: - dynamic_cast*>(m_ForestHandler)->SetDwis(m_SelectedDiffImages); - dynamic_cast*>(m_ForestHandler)->SetTractograms(m_SelectedFB); - dynamic_cast*>(m_ForestHandler)->SetMaskImages(m_MaskImages); - dynamic_cast*>(m_ForestHandler)->SetWhiteMatterImages(m_WhiteMatterImages); - dynamic_cast*>(m_ForestHandler)->SetNumTrees(m_Controls->m_NumTreesBox->value()); - dynamic_cast*>(m_ForestHandler)->SetMaxTreeDepth(m_Controls->m_MaxDepthBox->value()); - dynamic_cast*>(m_ForestHandler)->SetGrayMatterSamplesPerVoxel(m_Controls->m_GmSamplingBox->value()); - dynamic_cast*>(m_ForestHandler)->SetSampleFraction(m_Controls->m_SampleFractionBox->value()); - dynamic_cast*>(m_ForestHandler)->SetFiberSamplingStep(m_Controls->m_TrainingStepSizeBox->value()); - dynamic_cast*>(m_ForestHandler)->SetNumPreviousDirections(m_Controls->m_NumPrevDirs->value()); - dynamic_cast*>(m_ForestHandler)->StartTraining(); - break; - default: - mitkThrow() << "Unknown feature type!"; - } -} - -bool QmitkMLBTView::IsTrainingInputValid(void) const -{ - for (auto widget : m_TrainingWidgets) - { - if (widget->GetImage().IsNull() || widget->GetFibers().IsNull()) - { - return false; - } - } - - return true; -} - -#ifdef _MSC_VER -# pragma warning(pop) -#endif diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkMLBTView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkMLBTView.h deleted file mode 100644 index 801ffd8948..0000000000 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkMLBTView.h +++ /dev/null @@ -1,88 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center, -Division of Medical and Biological Informatics. -All rights reserved. - -This software is distributed WITHOUT ANY WARRANTY; without -even the implied warranty of MERCHANTABILITY or FITNESS FOR -A PARTICULAR PURPOSE. - -See LICENSE.txt or http://www.mitk.org for details. - -===================================================================*/ - -#include -#include "ui_QmitkMLBTViewControls.h" - -#ifndef Q_MOC_RUN -#include "mitkDataStorage.h" -#include "mitkDataStorageSelection.h" -#include -#include -#endif - -#include -#include -#include -#include -#include -#include - -/*! -\brief View to perform machine learning based fiber tractography. Includes training of the random forst classifier as well as the actual tractography. -*/ - -// Forward Qt class declarations - - -class QmitkMLBTView : public QmitkAbstractView -{ - - - // 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; - - typedef itk::Image ItkUcharImgType; - typedef itk::StreamlineTrackingFilter TrackerType; - - QmitkMLBTView(); - virtual ~QmitkMLBTView(); - virtual void CreateQtPartControl(QWidget *parent) override; - - /// - /// Sets the focus to an internal widget. - /// - virtual void SetFocus() override; - -protected slots: - - void StartTrainingThread(); - void OnTrainingThreadStop(); - void AddTrainingWidget(); - void RemoveTrainingWidget(); - -protected: - - void StartTraining(); - - Ui::QmitkMLBTViewControls* m_Controls; - mitk::TrackingDataHandler* m_ForestHandler; - QFutureWatcher m_TrainingWatcher; - std::vector< std::shared_ptr > m_TrainingWidgets; - -private: - - bool IsTrainingInputValid(void) const; - -}; - - - diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkMLBTViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkMLBTViewControls.ui deleted file mode 100644 index eea4153d87..0000000000 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkMLBTViewControls.ui +++ /dev/null @@ -1,444 +0,0 @@ - - - QmitkMLBTViewControls - - - - 0 - 0 - 476 - 548 - - - - Form - - - QCommandLinkButton:disabled { - border: none; -} - -QGroupBox { - background-color: transparent; -} - - - - - - QFrame::NoFrame - - - QFrame::Raised - - - - 0 - - - 0 - - - 0 - - - 0 - - - - - QFrame::NoFrame - - - QFrame::Raised - - - - 0 - - - 0 - - - 0 - - - 0 - - - - - Remove training data pair - - - ... - - - - :/org_mitk_icons/icons/tango/scalable/actions/list-remove.svg -:/org_mitk_icons/icons/tango/scalable/actions/list-remove.svg - - - - - - - Add additional training data pair - - - ... - - - - :/org_mitk_icons/icons/tango/scalable/actions/list-add.svg -:/org_mitk_icons/icons/tango/scalable/actions/list-add.svg - - - - - - - - - - <html><head/><body><p><span style=" font-weight:600; color:#ff0000;">It is recommended to use the command line application 'RfTraining' instead of this graphical user interface for the training process. Additional feature images are not supported here.</span></p><p><span style=" font-weight:600; color:#ff0000;">The GUI is intended for testing using small examples.</span></p></body></html> - - - Qt::RichText - - - Qt::AlignCenter - - - true - - - 5 - - - - - - - QFrame::NoFrame - - - QFrame::Raised - - - - 0 - - - 0 - - - 0 - - - 0 - - - 0 - - - 6 - - - - - QFrame::NoFrame - - - QFrame::Raised - - - - 0 - - - 0 - - - 0 - - - 0 - - - 0 - - - - - Input DWI: - - - Qt::AlignCenter - - - - - - - Reference Tractogram: - - - Qt::AlignCenter - - - - - - - Mask: - - - Qt::AlignCenter - - - - - - - WM: - - - Qt::AlignCenter - - - - - - - - - - - - - QFrame::NoFrame - - - QFrame::Raised - - - - 0 - - - 0 - - - 0 - - - 0 - - - - - Fraction of samples used to train each tree. - - - 3 - - - 1.000000000000000 - - - 0.100000000000000 - - - 0.700000000000000 - - - - - - - Num. Trees: - - - - - - - Sample Fraction: - - - - - - - Max. Depth: - - - - - - - Fiber sampling in mm. Determines the number of white-matter samples (-1 = auto). - - - 3 - - - -1.000000000000000 - - - 999.000000000000000 - - - 0.100000000000000 - - - -1.000000000000000 - - - - - - - Number of tress in the final random forest. - - - 1 - - - 999999999 - - - 30 - - - - - - - Non-WM Sampling Points: - - - - - - - Fiber Sampling: - - - - - - - Maximum tree depth. - - - 1 - - - 999999999 - - - 25 - - - - - - - Number of sampling points outside of the white-matter (-1 = automatic estimation). - - - -1 - - - 999999999 - - - -1 - - - - - - - Num. previous directions: - - - - - - - Number of previous fiber directions used as features. - - - - - - 0 - - - 50 - - - 1 - - - - - - - dMRI Features: - - - - - - - - Spherical Harmonics Coefficients - - - - - Raw Data - - - - - - - - - - - - - - Start Training. This can take up to a couple of hours. - - - Start Training - - - - - - - Qt::Vertical - - - - 20 - 40 - - - - - - - - - - - \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp index bbac01b9d5..e00e0b7373 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp @@ -1,1008 +1,1060 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkStreamlineTrackingView.h" #include "QmitkStdMultiWidget.h" // Qt #include +#include // MITK #include #include #include #include #include #include #include #include #include #include #include #include #include #include // VTK #include #include #include #include #include #include #include #include #include #include const std::string QmitkStreamlineTrackingView::VIEW_ID = "org.mitk.views.streamlinetracking"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace berry; QmitkStreamlineTrackingWorker::QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view) : m_View(view) { } void QmitkStreamlineTrackingWorker::run() { m_View->m_Tracker->Update(); m_View->m_TrackingThread.quit(); } QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : m_TrackingWorker(this) , m_Controls(nullptr) , m_FirstTensorProbRun(true) , m_FirstInteractiveRun(true) , m_TrackingHandler(nullptr) , m_ThreadIsRunning(false) , m_DeleteTrackingHandler(false) , m_Visible(false) , m_LastPrior(nullptr) , m_TrackingPriorHandler(nullptr) { m_TrackingWorker.moveToThread(&m_TrackingThread); connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); m_TrackingTimer = new QTimer(this); } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); m_TrackingThread.wait(); } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStreamlineTrackingViewControls; m_Controls->setupUi( parent ); m_Controls->m_FaImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_SeedImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_TargetImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_PriorImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_StopImageSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_ForestSelectionWidget->SetDataStorage(this->GetDataStorage()); m_Controls->m_ExclusionImageSelectionWidget->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isPeakImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTractographyForest = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New( isBinaryPredicate ); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New( isImagePredicate, isNotBinaryPredicate ); mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); m_Controls->m_ForestSelectionWidget->SetNodePredicate(isTractographyForest); m_Controls->m_FaImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_FaImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_FaImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_SeedImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_SeedImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_SeedImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_MaskImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_MaskImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_MaskImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_StopImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_StopImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_StopImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_TargetImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_TargetImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_TargetImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_PriorImageSelectionWidget->SetNodePredicate( isPeakImagePredicate ); m_Controls->m_PriorImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_PriorImageSelectionWidget->SetSelectionIsOptional(true); m_Controls->m_ExclusionImageSelectionWidget->SetNodePredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_ExclusionImageSelectionWidget->SetEmptyInfo("--"); m_Controls->m_ExclusionImageSelectionWidget->SetSelectionIsOptional(true); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); + connect( m_Controls->m_SaveParametersButton, SIGNAL(clicked()), this, SLOT(SaveParameters()) ); connect( m_Controls->commandLinkButton_2, SIGNAL(clicked()), this, SLOT(StopTractography()) ); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_InteractiveBox, SIGNAL(stateChanged(int)), this, SLOT(ToggleInteractive()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_FaImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OutputStyleSwitched()) ); connect( m_Controls->m_SeedImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_TargetImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_PriorImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ExclusionImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FaImageSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ForestSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(ForestSwitched()) ); connect( m_Controls->m_ForestSelectionWidget, SIGNAL(CurrentSelectionChanged(QList)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedsPerVoxelBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumFibersBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ScalarThresholdBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_OdfCutoffBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StepSizeBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SamplingDistanceBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_AngularThresholdBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MinTractLengthBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_fBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_gBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumSamplesBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedRadiusBox, SIGNAL(editingFinished()), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_NumSeedsBox, SIGNAL(editingFinished()), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SharpenOdfsBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_InterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskInterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipXBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipYBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipZBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FrontalSamplesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopVotesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_LoopCheckBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_TrialsPerSeedBox, SIGNAL(editingFinished()), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_EpConstraintsBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); m_Controls->m_SeedsPerVoxelBox->editingFinished(); m_Controls->m_NumFibersBox->editingFinished(); m_Controls->m_ScalarThresholdBox->editingFinished(); m_Controls->m_OdfCutoffBox->editingFinished(); m_Controls->m_StepSizeBox->editingFinished(); m_Controls->m_SamplingDistanceBox->editingFinished(); m_Controls->m_AngularThresholdBox->editingFinished(); m_Controls->m_MinTractLengthBox->editingFinished(); m_Controls->m_fBox->editingFinished(); m_Controls->m_gBox->editingFinished(); m_Controls->m_NumSamplesBox->editingFinished(); m_Controls->m_SeedRadiusBox->editingFinished(); m_Controls->m_NumSeedsBox->editingFinished(); m_Controls->m_LoopCheckBox->editingFinished(); m_Controls->m_TrialsPerSeedBox->editingFinished(); StartStopTrackingGui(false); } + m_ParameterFile = QDir::currentPath()+"/param.stp"; UpdateGui(); } +std::shared_ptr QmitkStreamlineTrackingView::GetParametersFromGui() +{ + std::shared_ptr params = std::make_shared(); + params->m_InteractiveRadius = m_Controls->m_SeedRadiusBox->value(); + params->m_MaxNumFibers = m_Controls->m_NumFibersBox->value(); + params->m_Cutoff = static_cast(m_Controls->m_ScalarThresholdBox->value()); + params->m_F = static_cast(m_Controls->m_fBox->value()); + params->m_G = static_cast(m_Controls->m_gBox->value()); + + params->m_OdfCutoff = static_cast(m_Controls->m_OdfCutoffBox->value()); + params->m_SharpenOdfs = m_Controls->m_SharpenOdfsBox->isChecked(); + + params->m_Weight = static_cast(m_Controls->m_PriorWeightBox->value()); + params->m_RestrictToPrior = m_Controls->m_PriorAsMaskBox->isChecked(); + params->m_NewDirectionsFromPrior = m_Controls->m_NewDirectionsFromPriorBox->isChecked(); + + + params->m_FlipX = m_Controls->m_FlipXBox->isChecked(); + params->m_FlipY = m_Controls->m_FlipYBox->isChecked(); + params->m_FlipZ = m_Controls->m_FlipZBox->isChecked(); + params->m_InterpolateTractographyData = m_Controls->m_InterpolationBox->isChecked(); + + + params->m_InterpolateRoiImages = m_Controls->m_MaskInterpolationBox->isChecked(); + params->m_SeedsPerVoxel = m_Controls->m_SeedsPerVoxelBox->value(); + params->SetStepSizeVox(m_Controls->m_StepSizeBox->value()); + params->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value()); + params->m_StopVotes = m_Controls->m_StopVotesBox->isChecked(); + params->m_OnlyForwardSamples = m_Controls->m_FrontalSamplesBox->isChecked(); + params->m_TrialsPerSeed = m_Controls->m_TrialsPerSeedBox->value(); + + params->m_NumSamples = m_Controls->m_NumSamplesBox->value(); + params->SetLoopCheckDeg(m_Controls->m_LoopCheckBox->value()); + params->SetAngularThresholdDeg(m_Controls->m_AngularThresholdBox->value()); + params->m_MinTractLength = m_Controls->m_MinTractLengthBox->value(); + params->m_OutputProbMap = m_Controls->m_OutputProbMap->isChecked(); + params->m_FixRandomSeed = m_Controls->m_FixSeedBox->isChecked(); + + switch (m_Controls->m_ModeBox->currentIndex()) + { + case 0: + params->m_Mode = mitk::TrackingDataHandler::MODE::DETERMINISTIC; + break; + case 1: + params->m_Mode = mitk::TrackingDataHandler::MODE::PROBABILISTIC; + break; + default: + params->m_Mode = mitk::TrackingDataHandler::MODE::DETERMINISTIC; + } + + switch (m_Controls->m_EpConstraintsBox->currentIndex()) + { + case 0: + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::NONE; + m_Tracker->SetTargetRegions(nullptr); + break; + case 1: + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET; + break; + case 2: + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF; + break; + case 3: + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET; + break; + case 4: + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET; + break; + case 5: + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET; + break; + case 6: + params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET; + break; + } + + return params; +} + +void QmitkStreamlineTrackingView::SaveParameters() +{ + QString filename = QFileDialog::getSaveFileName( + 0, + tr("Save Parameters"), + m_ParameterFile, + tr("Streamline Tractography Parameters (*.stp)") ); + + m_ParameterFile = filename; + + auto params = GetParametersFromGui(); + params->SaveParameters(m_ParameterFile.toStdString()); +} + void QmitkStreamlineTrackingView::StopTractography() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); } void QmitkStreamlineTrackingView::TimerUpdate() { if (m_Tracker.IsNull()) return; QString status_text(m_Tracker->GetStatusText().c_str()); m_Controls->m_StatusTextBox->setText(status_text); } void QmitkStreamlineTrackingView::BeforeThread() { m_TrackingTimer->start(1000); } void QmitkStreamlineTrackingView::AfterThread() { + auto params = m_Tracker->GetParameters(); m_TrackingTimer->stop(); - if (!m_Tracker->GetUseOutputProbabilityMap()) + if (!params->m_OutputProbMap) { vtkSmartPointer fiberBundle = m_Tracker->GetFiberPolyData(); if (!m_Controls->m_InteractiveBox->isChecked() && fiberBundle->GetNumberOfLines() == 0) { QMessageBox warnBox; warnBox.setWindowTitle("Warning"); warnBox.setText("No fiberbundle was generated!"); warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some images feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter.\n- A small step sizes also means many steps to go wrong. Especially in the case of probabilistic tractography. Try to adjust the angular threshold."); warnBox.setIcon(QMessageBox::Warning); warnBox.exec(); if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); return; } mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); fib->SetReferenceGeometry(dynamic_cast(m_ParentNode->GetData())->GetGeometry()); if (m_Controls->m_ResampleFibersBox->isChecked() && fiberBundle->GetNumberOfLines()>0) fib->Compress(m_Controls->m_FiberErrorBox->value()); fib->ColorFibersByOrientation(); m_Tracker->SetDicomProperties(fib); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(fib); m_InteractiveNode->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(fib); QString name("FiberBundle_"); name += m_ParentNode->GetName().c_str(); name += "_Streamline"; node->SetName(name.toStdString()); node->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); GetDataStorage()->Add(node, m_ParentNode); } } else { TrackerType::ItkDoubleImgType::Pointer outImg = m_Tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(img); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); m_InteractiveNode->SetProperty("LookupTable", lut_prop); m_InteractiveNode->SetProperty("opacity", mitk::FloatProperty::New(0.5)); m_InteractiveNode->SetFloatProperty("Fiber2DSliceThickness", m_Tracker->GetMinVoxelSize()/2); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); QString name("ProbabilityMap_"); name += m_ParentNode->GetName().c_str(); node->SetName(name.toStdString()); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); node->SetProperty("LookupTable", lut_prop); node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); GetDataStorage()->Add(node, m_ParentNode); } } if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); } void QmitkStreamlineTrackingView::InteractiveSeedChanged(bool posChanged) { if(!CheckAndStoreLastParams(sender()) && !posChanged) return; if (m_ThreadIsRunning || !m_Visible) return; if (!posChanged && (!m_Controls->m_InteractiveBox->isChecked() || !m_Controls->m_ParamUpdateBox->isChecked()) ) return; std::srand(std::time(0)); m_SeedPoints.clear(); itk::Point world_pos = this->GetRenderWindowPart()->GetSelectedPosition(); m_SeedPoints.push_back(world_pos); float radius = m_Controls->m_SeedRadiusBox->value(); int num = m_Controls->m_NumSeedsBox->value(); mitk::PointSet::Pointer pointset = mitk::PointSet::New(); pointset->InsertPoint(0, world_pos); m_InteractivePointSetNode->SetProperty("pointsize", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetProperty("point 2D size", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetData(pointset); for (int i=1; i p; p[0] = rand()%1000-500; p[1] = rand()%1000-500; p[2] = rand()%1000-500; p.Normalize(); p *= radius; m_SeedPoints.push_back(world_pos+p); } m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,0,0)); DoFiberTracking(); } bool QmitkStreamlineTrackingView::CheckAndStoreLastParams(QObject* obj) { if (obj!=nullptr) { std::string new_val = ""; if(qobject_cast(obj)!=nullptr) new_val = boost::lexical_cast(qobject_cast(obj)->value()); else if (qobject_cast(obj)!=nullptr) new_val = boost::lexical_cast(qobject_cast(obj)->value()); else return true; if (m_LastTractoParams.find(obj->objectName())==m_LastTractoParams.end()) { m_LastTractoParams[obj->objectName()] = new_val; return false; } else if (m_LastTractoParams.at(obj->objectName()) != new_val) { m_LastTractoParams[obj->objectName()] = new_val; return true; } else if (m_LastTractoParams.at(obj->objectName()) == new_val) return false; } return true; } void QmitkStreamlineTrackingView::OnParameterChanged() { UpdateGui(); if(!CheckAndStoreLastParams(sender())) return; if (m_Controls->m_InteractiveBox->isChecked() && m_Controls->m_ParamUpdateBox->isChecked()) DoFiberTracking(); } void QmitkStreamlineTrackingView::ToggleInteractive() { UpdateGui(); m_Controls->m_SeedsPerVoxelBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedsPerVoxelLabel->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedImageSelectionWidget->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_6->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); if ( m_Controls->m_InteractiveBox->isChecked() ) { if (m_FirstInteractiveRun) { QMessageBox::information(nullptr, "Information", "Place and move a spherical seed region anywhere in the image by left-clicking and dragging. If the seed region is colored red, tracking is in progress. If the seed region is colored white, tracking is finished.\nPlacing the seed region for the first time in a newly selected dataset might cause a short delay, since the tracker needs to be initialized."); m_FirstInteractiveRun = false; } QApplication::setOverrideCursor(Qt::PointingHandCursor); QApplication::processEvents(); m_InteractivePointSetNode = mitk::DataNode::New(); m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); m_InteractivePointSetNode->SetName("InteractiveSeedRegion"); mitk::PointSetShapeProperty::Pointer shape_prop = mitk::PointSetShapeProperty::New(); shape_prop->SetValue(mitk::PointSetShapeProperty::PointSetShape::CIRCLE); m_InteractivePointSetNode->SetProperty("Pointset.2D.shape", shape_prop); GetDataStorage()->Add(m_InteractivePointSetNode); m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } else { QApplication::restoreOverrideCursor(); QApplication::processEvents(); m_InteractiveNode = nullptr; m_InteractivePointSetNode = nullptr; m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); disconnect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } } void QmitkStreamlineTrackingView::Activated() { } void QmitkStreamlineTrackingView::Deactivated() { } void QmitkStreamlineTrackingView::Visible() { m_Visible = true; } void QmitkStreamlineTrackingView::Hidden() { m_Visible = false; m_Controls->m_InteractiveBox->setChecked(false); ToggleInteractive(); } void QmitkStreamlineTrackingView::OnSliceChanged() { InteractiveSeedChanged(true); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::DeleteTrackingHandler() { if (!m_ThreadIsRunning && m_TrackingHandler != nullptr) { if (m_TrackingPriorHandler != nullptr) { delete m_TrackingPriorHandler; m_TrackingPriorHandler = nullptr; } delete m_TrackingHandler; m_TrackingHandler = nullptr; m_DeleteTrackingHandler = false; m_LastPrior = nullptr; } else if (m_ThreadIsRunning) { m_DeleteTrackingHandler = true; } } void QmitkStreamlineTrackingView::ForestSwitched() { DeleteTrackingHandler(); } void QmitkStreamlineTrackingView::OutputStyleSwitched() { if (m_InteractiveNode.IsNotNull()) GetDataStorage()->Remove(m_InteractiveNode); m_InteractiveNode = nullptr; } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer , const QList& nodes ) { std::vector< mitk::DataNode::Pointer > last_nodes = m_InputImageNodes; m_InputImageNodes.clear(); m_AdditionalInputImages.clear(); bool retrack = false; for( auto node : nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { if( dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(node->GetData()))) { m_InputImageNodes.push_back(node); retrack = true; } else { mitk::Image* img = dynamic_cast(node->GetData()); if (img!=nullptr && img->GetDimension()==3) m_AdditionalInputImages.push_back(dynamic_cast(node->GetData())); } } } // sometimes the OnSelectionChanged event is sent twice and actually no selection has changed for the first event. We need to catch that. if (last_nodes.size() == m_InputImageNodes.size()) { bool same_nodes = true; for (unsigned int i=0; im_TensorImageLabel->setText("select in data-manager"); m_Controls->m_fBox->setEnabled(false); m_Controls->m_fLabel->setEnabled(false); m_Controls->m_gBox->setEnabled(false); m_Controls->m_gLabel->setEnabled(false); m_Controls->m_FaImageSelectionWidget->setEnabled(true); m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_OdfCutoffBox->setEnabled(false); m_Controls->m_OdfCutoffLabel->setEnabled(false); m_Controls->m_SharpenOdfsBox->setEnabled(false); m_Controls->m_ForestSelectionWidget->setVisible(false); m_Controls->m_ForestLabel->setVisible(false); m_Controls->commandLinkButton->setEnabled(false); m_Controls->m_TrialsPerSeedBox->setEnabled(false); m_Controls->m_TrialsPerSeedLabel->setEnabled(false); m_Controls->m_TargetImageSelectionWidget->setEnabled(false); m_Controls->m_TargetImageLabel->setEnabled(false); if (m_Controls->m_InteractiveBox->isChecked()) { m_Controls->m_InteractiveSeedingFrame->setVisible(true); m_Controls->m_StaticSeedingFrame->setVisible(false); m_Controls->commandLinkButton_2->setVisible(false); m_Controls->commandLinkButton->setVisible(false); } else { m_Controls->m_InteractiveSeedingFrame->setVisible(false); m_Controls->m_StaticSeedingFrame->setVisible(true); m_Controls->commandLinkButton_2->setVisible(m_ThreadIsRunning); m_Controls->commandLinkButton->setVisible(!m_ThreadIsRunning); } if (m_Controls->m_EpConstraintsBox->currentIndex()>0) { m_Controls->m_TargetImageSelectionWidget->setEnabled(true); m_Controls->m_TargetImageLabel->setEnabled(true); } // trials per seed are only important for probabilistic tractography if (m_Controls->m_ModeBox->currentIndex()==1) { m_Controls->m_TrialsPerSeedBox->setEnabled(true); m_Controls->m_TrialsPerSeedLabel->setEnabled(true); } if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText( ( std::to_string(m_InputImageNodes.size()) + " images selected").c_str() ); else m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked() && !m_ThreadIsRunning); m_Controls->m_ScalarThresholdBox->setEnabled(true); m_Controls->m_FaThresholdLabel->setEnabled(true); if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setEnabled(true); m_Controls->m_fLabel->setEnabled(true); m_Controls->m_gBox->setEnabled(true); m_Controls->m_gLabel->setEnabled(true); } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) || dynamic_cast(m_InputImageNodes.at(0)->GetData())) { m_Controls->m_OdfCutoffBox->setEnabled(true); m_Controls->m_OdfCutoffLabel->setEnabled(true); m_Controls->m_SharpenOdfsBox->setEnabled(true); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { m_Controls->m_ForestSelectionWidget->setVisible(true); m_Controls->m_ForestLabel->setVisible(true); m_Controls->m_ScalarThresholdBox->setEnabled(false); m_Controls->m_FaThresholdLabel->setEnabled(false); } } } void QmitkStreamlineTrackingView::StartStopTrackingGui(bool start) { m_ThreadIsRunning = start; if (!m_Controls->m_InteractiveBox->isChecked()) { m_Controls->commandLinkButton_2->setVisible(start); m_Controls->commandLinkButton->setVisible(!start); m_Controls->m_InteractiveBox->setEnabled(!start); m_Controls->m_StatusTextBox->setVisible(start); } } void QmitkStreamlineTrackingView::DoFiberTracking() { + auto params = GetParametersFromGui(); + if (m_InputImageNodes.empty()) { QMessageBox::information(nullptr, "Information", "Please select an input image in the datamaneger (tensor, ODF, peak or dMRI image)!"); return; } if (m_ThreadIsRunning || !m_Visible) return; if (m_Controls->m_InteractiveBox->isChecked() && m_SeedPoints.empty()) return; StartStopTrackingGui(true); m_Tracker = TrackerType::New(); if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_Controls->m_ModeBox->currentIndex()==1) { if (m_InputImageNodes.size()>1) { QMessageBox::information(nullptr, "Information", "Probabilistic tensor tractography is only implemented for single-tensor mode!"); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerOdf(); typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( mitk::convert::GetItkTensorFromTensorImage(dynamic_cast(m_InputImageNodes.at(0)->GetData())) ); filter->Update(); dynamic_cast(m_TrackingHandler)->SetOdfImage(filter->GetOutput()); if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } - dynamic_cast(m_TrackingHandler)->SetGfaThreshold(static_cast(m_Controls->m_ScalarThresholdBox->value())); - dynamic_cast(m_TrackingHandler)->SetOdfThreshold(0); - dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(true); dynamic_cast(m_TrackingHandler)->SetIsOdfFromTensor(true); } else { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerTensor(); for (unsigned int i=0; i(m_TrackingHandler)->AddTensorImage(mitk::convert::GetItkTensorFromTensorImage(dynamic_cast(m_InputImageNodes.at(i)->GetData())).GetPointer()); if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetFaImage(itkImg); } } - - dynamic_cast(m_TrackingHandler)->SetFaThreshold(static_cast(m_Controls->m_ScalarThresholdBox->value())); - dynamic_cast(m_TrackingHandler)->SetF(static_cast(m_Controls->m_fBox->value())); - dynamic_cast(m_TrackingHandler)->SetG(static_cast(m_Controls->m_gBox->value())); } } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) || dynamic_cast(m_InputImageNodes.at(0)->GetData())) { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerOdf(); if (dynamic_cast(m_InputImageNodes.at(0)->GetData())) dynamic_cast(m_TrackingHandler)->SetOdfImage(mitk::convert::GetItkOdfFromShImage(dynamic_cast(m_InputImageNodes.at(0)->GetData()))); else dynamic_cast(m_TrackingHandler)->SetOdfImage(mitk::convert::GetItkOdfFromOdfImage(dynamic_cast(m_InputImageNodes.at(0)->GetData()))); if (m_Controls->m_FaImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageSelectionWidget->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } - - dynamic_cast(m_TrackingHandler)->SetGfaThreshold(static_cast(m_Controls->m_ScalarThresholdBox->value())); - dynamic_cast(m_TrackingHandler)->SetOdfThreshold(static_cast(m_Controls->m_OdfCutoffBox->value())); - dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(m_Controls->m_SharpenOdfsBox->isChecked()); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { if ( m_Controls->m_ForestSelectionWidget->GetSelectedNode().IsNull() ) { QMessageBox::information(nullptr, "Information", "Not random forest for machine learning based tractography (raw dMRI tractography) selected. Did you accidentally select the raw diffusion-weighted image in the datamanager?"); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { mitk::TractographyForest::Pointer forest = dynamic_cast(m_Controls->m_ForestSelectionWidget->GetSelectedNode()->GetData()); mitk::Image::Pointer dwi = dynamic_cast(m_InputImageNodes.at(0)->GetData()); std::vector< std::vector< ItkFloatImageType::Pointer > > additionalFeatureImages; additionalFeatureImages.push_back(std::vector< ItkFloatImageType::Pointer >()); for (auto img : m_AdditionalInputImages) { ItkFloatImageType::Pointer itkimg = ItkFloatImageType::New(); mitk::CastToItkImage(img, itkimg); additionalFeatureImages.at(0).push_back(itkimg); } bool forest_valid = false; if (forest->GetNumFeatures()>=100) { - unsigned int num_previous_directions = static_cast((forest->GetNumFeatures() - (100 + additionalFeatureImages.at(0).size()))/3); + params->m_NumPreviousDirections = static_cast((forest->GetNumFeatures() - (100 + additionalFeatureImages.at(0).size()))/3); m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 100>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); - dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } else { - unsigned int num_previous_directions = static_cast((forest->GetNumFeatures() - (28 + additionalFeatureImages.at(0).size()))/3); + params->m_NumPreviousDirections = static_cast((forest->GetNumFeatures() - (28 + additionalFeatureImages.at(0).size()))/3); m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 28>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); - dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } if (!forest_valid) { QMessageBox::information(nullptr, "Information", "Random forest is invalid. The forest signatue does not match the parameters of TrackingHandlerRandomForest."); StartStopTrackingGui(false); return; } } } else { + if (m_Controls->m_ModeBox->currentIndex()==1) + { + QMessageBox::information(nullptr, "Information", "Probabilstic tractography is not implemented for peak images."); + StartStopTrackingGui(false); + return; + } if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingHandler)->SetPeakImage(mitk::convert::GetItkPeakFromPeakImage(dynamic_cast(m_InputImageNodes.at(0)->GetData()))); } - - dynamic_cast(m_TrackingHandler)->SetPeakThreshold(static_cast(m_Controls->m_ScalarThresholdBox->value())); - } - - m_TrackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); - m_TrackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); - m_TrackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); - m_TrackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); - switch (m_Controls->m_ModeBox->currentIndex()) - { - case 0: - m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); - break; - case 1: - m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); - break; - default: - m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); } if (m_Controls->m_InteractiveBox->isChecked()) { m_Tracker->SetSeedPoints(m_SeedPoints); } else if (m_Controls->m_SeedImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetSeedImage(mask); } if (m_Controls->m_MaskImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetMaskImage(mask); } if (m_Controls->m_StopImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetStoppingRegions(mask); } if (m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_TargetImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetTargetRegions(mask); } if (m_Controls->m_PriorImageSelectionWidget->GetSelectedNode().IsNotNull()) { + auto prior_params = GetParametersFromGui(); if (m_LastPrior!=m_Controls->m_PriorImageSelectionWidget->GetSelectedNode() || m_TrackingPriorHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(dynamic_cast(m_Controls->m_PriorImageSelectionWidget->GetSelectedNode()->GetData())); caster->SetCopyMemFlag(true); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); m_TrackingPriorHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingPriorHandler)->SetPeakImage(itkImg); - dynamic_cast(m_TrackingPriorHandler)->SetPeakThreshold(0.0); m_LastPrior = m_Controls->m_PriorImageSelectionWidget->GetSelectedNode(); } - - m_TrackingPriorHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); switch (m_Controls->m_ModeBox->currentIndex()) { case 0: m_TrackingPriorHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); break; case 1: m_TrackingPriorHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); break; default: m_TrackingPriorHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); } m_TrackingPriorHandler->SetFlipX(m_Controls->m_PriorFlipXBox->isChecked()); m_TrackingPriorHandler->SetFlipY(m_Controls->m_PriorFlipYBox->isChecked()); m_TrackingPriorHandler->SetFlipZ(m_Controls->m_PriorFlipZBox->isChecked()); + m_TrackingPriorHandler->SetParameters(prior_params); m_Tracker->SetTrackingPriorHandler(m_TrackingPriorHandler); - m_Tracker->SetTrackingPriorWeight(m_Controls->m_PriorWeightBox->value()); - m_Tracker->SetTrackingPriorAsMask(m_Controls->m_PriorAsMaskBox->isChecked()); - m_Tracker->SetIntroduceDirectionsFromPrior(m_Controls->m_NewDirectionsFromPriorBox->isChecked()); } else if (m_Controls->m_PriorImageSelectionWidget->GetSelectedNode().IsNull()) m_Tracker->SetTrackingPriorHandler(nullptr); if (m_Controls->m_ExclusionImageSelectionWidget->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_ExclusionImageSelectionWidget->GetSelectedNode()->GetData()), mask); m_Tracker->SetExclusionRegions(mask); } - // Endpoint constraints - switch (m_Controls->m_EpConstraintsBox->currentIndex()) - { - case 0: - m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NONE); + if (params->m_EpConstraints == itk::StreamlineTrackingFilter::EndpointConstraints::NONE) m_Tracker->SetTargetRegions(nullptr); break; case 1: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET); break; case 2: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF); break; case 3: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET); break; case 4: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET); break; case 5: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET); break; case 6: m_Tracker->SetEndpointConstraint(itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET); break; } - if (m_Tracker->GetEndpointConstraint()!=itk::StreamlineTrackingFilter::EndpointConstraints::NONE && m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNull()) + if (params->m_EpConstraints!=itk::StreamlineTrackingFilter::EndpointConstraints::NONE && m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Error", "Endpoint constraints are used but no target image is set!"); StartStopTrackingGui(false); return; } - else if (m_Tracker->GetEndpointConstraint()==itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET + else if (params->m_EpConstraints==itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET && (m_Controls->m_SeedImageSelectionWidget->GetSelectedNode().IsNull()|| m_Controls->m_TargetImageSelectionWidget->GetSelectedNode().IsNull()) ) { QMessageBox::information(nullptr, "Error", "Endpoint constraint EPS_IN_SEED_AND_TARGET is used but no target or no seed image is set!"); StartStopTrackingGui(false); return; } - m_Tracker->SetInterpolateMasks(m_Controls->m_MaskInterpolationBox->isChecked()); - m_Tracker->SetVerbose(!m_Controls->m_InteractiveBox->isChecked()); - m_Tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); - m_Tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); - m_Tracker->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value()); - m_Tracker->SetUseStopVotes(m_Controls->m_StopVotesBox->isChecked()); - m_Tracker->SetOnlyForwardSamples(m_Controls->m_FrontalSamplesBox->isChecked()); - m_Tracker->SetTrialsPerSeed(m_Controls->m_TrialsPerSeedBox->value()); - m_Tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); - m_Tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); m_Tracker->SetTrackingHandler(m_TrackingHandler); - m_Tracker->SetLoopCheck(m_Controls->m_LoopCheckBox->value()); - m_Tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); - m_Tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); - m_Tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); - m_Tracker->SetRandom(!m_Controls->m_FixSeedBox->isChecked()); + m_Tracker->SetVerbose(!m_Controls->m_InteractiveBox->isChecked()); m_ParentNode = m_InputImageNodes.at(0); m_TrackingThread.start(QThread::LowestPriority); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h index f2d3bcb433..3e392063b0 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h @@ -1,154 +1,158 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef QmitkStreamlineTrackingView_h #define QmitkStreamlineTrackingView_h #include #include "ui_QmitkStreamlineTrackingViewControls.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include class QmitkStreamlineTrackingView; class QmitkStreamlineTrackingWorker : public QObject { Q_OBJECT public: QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view); public slots: void run(); private: QmitkStreamlineTrackingView* m_View; }; /*! \brief View for tensor based deterministic streamline fiber tracking. */ class QmitkStreamlineTrackingView : 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; typedef itk::Image< unsigned int, 3 > ItkUintImgType; typedef itk::Image< unsigned char, 3 > ItkUCharImageType; typedef itk::Image< float, 3 > ItkFloatImageType; typedef itk::StreamlineTrackingFilter TrackerType; QmitkStreamlineTrackingView(); virtual ~QmitkStreamlineTrackingView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; TrackerType::Pointer m_Tracker; QmitkStreamlineTrackingWorker m_TrackingWorker; QThread m_TrackingThread; virtual void Activated() override; virtual void Deactivated() override; virtual void Visible() override; virtual void Hidden() override; protected slots: void DoFiberTracking(); ///< start fiber tracking void UpdateGui(); void ToggleInteractive(); void DeleteTrackingHandler(); void OnParameterChanged(); void InteractiveSeedChanged(bool posChanged=false); void ForestSwitched(); void OutputStyleSwitched(); void AfterThread(); ///< update gui etc. after tracking has finished void BeforeThread(); ///< start timer etc. void TimerUpdate(); void StopTractography(); void OnSliceChanged(); + void SaveParameters(); protected: /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkStreamlineTrackingViewControls* m_Controls; protected slots: private: bool CheckAndStoreLastParams(QObject* obj); void StartStopTrackingGui(bool start); + std::shared_ptr< mitk::StreamlineTractographyParameters > GetParametersFromGui(); std::vector< itk::Point > m_SeedPoints; mitk::DataNode::Pointer m_ParentNode; mitk::DataNode::Pointer m_InteractiveNode; mitk::DataNode::Pointer m_InteractivePointSetNode; std::vector< mitk::DataNode::Pointer > m_InputImageNodes; ///< input image nodes std::vector< mitk::Image::ConstPointer > m_AdditionalInputImages; bool m_FirstTensorProbRun; bool m_FirstInteractiveRun; mitk::TrackingDataHandler* m_TrackingHandler; bool m_ThreadIsRunning; QTimer* m_TrackingTimer; bool m_DeleteTrackingHandler; QmitkSliceNavigationListener m_SliceChangeListener; bool m_Visible; mitk::DataNode::Pointer m_LastPrior; mitk::TrackingDataHandler* m_TrackingPriorHandler; std::map< QString, std::string > m_LastTractoParams; + QString m_ParameterFile; }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui index 4ac04a83eb..119a845d75 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui @@ -1,1566 +1,1579 @@ QmitkStreamlineTrackingViewControls 0 0 453 859 0 0 QmitkTemplate QCommandLinkButton:disabled { border: none; } QGroupBox { background-color: transparent; } 3 3 0 40 QFrame::NoFrame QFrame::Raised 0 15 0 0 6 15 true 0 0 true QFrame::NoFrame QFrame::Raised 0 0 0 0 Input Image. ODF, tensor and peak images are currently supported. Input Image: Input Image. ODF, tensor, peak, and, in case of ML tractography, raw diffusion-weighted images are currently supported. <html><head/><body><p><span style=" color:#ff0000;">select image in data-manager</span></p></body></html> true Tractography Forest: true Stop tractography and return all fibers reconstructed until now. Stop Tractography false Start Tractography + + + + true + + + Stop tractography and return all fibers reconstructed until now. + + + Save Parameters + + + 0 0 0 0 75 true 0 0 0 421 254 Seeding Specify how, where and how many tractography seed points are placed. QFrame::NoFrame QFrame::Raised 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 Number of seed points equally distributed around selected position. 1 9999999 50 Radius: Seedpoints are equally distributed within a sphere centered at the selected position with the specified radius (in mm). 2 50.000000000000000 0.100000000000000 2.000000000000000 Num. Seeds: true When checked, parameter changes cause instant retracking while in interactive mode. Update on Parameter Change true QFrame::NoFrame QFrame::Raised 0 0 0 0 Try each seed N times until a valid streamline is obtained (only for probabilistic tractography). Minimum fiber length (in mm) 1 999 10 Trials Per Seed: Max. Num. Fibers: Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed (-1 means no limit). -1 999999999 -1 QFrame::NoFrame QFrame::Raised 0 0 0 0 Seeds per Voxel: Seed Image: Number of seed points placed in each voxel. 1 9999999 true Dynamically pick a seed location by click into image. Enable Interactive Tractography Qt::Vertical 20 40 0 0 435 232 ROI Constraints Specify various ROI and mask images to constrain the tractography process. Mask Image: Select which fibers should be accepted or rejected based on the location of their endpoints. QComboBox::AdjustToMinimumContentsLength No Constraints on EP locations Both EPs in Target Image Both EPs in Target Image But Different Label One EP in Seed Image and One EP in Target Image At Least One EP in Target Image Exactly One EP in Target Image No EP in Target Image Endpoint Constraints: Stop ROI Image: Exclusion ROI Image: Target ROI Image: Qt::Vertical 20 40 0 0 421 366 Tractography Parameters Specify the behavior of the tractography at each streamline integration step (step size, deterministic/probabilistic, ...). Angular threshold between two steps (in degree). Default: 90° * step_size -1 90 1 -1 FA/GFA Image: ODF Cutoff: Always produce the same random numbers. Qt::Vertical 20 40 f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 0.000000000000000 Minimum tract length in mm. Shorter fibers are discarded. Minimum fiber length (in mm) 1 999.000000000000000 1.000000000000000 20.000000000000000 Sharpen ODFs: Threshold on peak magnitude, FA, GFA, ... 5 1.000000000000000 0.100000000000000 0.100000000000000 g: Loop Check: 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 naturally much sharper. f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 1.000000000000000 Maximum allowed angular SDTEV over 4 voxel lengths. Default: no loop check. -1 180 -1 Min. Tract Length: f parameter of tensor tractography. f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). f: Angular Threshold: Step size (in voxels) 2 0.010000000000000 10.000000000000000 0.100000000000000 0.500000000000000 Toggle between deterministic and probabilistic tractography. Some modes might not be available for all types of tractography. Deterministic Probabilistic Mode: Step Size: Cutoff: Additional threshold on the ODF magnitude. This is useful in case of CSD fODF tractography. For fODFs a good default value is 0.1, for normalized dODFs, e.g. Q-ball ODFs, this threshold should be very low (0.00025) or 0. 5 1.000000000000000 0.100000000000000 0.000250000000000 Fix Random Seed: 0 0 435 232 Tractography Prior Weight: Weighting factor between prior and data. 1.000000000000000 0.100000000000000 0.500000000000000 Peak Image: Qt::Vertical 20 40 If unchecked, the prior cannot create directions where there are none in the data. true New Directions from Prior: Restrict to Prior: Restrict tractography to regions where the prior is valid. true QFrame::NoFrame QFrame::Raised 0 0 0 0 0 y x z Flip Directions: 0 0 435 232 Neighborhood Sampling Specify if and how information about the current streamline neighborhood should be used. Only neighborhood samples in front of the current streamline position are considered. Use Only Frontal Samples false If checked, the majority of sampling points has to place a stop-vote for the streamline to terminate. If not checked, all sampling positions have to vote for a streamline termination. Use Stop-Votes false QFrame::NoFrame QFrame::Raised 0 0 0 0 Num. Samples: Number of neighborhood samples that are used to determine the next fiber progression direction. 50 Sampling Distance: Sampling distance (in voxels) 2 10.000000000000000 0.100000000000000 0.250000000000000 Qt::Vertical 20 40 0 0 435 232 Data Handling Specify interpolation and direction flips. QFrame::NoFrame QFrame::Raised 0 0 0 0 Trilinearly interpolate the input image used for tractography. Interpolate Tractography Data true Trilinearly interpolate the ROI images used to constrain the tractography. Interpolate ROI Images true QFrame::NoFrame QFrame::Raised 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 Internally flips progression directions. This might be necessary depending on the input data. x Internally flips progression directions. This might be necessary depending on the input data. y Internally flips progression directions. This might be necessary depending on the input data. z Flip directions: Qt::Vertical 20 40 0 0 435 232 Output and Postprocessing Specify the tractography output (streamlines or probability maps) and postprocessing steps. QFrame::NoFrame QFrame::Raised 0 0 0 0 Compress fibers using the specified error constraint. Compress Fibers true Qt::StrongFocus Lossy fiber compression. Recommended for large tractograms. Maximum error in mm. 3 10.000000000000000 0.010000000000000 0.100000000000000 Output map with voxel-wise visitation counts instead of streamlines. Output Probability Map false Qt::Vertical 20 40 QmitkSingleNodeSelectionWidget QWidget
QmitkSingleNodeSelectionWidget.h
1 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/mitkPluginActivator.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/mitkPluginActivator.cpp index d5ad628635..9925065c85 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/mitkPluginActivator.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/mitkPluginActivator.cpp @@ -1,49 +1,47 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPluginActivator.h" #include "src/internal/Perspectives/QmitkGibbsTractographyPerspective.h" #include "src/internal/Perspectives/QmitkStreamlineTractographyPerspective.h" #include "src/internal/QmitkGibbsTrackingView.h" #include "src/internal/QmitkStreamlineTrackingView.h" -#include "src/internal/QmitkMLBTView.h" ctkPluginContext* mitk::PluginActivator::m_Context = nullptr; ctkPluginContext* mitk::PluginActivator::GetContext() { return m_Context; } void mitk::PluginActivator::start(ctkPluginContext* context) { BERRY_REGISTER_EXTENSION_CLASS(QmitkGibbsTractographyPerspective, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkStreamlineTractographyPerspective, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkGibbsTrackingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkStreamlineTrackingView, context) - BERRY_REGISTER_EXTENSION_CLASS(QmitkMLBTView, context) m_Context = context; } void mitk::PluginActivator::stop(ctkPluginContext* context) { Q_UNUSED(context) m_Context = nullptr; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimagingapp/resources/welcome/pics/3D_4_2.png b/Plugins/org.mitk.gui.qt.diffusionimagingapp/resources/welcome/pics/3D_4_2.png new file mode 100644 index 0000000000..b4ddce91d3 Binary files /dev/null and b/Plugins/org.mitk.gui.qt.diffusionimagingapp/resources/welcome/pics/3D_4_2.png differ