diff --git a/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp b/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp index 346a553..8bd29d8 100644 --- a/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp +++ b/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp @@ -1,134 +1,134 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, 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 "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include /*! \brief Perform global fiber tractography (Gibbs tractography) */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Gibbs Tracking"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform global fiber tractography (Gibbs tractography)"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image (tensor, ODF or SH-coefficient image)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output tractogram", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("parameters", "", mitkDiffusionCommandLineParser::String, "Parameters:", "parameter file (.gtp)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("mask", "", mitkDiffusionCommandLineParser::String, "Mask:", "binary mask image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string paramFileName = us::any_cast(parsedArgs["parameters"]); std::string outFileName = us::any_cast(parsedArgs["o"]); try { // instantiate gibbs tracker typedef itk::Vector OdfVectorType; typedef itk::Image ItkOdfImageType; typedef itk::GibbsTrackingFilter GibbsTrackingFilterType; GibbsTrackingFilterType::Pointer gibbsTracker = GibbsTrackingFilterType::New(); // load input image mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image"}, {}); mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(inFileName, &functor); // try to cast to Odf image if( dynamic_cast(mitkImage.GetPointer()) ) { mitk::OdfImage::Pointer mitkOdfImage = dynamic_cast(mitkImage.GetPointer()); ItkOdfImageType::Pointer itk_odf = ItkOdfImageType::New(); mitk::CastToItkImage(mitkOdfImage, itk_odf); gibbsTracker->SetOdfImage(itk_odf.GetPointer()); } else if( dynamic_cast(mitkImage.GetPointer()) ) { typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(mitkImage.GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); gibbsTracker->SetTensorImage(itk_dti); } else if ( dynamic_cast(mitkImage.GetPointer()) ) { MITK_INFO << "Assuming MITK/MRtrix style SH convention!"; mitk::Image::Pointer shImage = dynamic_cast(mitkImage.GetPointer()); gibbsTracker->SetOdfImage(mitk::convert::GetItkOdfFromShImage(shImage)); } else return EXIT_FAILURE; // global tracking if (parsedArgs.count("mask")) { typedef itk::Image MaskImgType; mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::Load(us::any_cast(parsedArgs["mask"])); MaskImgType::Pointer itk_mask = MaskImgType::New(); mitk::CastToItkImage(mitkMaskImage, itk_mask); gibbsTracker->SetMaskImage(itk_mask); } gibbsTracker->SetDuplicateImage(false); gibbsTracker->SetLoadParameterFile( paramFileName ); // gibbsTracker->SetLutPath( "" ); gibbsTracker->Update(); mitk::FiberBundle::Pointer mitkFiberBundle = mitk::FiberBundle::New(gibbsTracker->GetFiberBundle()); - mitkFiberBundle->SetReferenceGeometry(mitkImage->GetGeometry()); + mitkFiberBundle->SetTrackVisHeader(mitkImage->GetGeometry()); mitk::IOUtil::Save(mitkFiberBundle, outFileName ); } catch (const itk::ExceptionObject& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (std::exception& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp b/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp index e6daff8..ce1204b 100644 --- a/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp @@ -1,569 +1,577 @@ /*=================================================================== 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[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform streamline tractography"); parser.setContributor("MIC"); // parameters fo all methods parser.setArgumentPrefix("--", "-"); parser.beginGroup("1. Mandatory arguments:"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::StringList, "Input:", "input image (multiple possible for 'DetTensor' algorithm)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output fiberbundle/probability map", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("type", "", mitkDiffusionCommandLineParser::String, "Type:", "which tracker to use (Peaks; Tensor; ODF; ODF-DIPY/FSL; RF)", us::Any(), false); parser.addArgument("probabilistic", "", mitkDiffusionCommandLineParser::Bool, "Probabilistic:", "Probabilistic tractography", us::Any(false)); parser.endGroup(); parser.beginGroup("2. Seeding:"); parser.addArgument("seeds", "", mitkDiffusionCommandLineParser::Int, "Seeds per voxel:", "number of seed points per voxel", 1); parser.addArgument("seed_image", "", mitkDiffusionCommandLineParser::String, "Seed image:", "mask image defining seed voxels", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("trials_per_seed", "", mitkDiffusionCommandLineParser::Int, "Max. trials per seed:", "try each seed N times until a valid streamline is obtained (only for probabilistic tractography)", 10); parser.addArgument("max_tracts", "", mitkDiffusionCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded", -1); parser.endGroup(); parser.beginGroup("3. Tractography constraints:"); parser.addArgument("tracking_mask", "", mitkDiffusionCommandLineParser::String, "Mask image:", "streamlines leaving the mask will stop immediately", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("stop_image", "", mitkDiffusionCommandLineParser::String, "Stop ROI image:", "streamlines entering the mask will stop immediately", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("exclusion_image", "", mitkDiffusionCommandLineParser::String, "Exclusion ROI image:", "streamlines entering the mask will be discarded", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("ep_constraint", "", mitkDiffusionCommandLineParser::String, "Endpoint constraint:", "determines which fibers are accepted based on their endpoint location - options are NONE, EPS_IN_TARGET, EPS_IN_TARGET_LABELDIFF, EPS_IN_SEED_AND_TARGET, MIN_ONE_EP_IN_TARGET, ONE_EP_IN_TARGET and NO_EP_IN_TARGET", us::Any()); parser.addArgument("target_image", "", mitkDiffusionCommandLineParser::String, "Target ROI image:", "effact depends on the chosen endpoint constraint (option ep_constraint)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.endGroup(); parser.beginGroup("4. Streamline integration parameters:"); parser.addArgument("sharpen_odfs", "", mitkDiffusionCommandLineParser::Bool, "SHarpen ODFs:", "if you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). this is not necessary for CSD fODFs, since they are narurally much sharper."); parser.addArgument("cutoff", "", mitkDiffusionCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 0.1); parser.addArgument("odf_cutoff", "", mitkDiffusionCommandLineParser::Float, "ODF Cutoff:", "threshold on the ODF magnitude. this is useful in case of CSD fODF tractography.", 0.0); parser.addArgument("step_size", "", mitkDiffusionCommandLineParser::Float, "Step size:", "step size (in voxels)", 0.5); parser.addArgument("min_tract_length", "", mitkDiffusionCommandLineParser::Float, "Min. tract length:", "minimum fiber length (in mm)", 20); parser.addArgument("angular_threshold", "", mitkDiffusionCommandLineParser::Float, "Angular threshold:", "angular threshold between two successive steps, (default: 90° * step_size, minimum 15°)"); parser.addArgument("loop_check", "", mitkDiffusionCommandLineParser::Float, "Check for loops:", "threshold on angular stdev over the last 4 voxel lengths"); parser.addArgument("peak_jitter", "", mitkDiffusionCommandLineParser::Float, "Peak jitter:", "important for probabilistic peak tractography and peak prior. actual jitter is drawn from a normal distribution with peak_jitter*fabs(direction_value) as standard deviation.", 0.01); parser.endGroup(); parser.beginGroup("5. Tractography prior:"); parser.addArgument("prior_image", "", mitkDiffusionCommandLineParser::String, "Peak prior:", "tractography prior in thr for of a peak image", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("prior_weight", "", mitkDiffusionCommandLineParser::Float, "Prior weight", "weighting factor between prior and data.", 0.5); parser.addArgument("dont_restrict_to_prior", "", mitkDiffusionCommandLineParser::Bool, "Don't restrict to prior:", "don't restrict tractography to regions where the prior is valid.", us::Any(false)); parser.addArgument("no_new_directions_from_prior", "", mitkDiffusionCommandLineParser::Bool, "No new directios from prior:", "the prior cannot create directions where there are none in the data.", us::Any(false)); parser.addArgument("prior_flip_x", "", mitkDiffusionCommandLineParser::Bool, "Prior Flip X:", "multiply x-coordinate of prior direction by -1"); parser.addArgument("prior_flip_y", "", mitkDiffusionCommandLineParser::Bool, "Prior Flip Y:", "multiply y-coordinate of prior direction by -1"); parser.addArgument("prior_flip_z", "", mitkDiffusionCommandLineParser::Bool, "Prior Flip Z:", "multiply z-coordinate of prior direction by -1"); parser.endGroup(); parser.beginGroup("6. Neighborhood sampling:"); parser.addArgument("num_samples", "", mitkDiffusionCommandLineParser::Int, "Num. neighborhood samples:", "number of neighborhood samples that are use to determine the next progression direction", 0); parser.addArgument("sampling_distance", "", mitkDiffusionCommandLineParser::Float, "Sampling distance:", "distance of neighborhood sampling points (in voxels)", 0.25); parser.addArgument("use_stop_votes", "", mitkDiffusionCommandLineParser::Bool, "Use stop votes:", "use stop votes"); parser.addArgument("use_only_forward_samples", "", mitkDiffusionCommandLineParser::Bool, "Use only forward samples:", "use only forward samples"); parser.endGroup(); parser.beginGroup("7. Tensor tractography specific:"); parser.addArgument("tend_f", "", mitkDiffusionCommandLineParser::Float, "Weight f", "weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0); parser.addArgument("tend_g", "", mitkDiffusionCommandLineParser::Float, "Weight g", "weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0); parser.endGroup(); parser.beginGroup("8. Random forest tractography specific:"); parser.addArgument("forest", "", mitkDiffusionCommandLineParser::String, "Forest:", "input random forest (HDF5 file)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("use_sh_features", "", mitkDiffusionCommandLineParser::Bool, "Use SH features:", "use SH features"); parser.endGroup(); parser.beginGroup("9. Additional input:"); parser.addArgument("additional_images", "", mitkDiffusionCommandLineParser::StringList, "Additional images:", "specify a list of float images that hold additional information (FA, GFA, additional features for RF tractography)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.endGroup(); parser.beginGroup("10. Misc:"); parser.addArgument("flip_x", "", mitkDiffusionCommandLineParser::Bool, "Flip X:", "multiply x-coordinate of direction proposal by -1"); parser.addArgument("flip_y", "", mitkDiffusionCommandLineParser::Bool, "Flip Y:", "multiply y-coordinate of direction proposal by -1"); parser.addArgument("flip_z", "", mitkDiffusionCommandLineParser::Bool, "Flip Z:", "multiply z-coordinate of direction proposal by -1"); parser.addArgument("no_data_interpolation", "", mitkDiffusionCommandLineParser::Bool, "Don't interpolate input data:", "don't interpolate input image values"); parser.addArgument("no_mask_interpolation", "", mitkDiffusionCommandLineParser::Bool, "Don't interpolate masks:", "don't interpolate mask image values"); parser.addArgument("compress", "", mitkDiffusionCommandLineParser::Float, "Compress:", "compress output fibers using the given error threshold (in mm)"); parser.addArgument("fix_seed", "", mitkDiffusionCommandLineParser::Bool, "Fix Random Seed:", "always use the same random numbers"); parser.addArgument("parameter_file", "", mitkDiffusionCommandLineParser::String, "Parameter File:", "load parameters from json file (svae using MITK Diffusion GUI). the parameters loaded form this file are overwritten by the manually set parameters.", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.endGroup(); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["i"]); std::string outFile = us::any_cast(parsedArgs["o"]); std::string type = us::any_cast(parsedArgs["type"]); std::shared_ptr< mitk::StreamlineTractographyParameters > params = std::make_shared(); if (parsedArgs.count("parameter_file")) { auto parameter_file = us::any_cast(parsedArgs["parameter_file"]); params->LoadParameters(parameter_file); } if (parsedArgs.count("probabilistic")) params->m_Mode = mitk::StreamlineTractographyParameters::MODE::PROBABILISTIC; else { params->m_Mode = mitk::StreamlineTractographyParameters::MODE::DETERMINISTIC; } std::string prior_image = ""; if (parsedArgs.count("prior_image")) prior_image = us::any_cast(parsedArgs["prior_image"]); if (parsedArgs.count("prior_weight")) params->m_Weight = us::any_cast(parsedArgs["prior_weight"]); if (parsedArgs.count("fix_seed")) params->m_FixRandomSeed = us::any_cast(parsedArgs["fix_seed"]); params->m_RestrictToPrior = true; if (parsedArgs.count("dont_restrict_to_prior")) params->m_RestrictToPrior = !us::any_cast(parsedArgs["dont_restrict_to_prior"]); params->m_NewDirectionsFromPrior = true; if (parsedArgs.count("no_new_directions_from_prior")) params->m_NewDirectionsFromPrior = !us::any_cast(parsedArgs["no_new_directions_from_prior"]); params->m_SharpenOdfs = false; if (parsedArgs.count("sharpen_odfs")) params->m_SharpenOdfs = us::any_cast(parsedArgs["sharpen_odfs"]); params->m_InterpolateTractographyData = true; if (parsedArgs.count("no_data_interpolation")) params->m_InterpolateTractographyData = !us::any_cast(parsedArgs["no_data_interpolation"]); params->m_InterpolateRoiImages = true; if (parsedArgs.count("no_mask_interpolation")) params->m_InterpolateRoiImages = !us::any_cast(parsedArgs["no_mask_interpolation"]); bool use_sh_features = false; if (parsedArgs.count("use_sh_features")) use_sh_features = us::any_cast(parsedArgs["use_sh_features"]); params->m_StopVotes = false; if (parsedArgs.count("use_stop_votes")) params->m_StopVotes = us::any_cast(parsedArgs["use_stop_votes"]); params->m_OnlyForwardSamples = false; if (parsedArgs.count("use_only_forward_samples")) params->m_OnlyForwardSamples = us::any_cast(parsedArgs["use_only_forward_samples"]); params->m_FlipX = false; if (parsedArgs.count("flip_x")) params->m_FlipX = us::any_cast(parsedArgs["flip_x"]); params->m_FlipY = false; if (parsedArgs.count("flip_y")) params->m_FlipY = us::any_cast(parsedArgs["flip_y"]); params->m_FlipZ = false; if (parsedArgs.count("flip_z")) params->m_FlipZ = us::any_cast(parsedArgs["flip_z"]); bool prior_flip_x = false; if (parsedArgs.count("prior_flip_x")) prior_flip_x = us::any_cast(parsedArgs["prior_flip_x"]); bool prior_flip_y = false; if (parsedArgs.count("prior_flip_y")) prior_flip_y = us::any_cast(parsedArgs["prior_flip_y"]); bool prior_flip_z = false; if (parsedArgs.count("prior_flip_z")) prior_flip_z = us::any_cast(parsedArgs["prior_flip_z"]); params->m_ApplyDirectionMatrix = false; if (parsedArgs.count("apply_image_rotation")) params->m_ApplyDirectionMatrix = us::any_cast(parsedArgs["apply_image_rotation"]); float compress = -1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); params->m_MinTractLengthMm = 20; if (parsedArgs.count("min_tract_length")) params->m_MinTractLengthMm = us::any_cast(parsedArgs["min_tract_length"]); params->SetLoopCheckDeg(-1); if (parsedArgs.count("loop_check")) params->SetLoopCheckDeg(us::any_cast(parsedArgs["loop_check"])); std::string forestFile; if (parsedArgs.count("forest")) forestFile = us::any_cast(parsedArgs["forest"]); std::string maskFile = ""; if (parsedArgs.count("tracking_mask")) maskFile = us::any_cast(parsedArgs["tracking_mask"]); std::string seedFile = ""; if (parsedArgs.count("seed_image")) seedFile = us::any_cast(parsedArgs["seed_image"]); std::string targetFile = ""; if (parsedArgs.count("target_image")) targetFile = us::any_cast(parsedArgs["target_image"]); std::string exclusionFile = ""; if (parsedArgs.count("exclusion_image")) exclusionFile = us::any_cast(parsedArgs["exclusion_image"]); std::string stopFile = ""; if (parsedArgs.count("stop_image")) stopFile = us::any_cast(parsedArgs["stop_image"]); std::string ep_constraint = "NONE"; if (parsedArgs.count("ep_constraint")) ep_constraint = us::any_cast(parsedArgs["ep_constraint"]); params->m_Cutoff = 0.1f; if (parsedArgs.count("cutoff")) params->m_Cutoff = us::any_cast(parsedArgs["cutoff"]); params->m_OdfCutoff = 0.0; if (parsedArgs.count("odf_cutoff")) params->m_OdfCutoff = us::any_cast(parsedArgs["odf_cutoff"]); params->m_PeakJitter = 0.01; if (parsedArgs.count("peak_jitter")) params->m_PeakJitter = us::any_cast(parsedArgs["peak_jitter"]); params->SetStepSizeVox(-1); if (parsedArgs.count("step_size")) params->SetStepSizeVox(us::any_cast(parsedArgs["step_size"])); params->SetSamplingDistanceVox(-1); if (parsedArgs.count("sampling_distance")) params->SetSamplingDistanceVox(us::any_cast(parsedArgs["sampling_distance"])); params->m_NumSamples = 0; if (parsedArgs.count("num_samples")) params->m_NumSamples = static_cast(us::any_cast(parsedArgs["num_samples"])); params->m_SeedsPerVoxel = 1; if (parsedArgs.count("seeds")) params->m_SeedsPerVoxel = us::any_cast(parsedArgs["seeds"]); params->m_TrialsPerSeed = 10; if (parsedArgs.count("trials_per_seed")) params->m_TrialsPerSeed = static_cast(us::any_cast(parsedArgs["trials_per_seed"])); params->m_F = 1; if (parsedArgs.count("tend_f")) params->m_F = us::any_cast(parsedArgs["tend_f"]); params->m_G = 0; if (parsedArgs.count("tend_g")) params->m_G = us::any_cast(parsedArgs["tend_g"]); params->SetAngularThresholdDeg(-1); if (parsedArgs.count("angular_threshold")) params->SetAngularThresholdDeg(us::any_cast(parsedArgs["angular_threshold"])); params->m_MaxNumFibers = -1; if (parsedArgs.count("max_tracts")) params->m_MaxNumFibers = us::any_cast(parsedArgs["max_tracts"]); std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); if (ext != ".fib" && ext != ".trk") { MITK_INFO << "Output file format not supported. Use one of .fib, .trk, .nii, .nii.gz, .nrrd"; return EXIT_FAILURE; } // LOAD DATASETS mitkDiffusionCommandLineParser::StringContainerType addFiles; if (parsedArgs.count("additional_images")) addFiles = us::any_cast(parsedArgs["additional_images"]); typedef itk::Image ItkFloatImgType; ItkFloatImgType::Pointer mask = nullptr; if (!maskFile.empty()) { MITK_INFO << "loading mask image"; mitk::Image::Pointer img = mitk::IOUtil::Load(maskFile); mask = ItkFloatImgType::New(); mitk::CastToItkImage(img, mask); } ItkFloatImgType::Pointer seed = nullptr; if (!seedFile.empty()) { MITK_INFO << "loading seed ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(seedFile); seed = ItkFloatImgType::New(); mitk::CastToItkImage(img, seed); } ItkFloatImgType::Pointer stop = nullptr; if (!stopFile.empty()) { MITK_INFO << "loading stop ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(stopFile); stop = ItkFloatImgType::New(); mitk::CastToItkImage(img, stop); } ItkFloatImgType::Pointer target = nullptr; if (!targetFile.empty()) { MITK_INFO << "loading target ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(targetFile); target = ItkFloatImgType::New(); mitk::CastToItkImage(img, target); } ItkFloatImgType::Pointer exclusion = nullptr; if (!exclusionFile.empty()) { MITK_INFO << "loading exclusion ROI image"; mitk::Image::Pointer img = mitk::IOUtil::Load(exclusionFile); exclusion = ItkFloatImgType::New(); mitk::CastToItkImage(img, exclusion); } MITK_INFO << "loading additional images"; std::vector< std::vector< ItkFloatImgType::Pointer > > addImages; addImages.push_back(std::vector< ItkFloatImgType::Pointer >()); for (auto file : addFiles) { mitk::Image::Pointer img = mitk::IOUtil::Load(file); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); addImages.at(0).push_back(itkimg); } // ////////////////////////////////////////////////////////////////// // omp_set_num_threads(1); typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); if (!prior_image.empty()) { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image"}, {}); mitk::PeakImage::Pointer priorImage = mitk::IOUtil::Load(prior_image, &functor); if (priorImage.IsNull()) { MITK_INFO << "Only peak images are supported as prior at the moment!"; return EXIT_FAILURE; } mitk::TrackingDataHandler* priorhandler = new mitk::TrackingHandlerPeaks(); typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(priorImage); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); std::shared_ptr< mitk::StreamlineTractographyParameters > prior_params = std::make_shared< mitk::StreamlineTractographyParameters >(*params); prior_params->m_FlipX = prior_flip_x; prior_params->m_FlipY = prior_flip_y; prior_params->m_FlipZ = prior_flip_z; prior_params->m_Cutoff = 0.0; dynamic_cast(priorhandler)->SetPeakImage(itkImg); priorhandler->SetParameters(prior_params); tracker->SetTrackingPriorHandler(priorhandler); } mitk::TrackingDataHandler* handler; + mitk::Image::Pointer reference_image; if (type == "RF") { mitk::TractographyForest::Pointer forest = mitk::IOUtil::Load(forestFile); if (forest.IsNull()) mitkThrow() << "Forest file " << forestFile << " could not be read."; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images"}, {}); auto input = mitk::IOUtil::Load(input_files.at(0), &functor); + reference_image = input; if (use_sh_features) { handler = new mitk::TrackingHandlerRandomForest<6,28>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } else { handler = new mitk::TrackingHandlerRandomForest<6,100>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } } else if (type == "Peaks") { handler = new mitk::TrackingHandlerPeaks(); MITK_INFO << "loading input peak image"; mitk::Image::Pointer mitkImage = mitk::IOUtil::Load(input_files.at(0)); + reference_image = mitkImage; mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = mitk::convert::GetItkPeakFromPeakImage(mitkImage); dynamic_cast(handler)->SetPeakImage(itkImg); } else if (type == "Tensor" && params->m_Mode == mitk::StreamlineTractographyParameters::MODE::DETERMINISTIC) { handler = new mitk::TrackingHandlerTensor(); MITK_INFO << "loading input tensor images"; std::vector< mitk::Image::Pointer > input_images; for (unsigned int i=0; i(input_files.at(i)); + reference_image = mitkImage; mitk::TensorImage::ItkTensorImageType::Pointer itkImg = mitk::convert::GetItkTensorFromTensorImage(mitkImage); dynamic_cast(handler)->AddTensorImage(itkImg.GetPointer()); } if (addImages.at(0).size()>0) dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); } else if (type == "ODF" || type == "ODF-DIPY/FSL" || (type == "Tensor" && params->m_Mode == mitk::StreamlineTractographyParameters::MODE::PROBABILISTIC)) { handler = new mitk::TrackingHandlerOdf(); mitk::OdfImage::ItkOdfImageType::Pointer itkImg = nullptr; if (type == "Tensor") { MITK_INFO << "Converting Tensor to ODF image"; auto input = mitk::IOUtil::Load(input_files.at(0)); + reference_image = input; itkImg = mitk::convert::GetItkOdfFromTensorImage(input); dynamic_cast(handler)->SetIsOdfFromTensor(true); } else { mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"SH Image", "ODF Image"}, {}); auto input = mitk::IOUtil::Load(input_files.at(0), &functor)[0]; + reference_image = dynamic_cast(input.GetPointer()); if (dynamic_cast(input.GetPointer())) { MITK_INFO << "Converting SH to ODF image"; mitk::ShImage::Pointer mitkShImage = dynamic_cast(input.GetPointer()); if (type == "ODF-DIPY/FSL") mitkShImage->SetShConvention(mitk::ShImage::SH_CONVENTION::FSL); mitk::Image::Pointer mitkImg = dynamic_cast(mitkShImage.GetPointer()); itkImg = mitk::convert::GetItkOdfFromShImage(mitkImg); } else if (dynamic_cast(input.GetPointer())) { mitk::Image::Pointer mitkImg = dynamic_cast(input.GetPointer()); itkImg = mitk::convert::GetItkOdfFromOdfImage(mitkImg); } else mitkThrow() << ""; } dynamic_cast(handler)->SetOdfImage(itkImg); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); } else { MITK_INFO << "Unknown tractography algorithm (" + type+"). Known types are Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF."; return EXIT_FAILURE; } if (ep_constraint=="NONE") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::NONE; else if (ep_constraint=="EPS_IN_TARGET") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET; else if (ep_constraint=="EPS_IN_TARGET_LABELDIFF") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF; else if (ep_constraint=="EPS_IN_SEED_AND_TARGET") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET; else if (ep_constraint=="MIN_ONE_EP_IN_TARGET") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET; else if (ep_constraint=="ONE_EP_IN_TARGET") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET; else if (ep_constraint=="NO_EP_IN_TARGET") params->m_EpConstraints = itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET; MITK_INFO << "Tractography algorithm: " << type; tracker->SetMaskImage(mask); tracker->SetSeedImage(seed); tracker->SetStoppingRegions(stop); tracker->SetTargetRegions(target); tracker->SetExclusionRegions(exclusion); tracker->SetTrackingHandler(handler); if (ext != ".fib" && ext != ".trk") params->m_OutputProbMap = true; tracker->SetParameters(params); tracker->Update(); if (ext == ".fib" || ext == ".trk") { vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); if (compress > 0) outFib->Compress(compress); + outFib->SetTrackVisHeader(reference_image->GetGeometry()); + mitk::IOUtil::Save(outFib, outFile); } else { TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (ext != ".nii" && ext != ".nii.gz" && ext != ".nrrd") outFile += ".nii.gz"; mitk::IOUtil::Save(img, outFile); } delete handler; return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp index d7794bb..3c27ede 100644 --- a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp +++ b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.cpp @@ -1,2705 +1,2756 @@ /*=================================================================== 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 "mitkFiberBundle.h" #include #include #include #include "mitkImagePixelReadAccessor.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs"; mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData ) : m_NumFibers(0) { m_TrackVisHeader.hdr_size = 0; m_FiberWeights = vtkSmartPointer::New(); m_FiberWeights->SetName("FIBER_WEIGHTS"); m_FiberPolyData = vtkSmartPointer::New(); if (fiberPolyData != nullptr) m_FiberPolyData = fiberPolyData; else { this->m_FiberPolyData->SetPoints(vtkSmartPointer::New()); this->m_FiberPolyData->SetLines(vtkSmartPointer::New()); } this->UpdateFiberGeometry(); this->GenerateFiberIds(); this->ColorFibersByOrientation(); } mitk::FiberBundle::~FiberBundle() { } mitk::FiberBundle::Pointer mitk::FiberBundle::GetDeepCopy() { mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData); newFib->SetFiberColors(this->m_FiberColors); newFib->SetFiberWeights(this->m_FiberWeights); - newFib->SetReferenceGeometry(this->GetReferenceGeometry()); + newFib->SetTrackVisHeader(this->GetTrackVisHeader()); return newFib; } vtkSmartPointer mitk::FiberBundle::GeneratePolyDataByIds(std::vector fiberIds, vtkSmartPointer weights) { vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); vtkSmartPointer newLineSet = vtkSmartPointer::New(); vtkSmartPointer newPointSet = vtkSmartPointer::New(); weights->SetNumberOfValues(fiberIds.size()); int counter = 0; auto finIt = fiberIds.begin(); while ( finIt != fiberIds.end() ) { if (*finIt>GetNumFibers()){ MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt; break; } vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); vtkSmartPointer fibPoints = fiber->GetPoints(); vtkSmartPointer newFiber = vtkSmartPointer::New(); newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); for(int i=0; iGetNumberOfPoints(); i++) { newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); } weights->InsertValue(counter, this->GetFiberWeight(*finIt)); newLineSet->InsertNextCell(newFiber); ++finIt; ++counter; } newFiberPolyData->SetPoints(newPointSet); newFiberPolyData->SetLines(newLineSet); return newFiberPolyData; } // merge two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundles(std::vector< mitk::FiberBundle::Pointer > fibs) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); auto num_weights = this->GetNumFibers(); for (auto fib : fibs) num_weights += fib->GetNumFibers(); weights->SetNumberOfValues(num_weights); unsigned int counter = 0; for (unsigned int i=0; iGetNumberOfCells(); ++i) { vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (unsigned int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, this->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } for (auto fib : fibs) { // add new fiber bundle for (unsigned int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (unsigned int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, fib->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); newFib->SetFiberWeights(weights); return newFib; } // merge two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib) { if (fib==nullptr) return this->GetDeepCopy(); MITK_INFO << "Adding fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers()); unsigned int counter = 0; for (unsigned int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (unsigned int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, this->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } // add new fiber bundle for (unsigned int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (unsigned int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, fib->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); newFib->SetFiberWeights(weights); return newFib; } // Only retain fibers with a weight larger than the specified threshold mitk::FiberBundle::Pointer mitk::FiberBundle::FilterByWeights(float weight_thr, bool invert) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); std::vector weights; for (unsigned int i=0; iGetNumFibers(); i++) { if ( (invert && this->GetFiberWeight(i)>weight_thr) || (!invert && this->GetFiberWeight(i)<=weight_thr)) continue; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); weights.push_back(this->GetFiberWeight(i)); } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); for (unsigned int i=0; iSetFiberWeight(i, weights.at(i)); return newFib; } // Only retain a subsample of the fibers mitk::FiberBundle::Pointer mitk::FiberBundle::SubsampleFibers(float factor, bool random_seed) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); unsigned int new_num_fibs = static_cast(std::round(this->GetNumFibers()*factor)); MITK_INFO << "Subsampling fibers with factor " << factor << "(" << new_num_fibs << "/" << this->GetNumFibers() << ")"; // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); weights->SetNumberOfValues(new_num_fibs); std::vector< unsigned int > ids; for (unsigned int i=0; iGetNumFibers(); i++) ids.push_back(i); if (random_seed) std::srand(static_cast(std::time(nullptr))); else std::srand(0); std::random_shuffle(ids.begin(), ids.end()); unsigned int counter = 0; for (unsigned int i=0; iGetCell(ids.at(i)); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, this->GetFiberWeight(ids.at(i))); vNewLines->InsertNextCell(container); counter++; } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); newFib->SetFiberWeights(weights); return newFib; } // subtract two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib) { if (fib==nullptr) return this->GetDeepCopy(); MITK_INFO << "Subtracting fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); std::vector< std::vector< itk::Point > > points1; for(unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; itk::Point start = mitk::imv::GetItkPoint(points->GetPoint(0)); itk::Point end = mitk::imv::GetItkPoint(points->GetPoint(numPoints-1)); points1.push_back( {start, end} ); } std::vector< std::vector< itk::Point > > points2; for(unsigned int i=0; iGetNumFibers(); i++ ) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; itk::Point start =mitk::imv::GetItkPoint(points->GetPoint(0)); itk::Point end =mitk::imv::GetItkPoint(points->GetPoint(numPoints-1)); points2.push_back( {start, end} ); } // int progress = 0; std::vector< int > ids; #pragma omp parallel for for (int i=0; i(points1.size()); i++) { bool match = false; for (unsigned int j=0; j(i)); auto v2 = points2.at(j); float dist=0; for (unsigned int c=0; cGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(points->GetPoint(j)); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } if(vNewLines->GetNumberOfCells()==0) return mitk::FiberBundle::New(); // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle return mitk::FiberBundle::New(vNewPolyData); } /* * set PolyData (additional flag to recompute fiber geometry, default = true) */ void mitk::FiberBundle::SetFiberPolyData(vtkSmartPointer fiberPD, bool updateGeometry) { if (fiberPD == nullptr) this->m_FiberPolyData = vtkSmartPointer::New(); else m_FiberPolyData->DeepCopy(fiberPD); m_NumFibers = static_cast(m_FiberPolyData->GetNumberOfLines()); if (updateGeometry) UpdateFiberGeometry(); GenerateFiberIds(); ColorFibersByOrientation(); } /* * return vtkPolyData */ vtkSmartPointer mitk::FiberBundle::GetFiberPolyData() const { return m_FiberPolyData; } void mitk::FiberBundle::ColorFibersByLength(bool opacity, bool normalize) { if (m_MaxFiberLength<=0) return; auto numOfPoints = this->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(numOfPoints * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); auto numOfFibers = m_FiberPolyData->GetNumberOfLines(); if (numOfFibers < 1) return; mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); unsigned int count = 0; for (unsigned int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); float l = m_FiberLengths.at(i)/m_MaxFiberLength; if (!normalize) { l = m_FiberLengths.at(i)/255.0f; if (l > 1.0f) l = 1.0; } for (int j=0; jGetColor(1.0 - static_cast(l), color); rgba[0] = static_cast(255.0 * color[0]); rgba[1] = static_cast(255.0 * color[1]); rgba[2] = static_cast(255.0 * color[2]); if (opacity) rgba[3] = static_cast(255.0f * l); else rgba[3] = static_cast(255.0); m_FiberColors->InsertTypedTuple(cell->GetPointId(j), rgba); count++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByOrientation() { //===== FOR WRITING A TEST ======================== // colorT size == tupelComponents * tupelElements // compare color results // to cover this code 100% also PolyData needed, where colorarray already exists // + one fiber with exactly 1 point // + one fiber with 0 points //================================================= vtkPoints* extrPoints = m_FiberPolyData->GetPoints(); vtkIdType numOfPoints = 0; if (extrPoints!=nullptr) numOfPoints = extrPoints->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(numOfPoints * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); auto numOfFibers = m_FiberPolyData->GetNumberOfLines(); if (numOfFibers < 1) return; /* extract single fibers of fiberBundle */ vtkCellArray* fiberList = m_FiberPolyData->GetLines(); fiberList->InitTraversal(); for (int fi=0; fiGetNextCell(pointsPerFiber, idList); /* single fiber checkpoints: is number of points valid */ if (pointsPerFiber > 1) { /* operate on points of single fiber */ for (int i=0; i 0) { /* The color value of the current point is influenced by the previous point and next point. */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; vnl_vector_fixed< double, 3 > diff; diff = (diff1 - diff2) / 2.0; diff.normalize(); rgba[0] = static_cast(255.0 * std::fabs(diff[0])); rgba[1] = static_cast(255.0 * std::fabs(diff[1])); rgba[2] = static_cast(255.0 * std::fabs(diff[2])); rgba[3] = static_cast(255.0); } else if (i==0) { /* First point has no previous point, therefore only diff1 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; diff1.normalize(); rgba[0] = static_cast(255.0 * std::fabs(diff1[0])); rgba[1] = static_cast(255.0 * std::fabs(diff1[1])); rgba[2] = static_cast(255.0 * std::fabs(diff1[2])); rgba[3] = static_cast(255.0); } else if (i==pointsPerFiber-1) { /* Last point has no next point, therefore only diff2 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; diff2.normalize(); rgba[0] = static_cast(255.0 * std::fabs(diff2[0])); rgba[1] = static_cast(255.0 * std::fabs(diff2[1])); rgba[2] = static_cast(255.0 * std::fabs(diff2[2])); rgba[3] = static_cast(255.0); } m_FiberColors->InsertTypedTuple(idList[i], rgba); } } else if (pointsPerFiber == 1) { /* a single point does not define a fiber (use vertex mechanisms instead */ continue; } else { MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; continue; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByCurvature(bool, bool normalize) { double window = 5; //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); std::vector< double > values; double min = 1; double max = 0; MITK_INFO << "Coloring fibers by curvature"; boost::progress_display disp(static_cast(m_FiberPolyData->GetNumberOfCells())); for (int i=0; iGetNumberOfCells(); i++) { ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures for (int j=0; j > vectors; vnl_vector_fixed< double, 3 > meanV; meanV.fill(0.0); while(dist1) { double p1[3]; points->GetPoint(c-1, p1); double p2[3]; points->GetPoint(c, p2); vnl_vector_fixed< double, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); meanV += v; c--; } c = j; dist = 0; while(distGetPoint(c, p1); double p2[3]; points->GetPoint(c+1, p2); vnl_vector_fixed< double, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); meanV += v; c++; } meanV.normalize(); double dev = 0; for (unsigned int c=0; c1.0) angle = 1.0; if (angle<-1.0) angle = -1.0; dev += acos(angle)*180/itk::Math::pi; } if (vectors.size()>0) dev /= vectors.size(); dev = 1.0-dev/180.0; values.push_back(dev); if (devmax) max = dev; } } unsigned int count = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); for (int j=0; j1) dev = 1; lookupTable->GetColor(dev, color); rgba[0] = static_cast(255.0 * color[0]); rgba[1] = static_cast(255.0 * color[1]); rgba[2] = static_cast(255.0 * color[2]); rgba[3] = static_cast(255.0); m_FiberColors->InsertTypedTuple(cell->GetPointId(j), rgba); count++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberOpacity(vtkDoubleArray* FAValArray) { for(long i=0; iGetNumberOfTuples(); i++) { double faValue = FAValArray->GetValue(i); faValue = faValue * 255.0; m_FiberColors->SetComponent(i,3, static_cast(faValue) ); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ResetFiberOpacity() { for(long i=0; iGetNumberOfTuples(); i++) m_FiberColors->SetComponent(i,3, 255.0 ); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByScalarMap(mitk::Image::Pointer FAimage, bool opacity, bool normalize) { mitkPixelTypeMultiplex3( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity, normalize ); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } template void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity, bool normalize) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::ImagePixelReadAccessor readimage(image, image->GetVolumeData(0)); unsigned char rgba[4] = {0,0,0,0}; vtkPoints* pointSet = m_FiberPolyData->GetPoints(); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); double min = 999999; double max = -999999; for(long i=0; iGetNumberOfPoints(); ++i) { Point3D px; px[0] = pointSet->GetPoint(i)[0]; px[1] = pointSet->GetPoint(i)[1]; px[2] = pointSet->GetPoint(i)[2]; auto pixelValue = static_cast(readimage.GetPixelByWorldCoordinates(px)); if (pixelValue>max) max = pixelValue; if (pixelValueGetNumberOfPoints(); ++i) { Point3D px; px[0] = pointSet->GetPoint(i)[0]; px[1] = pointSet->GetPoint(i)[1]; px[2] = pointSet->GetPoint(i)[2]; auto pixelValue = static_cast(readimage.GetPixelByWorldCoordinates(px)); if (normalize) pixelValue = (pixelValue-min)/(max-min); else if (pixelValue>1) pixelValue = 1; double color[3]; lookupTable->GetColor(1-pixelValue, color); rgba[0] = static_cast(255.0 * color[0]); rgba[1] = static_cast(255.0 * color[1]); rgba[2] = static_cast(255.0 * color[2]); if (opacity) rgba[3] = static_cast(255.0 * pixelValue); else rgba[3] = static_cast(255.0); m_FiberColors->InsertTypedTuple(i, rgba); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByFiberWeights(bool opacity, bool normalize) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); unsigned char rgba[4] = {0,0,0,0}; unsigned int counter = 0; float max = -999999; float min = 999999; for (unsigned int i=0; iGetFiberWeight(i); if (weight>max) max = weight; if (weightGetCell(i); auto numPoints = cell->GetNumberOfPoints(); auto weight = this->GetFiberWeight(i); for (int j=0; j1) v = 1; double color[3]; lookupTable->GetColor(static_cast(1-v), color); rgba[0] = static_cast(255.0 * color[0]); rgba[1] = static_cast(255.0 * color[1]); rgba[2] = static_cast(255.0 * color[2]); if (opacity) rgba[3] = static_cast(255.0f * v); else rgba[3] = static_cast(255.0); m_FiberColors->InsertTypedTuple(counter, rgba); counter++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberColors(float r, float g, float b, float alpha) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); unsigned char rgba[4] = {0,0,0,0}; for(long i=0; iGetNumberOfPoints(); ++i) { rgba[0] = static_cast(r); rgba[1] = static_cast(g); rgba[2] = static_cast(b); rgba[3] = static_cast(alpha); m_FiberColors->InsertTypedTuple(i, rgba); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::GenerateFiberIds() { if (m_FiberPolyData == nullptr) return; vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); idFiberFilter->SetInputData(m_FiberPolyData); idFiberFilter->CellIdsOn(); // idFiberFilter->PointIdsOn(); // point id's are not needed idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); idFiberFilter->FieldDataOn(); idFiberFilter->Update(); m_FiberIdDataSet = idFiberFilter->GetOutput(); } float mitk::FiberBundle::GetNumEpFractionInMask(ItkUcharImgType* mask, bool different_label) { vtkSmartPointer PolyData = m_FiberPolyData; MITK_INFO << "Calculating EP-Fraction"; boost::progress_display disp(m_NumFibers); unsigned int in_mask = 0; for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); itk::Point startVertex =mitk::imv::GetItkPoint(points->GetPoint(0)); itk::Index<3> startIndex; mask->TransformPhysicalPointToIndex(startVertex, startIndex); itk::Point endVertex =mitk::imv::GetItkPoint(points->GetPoint(numPoints-1)); itk::Index<3> endIndex; mask->TransformPhysicalPointToIndex(endVertex, endIndex); if (mask->GetLargestPossibleRegion().IsInside(startIndex) && mask->GetLargestPossibleRegion().IsInside(endIndex)) { float v1 = mask->GetPixel(startIndex); if (v1 < 0.5f) continue; float v2 = mask->GetPixel(startIndex); if (v2 < 0.5f) continue; if (!different_label) ++in_mask; else if (fabs(v1-v2)>0.00001f) ++in_mask; } } return float(in_mask)/m_NumFibers; } std::tuple mitk::FiberBundle::GetDirectionalOverlap(ItkUcharImgType* mask, mitk::PeakImage::ItkPeakImageType* peak_image) { vtkSmartPointer PolyData = m_FiberPolyData; MITK_INFO << "Calculating overlap"; auto spacing = mask->GetSpacing(); boost::progress_display disp(m_NumFibers); double length_sum = 0; double in_mask_length = 0; double aligned_length = 0; for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; j startVertex =mitk::imv::GetItkPoint(points->GetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; mask->TransformPhysicalPointToIndex(startVertex, startIndex); mask->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); itk::Point endVertex =mitk::imv::GetItkPoint(points->GetPoint(j + 1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; mask->TransformPhysicalPointToIndex(endVertex, endIndex); mask->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); vnl_vector_fixed< float, 3 > fdir; fdir[0] = endVertex[0] - startVertex[0]; fdir[1] = endVertex[1] - startVertex[1]; fdir[2] = endVertex[2] - startVertex[2]; fdir.normalize(); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > segment : segments) { if ( mask->GetLargestPossibleRegion().IsInside(segment.first) && mask->GetPixel(segment.first) > 0 ) { in_mask_length += segment.second; mitk::PeakImage::ItkPeakImageType::IndexType idx4; idx4[0] = segment.first[0]; idx4[1] = segment.first[1]; idx4[2] = segment.first[2]; vnl_vector_fixed< float, 3 > peak; idx4[3] = 0; peak[0] = peak_image->GetPixel(idx4); idx4[3] = 1; peak[1] = peak_image->GetPixel(idx4); idx4[3] = 2; peak[2] = peak_image->GetPixel(idx4); if (std::isnan(peak[0]) || std::isnan(peak[1]) || std::isnan(peak[2]) || peak.magnitude()<0.0001f) continue; peak.normalize(); double f = 1.0 - std::acos(std::fabs(static_cast(dot_product(fdir, peak)))) * 2.0/itk::Math::pi; aligned_length += segment.second * f; } length_sum += segment.second; } } } if (length_sum<=0.0001) { MITK_INFO << "Fiber length sum is zero!"; return std::make_tuple(0,0); } return std::make_tuple(aligned_length/length_sum, in_mask_length/length_sum); } float mitk::FiberBundle::GetOverlap(ItkUcharImgType* mask) { vtkSmartPointer PolyData = m_FiberPolyData; MITK_INFO << "Calculating overlap"; auto spacing = mask->GetSpacing(); boost::progress_display disp(m_NumFibers); double length_sum = 0; double in_mask_length = 0; for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; j startVertex =mitk::imv::GetItkPoint(points->GetPoint(j)); itk::Index<3> startIndex; itk::ContinuousIndex startIndexCont; mask->TransformPhysicalPointToIndex(startVertex, startIndex); mask->TransformPhysicalPointToContinuousIndex(startVertex, startIndexCont); itk::Point endVertex =mitk::imv::GetItkPoint(points->GetPoint(j + 1)); itk::Index<3> endIndex; itk::ContinuousIndex endIndexCont; mask->TransformPhysicalPointToIndex(endVertex, endIndex); mask->TransformPhysicalPointToContinuousIndex(endVertex, endIndexCont); std::vector< std::pair< itk::Index<3>, double > > segments = mitk::imv::IntersectImage(spacing, startIndex, endIndex, startIndexCont, endIndexCont); for (std::pair< itk::Index<3>, double > segment : segments) { if ( mask->GetLargestPossibleRegion().IsInside(segment.first) && mask->GetPixel(segment.first) > 0 ) in_mask_length += segment.second; length_sum += segment.second; } } } if (length_sum<=0.000001) { MITK_INFO << "Fiber length sum is zero!"; return 0; } return static_cast(in_mask_length/length_sum); } mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); std::vector< float > fib_weights; MITK_INFO << "Cutting fibers"; boost::progress_display disp(m_NumFibers); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); int newNumPoints = 0; if (numPoints>1) { for (int j=0; j itkP =mitk::imv::GetItkPoint(points->GetPoint(j)); itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); bool inside = false; if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx)!=0 ) inside = true; if (inside && !invert) { vtkIdType id = vtkNewPoints->InsertNextPoint(itkP.GetDataPointer()); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if ( !inside && invert ) { vtkIdType id = vtkNewPoints->InsertNextPoint(itkP.GetDataPointer()); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if (newNumPoints>1) { fib_weights.push_back(this->GetFiberWeight(i)); vtkNewCells->InsertNextCell(container); newNumPoints = 0; container = vtkSmartPointer::New(); } else { newNumPoints = 0; container = vtkSmartPointer::New(); } } if (newNumPoints>1) { fib_weights.push_back(this->GetFiberWeight(i)); vtkNewCells->InsertNextCell(container); } } } vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(static_cast(fib_weights.size())); if (vtkNewCells->GetNumberOfCells()<=0) return nullptr; for (unsigned int i=0; iGetNumberOfValues(); i++) newFiberWeights->SetValue(i, fib_weights.at(i)); // vtkSmartPointer newFiberColors = vtkSmartPointer::New(); // newFiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); // newFiberColors->SetNumberOfComponents(4); // newFiberColors->SetName("FIBER_COLORS"); // unsigned char rgba[4] = {0,0,0,0}; // for(long i=0; iGetNumberOfPoints(); ++i) // { // rgba[0] = (unsigned char) r; // rgba[1] = (unsigned char) g; // rgba[2] = (unsigned char) b; // rgba[3] = (unsigned char) alpha; // m_FiberColors->InsertTypedTuple(i, rgba); // } vtkSmartPointer newPolyData = vtkSmartPointer::New(); newPolyData->SetPoints(vtkNewPoints); newPolyData->SetLines(vtkNewCells); mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData); newFib->SetFiberWeights(newFiberWeights); // newFib->Compress(0.1); return newFib; } mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(DataNode* roi, DataStorage* storage) { if (roi==nullptr || !(dynamic_cast(roi->GetData()) || dynamic_cast(roi->GetData())) ) return nullptr; std::vector tmp = ExtractFiberIdSubset(roi, storage); if (tmp.size()<=0) return mitk::FiberBundle::New(); vtkSmartPointer weights = vtkSmartPointer::New(); vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp, weights); mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); fib->SetFiberWeights(weights); return fib; } std::vector mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage) { std::vector result; if (roi==nullptr || roi->GetData()==nullptr) return result; mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast(roi->GetData()); if (!pfc.IsNull()) // handle composite { DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi); if (children->size()==0) return result; switch (pfc->getOperationType()) { case 0: // AND { MITK_INFO << "AND"; result = this->ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { std::vector inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(std::min(result.size(),inRoi.size())); it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( static_cast(it - rest.begin()) ); result = rest; } break; } case 1: // OR { MITK_INFO << "OR"; result = ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { it = result.end(); std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); result.insert(it, inRoi.begin(), inRoi.end()); } // remove duplicates sort(result.begin(), result.end()); it = unique(result.begin(), result.end()); result.resize( static_cast(it - result.begin()) ); break; } case 2: // NOT { MITK_INFO << "NOT"; for(unsigned int i=0; iGetNumFibers(); i++) result.push_back(i); std::vector::iterator it; for (unsigned int i=0; iSize(); ++i) { std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(result.size()-inRoi.size()); it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( static_cast(it - rest.begin()) ); result = rest; } break; } } } else if ( dynamic_cast(roi->GetData()) ) // actual extraction { if ( dynamic_cast(roi->GetData()) ) { mitk::PlanarFigure::Pointer planarPoly = dynamic_cast(roi->GetData()); //create vtkPolygon using controlpoints from planarFigure polygon vtkSmartPointer polygonVtk = vtkSmartPointer::New(); for (unsigned int i=0; iGetNumberOfControlPoints(); ++i) { itk::Point p = planarPoly->GetWorldControlPoint(i); vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] ); polygonVtk->GetPointIds()->InsertNextId(id); } MITK_INFO << "Extracting with polygon"; boost::progress_display disp(m_NumFibers); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p1); double p2[3] = {0,0,0}; points->GetPoint(j+1, p2); double tolerance = 0.001; // Outputs double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) double x[3] = {0,0,0}; // The coordinate of the intersection double pcoords[3] = {0,0,0}; int subId = 0; int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId); if (iD!=0) { result.push_back(i); break; } } } } else if ( dynamic_cast(roi->GetData()) ) { mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(roi->GetData()); Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal(); planeNormal.Normalize(); //calculate circle radius mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint mitk::Point3D V2w = planarFigure->GetWorldControlPoint(1); //radiusPoint double radius = V1w.EuclideanDistanceTo(V2w); radius *= radius; MITK_INFO << "Extracting with circle"; boost::progress_display disp(m_NumFibers); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p1); double p2[3] = {0,0,0}; points->GetPoint(j+1, p2); // Outputs double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) double x[3] = {0,0,0}; // The coordinate of the intersection int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x); if (iD!=0) { double dist = (x[0]-V1w[0])*(x[0]-V1w[0])+(x[1]-V1w[1])*(x[1]-V1w[1])+(x[2]-V1w[2])*(x[2]-V1w[2]); if( dist <= radius) { result.push_back(i); break; } } } } } return result; } return result; } void mitk::FiberBundle::UpdateFiberGeometry() { vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInputData(m_FiberPolyData); cleaner->PointMergingOff(); cleaner->Update(); m_FiberPolyData = cleaner->GetOutput(); m_FiberLengths.clear(); m_MeanFiberLength = 0; m_MedianFiberLength = 0; m_LengthStDev = 0; m_NumFibers = static_cast(m_FiberPolyData->GetNumberOfCells()); if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints()) this->ColorFibersByOrientation(); if (m_FiberWeights->GetNumberOfValues()!=m_NumFibers) { m_FiberWeights = vtkSmartPointer::New(); m_FiberWeights->SetName("FIBER_WEIGHTS"); m_FiberWeights->SetNumberOfValues(m_NumFibers); this->SetFiberWeights(1); } if (m_NumFibers<=0) // no fibers present; apply default geometry { m_MinFiberLength = 0; m_MaxFiberLength = 0; mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetImageGeometry(false); float b[] = {0, 1, 0, 1, 0, 1}; geometry->SetFloatBounds(b); SetGeometry(geometry); return; } double b[6]; m_FiberPolyData->GetBounds(b); // calculate statistics for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); auto p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float length = 0; for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2])); length += static_cast(dist); } m_FiberLengths.push_back(length); m_MeanFiberLength += length; if (i==0) { m_MinFiberLength = length; m_MaxFiberLength = length; } else { if (lengthm_MaxFiberLength) m_MaxFiberLength = length; } } m_MeanFiberLength /= m_NumFibers; std::vector< float > sortedLengths = m_FiberLengths; std::sort(sortedLengths.begin(), sortedLengths.end()); for (unsigned int i=0; i1) m_LengthStDev /= (m_NumFibers-1); else m_LengthStDev = 0; m_LengthStDev = std::sqrt(m_LengthStDev); m_MedianFiberLength = sortedLengths.at(m_NumFibers/2); mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetFloatBounds(b); this->SetGeometry(geometry); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } float mitk::FiberBundle::GetFiberWeight(unsigned int fiber) const { return m_FiberWeights->GetValue(fiber); } void mitk::FiberBundle::SetFiberWeights(float newWeight) { for (int i=0; iGetNumberOfValues(); i++) m_FiberWeights->SetValue(i, newWeight); } void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer weights) { if (m_NumFibers!=weights->GetNumberOfValues()) { MITK_INFO << "Weights array not equal to number of fibers! " << weights->GetNumberOfValues() << " vs " << m_NumFibers; return; } for (int i=0; iGetNumberOfValues(); i++) m_FiberWeights->SetValue(i, weights->GetValue(i)); m_FiberWeights->SetName("FIBER_WEIGHTS"); } void mitk::FiberBundle::SetFiberWeight(unsigned int fiber, float weight) { m_FiberWeights->SetValue(fiber, weight); } void mitk::FiberBundle::SetFiberColors(vtkSmartPointer fiberColors) { for(long i=0; iGetNumberOfPoints(); ++i) { unsigned char source[4] = {0,0,0,0}; fiberColors->GetTypedTuple(i, source); unsigned char target[4] = {0,0,0,0}; target[0] = source[0]; target[1] = source[1]; target[2] = source[2]; target[3] = source[3]; m_FiberColors->InsertTypedTuple(i, target); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } itk::Matrix< double, 3, 3 > mitk::FiberBundle::TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz) { rx = rx*itk::Math::pi/180; ry = ry*itk::Math::pi/180; rz = rz*itk::Math::pi/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; m = rot*m; return m; } void mitk::FiberBundle::TransformFibers(itk::ScalableAffineTransform< mitk::ScalarType >::Pointer transform) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; j p =mitk::imv::GetItkPoint(points->GetPoint(j)); p = transform->TransformPoint(p); vtkIdType id = vtkNewPoints->InsertNextPoint(p.GetDataPointer()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz) { vnl_matrix_fixed< double, 3, 3 > rot = mitk::imv::GetRotationMatrixVnl(rx, ry, rz); mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vnl_vector_fixed< double, 3 > dir; dir[0] = p[0]-center[0]; dir[1] = p[1]-center[1]; dir[2] = p[2]-center[2]; dir = rot*dir; dir[0] += center[0]+tx; dir[1] += center[1]+ty; dir[2] += center[2]+tz; vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RotateAroundAxis(double x, double y, double z) { x = x*itk::Math::pi/180; y = y*itk::Math::pi/180; z = z*itk::Math::pi/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vnl_vector_fixed< double, 3 > dir; dir[0] = p[0]-center[0]; dir[1] = p[1]-center[1]; dir[2] = p[2]-center[2]; dir = rotZ*rotY*rotX*dir; dir[0] += center[0]; dir[1] += center[1]; dir[2] += center[2]; vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ScaleFibers(double x, double y, double z, bool subtractCenter) { MITK_INFO << "Scaling fibers"; boost::progress_display disp(m_NumFibers); mitk::BaseGeometry* geom = this->GetGeometry(); mitk::Point3D c = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); if (subtractCenter) { p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2]; } p[0] *= x; p[1] *= y; p[2] *= z; if (subtractCenter) { p[0] += c[0]; p[1] += c[1]; p[2] += c[2]; } vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::TranslateFibers(double x, double y, double z) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[0] += x; p[1] += y; p[2] += z; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::MirrorFibers(unsigned int axis) { if (axis>2) return; MITK_INFO << "Mirroring fibers"; boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[axis] = -p[axis]; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RemoveDir(vnl_vector_fixed dir, double threshold) { dir.normalize(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); boost::progress_display disp(static_cast(m_FiberPolyData->GetNumberOfCells())); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); bool discard = false; for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); vnl_vector_fixed< double, 3 > v1; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; if (v1.magnitude()>0.001) { v1.normalize(); if (fabs(dot_product(v1,dir))>threshold) { discard = true; break; } } } if (!discard) { for (int j=0; jGetPoint(j, p1); vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); // UpdateColorCoding(); // UpdateFiberGeometry(); } bool mitk::FiberBundle::ApplyCurvatureThreshold(float minRadius, bool deleteFibers) { if (minRadius<0) return true; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Applying curvature threshold"; boost::progress_display disp(static_cast(m_FiberPolyData->GetNumberOfCells())); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); double p3[3]; points->GetPoint(j+2, p3); vnl_vector_fixed< float, 3 > v1, v2, v3; v1[0] = static_cast(p2[0]-p1[0]); v1[1] = static_cast(p2[1]-p1[1]); v1[2] = static_cast(p2[2]-p1[2]); v2[0] = static_cast(p3[0]-p2[0]); v2[1] = static_cast(p3[1]-p2[1]); v2[2] = static_cast(p3[2]-p2[2]); v3[0] = static_cast(p1[0]-p3[0]); v3[1] = static_cast(p1[1]-p3[1]); v3[2] = static_cast(p1[2]-p3[2]); float a = v1.magnitude(); float b = v2.magnitude(); float c = v3.magnitude(); float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle) vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); if (deleteFibers && rInsertNextCell(container); container = vtkSmartPointer::New(); } else if (j==numPoints-3) { id = vtkNewPoints->InsertNextPoint(p2); container->GetPointIds()->InsertNextId(id); id = vtkNewPoints->InsertNextPoint(p3); container->GetPointIds()->InsertNextId(id); vtkNewCells->InsertNextCell(container); } } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } bool mitk::FiberBundle::RemoveShortFibers(float lengthInMM) { MITK_INFO << "Removing short fibers"; if (lengthInMM<=0 || lengthInMMm_MaxFiberLength) // can't remove all fibers { MITK_WARN << "Process aborted. No fibers would be left!"; return false; } vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); float min = m_MaxFiberLength; boost::progress_display disp(m_NumFibers); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)>=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); if (m_FiberLengths.at(i)GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } bool mitk::FiberBundle::RemoveLongFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength) return true; if (lengthInMM vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Removing long fibers"; boost::progress_display disp(m_NumFibers); for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)<=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } void mitk::FiberBundle::ResampleSpline(float pointDistance, double tension, double continuity, double bias ) { if (pointDistance<=0) return; vtkSmartPointer vtkSmoothPoints = vtkSmartPointer::New(); //in smoothpoints the interpolated points representing a fiber are stored. //in vtkcells all polylines are stored, actually all id's of them are stored vtkSmartPointer vtkSmoothCells = vtkSmartPointer::New(); //cellcontainer for smoothed lines MITK_INFO << "Smoothing fibers"; vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); std::vector< vtkSmartPointer > resampled_streamlines; resampled_streamlines.resize(m_NumFibers); boost::progress_display disp(m_NumFibers); #pragma omp parallel for for (int i=0; i(m_NumFibers); i++) { vtkSmartPointer newPoints = vtkSmartPointer::New(); float length = 0; #pragma omp critical { length = m_FiberLengths.at(static_cast(i)); ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jInsertNextPoint(points->GetPoint(j)); } int sampling = static_cast(std::ceil(length/pointDistance)); vtkSmartPointer xSpline = vtkSmartPointer::New(); vtkSmartPointer ySpline = vtkSmartPointer::New(); vtkSmartPointer zSpline = vtkSmartPointer::New(); xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity); ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity); zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity); vtkSmartPointer spline = vtkSmartPointer::New(); spline->SetXSpline(xSpline); spline->SetYSpline(ySpline); spline->SetZSpline(zSpline); spline->SetPoints(newPoints); vtkSmartPointer functionSource = vtkSmartPointer::New(); functionSource->SetParametricFunction(spline); functionSource->SetUResolution(sampling); functionSource->SetVResolution(sampling); functionSource->SetWResolution(sampling); functionSource->Update(); vtkPolyData* outputFunction = functionSource->GetOutput(); vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber vtkSmartPointer smoothLine = vtkSmartPointer::New(); #pragma omp critical { for (int j=0; jGetNumberOfPoints(); j++) { vtkIdType id = vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j)); smoothLine->GetPointIds()->InsertNextId(id); } resampled_streamlines[static_cast(i)] = smoothLine; } } for (auto container : resampled_streamlines) { vtkSmoothCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkSmoothPoints); m_FiberPolyData->SetLines(vtkSmoothCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ResampleSpline(float pointDistance) { ResampleSpline(pointDistance, 0, 0, 0 ); } unsigned int mitk::FiberBundle::GetNumberOfPoints() const { unsigned int points = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); points += cell->GetNumberOfPoints(); } return points; } void mitk::FiberBundle::Compress(float error) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Compressing fibers"; unsigned int numRemovedPoints = 0; boost::progress_display disp(static_cast(m_FiberPolyData->GetNumberOfCells())); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); #pragma omp parallel for for (int i=0; i(m_FiberPolyData->GetNumberOfCells()); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; #pragma omp critical { ++disp; weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } // calculate curvatures auto numPoints = vertices.size(); std::vector< int > removedPoints; removedPoints.resize(numPoints, 0); removedPoints[0]=-1; removedPoints[numPoints-1]=-1; vtkSmartPointer container = vtkSmartPointer::New(); unsigned int remCounter = 0; bool pointFound = true; while (pointFound) { pointFound = false; double minError = static_cast(error); unsigned int removeIndex = 0; for (unsigned int j=0; j candV = vertices.at(j); int validP = -1; vnl_vector_fixed< double, 3 > pred; for (int k=static_cast(j)-1; k>=0; k--) if (removedPoints[static_cast(k)]<=0) { pred = vertices.at(static_cast(k)); validP = k; break; } int validS = -1; vnl_vector_fixed< double, 3 > succ; for (unsigned int k=j+1; k(k); break; } if (validP>=0 && validS>=0) { double a = (candV-pred).magnitude(); double b = (candV-succ).magnitude(); double c = (pred-succ).magnitude(); double s=0.5*(a+b+c); double hc=(2.0/c)*sqrt(fabs(s*(s-a)*(s-b)*(s-c))); if (hcInsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } #pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); numRemovedPoints += remCounter; vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()>0) { MITK_INFO << "Removed points: " << numRemovedPoints; SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } void mitk::FiberBundle::ResampleToNumPoints(unsigned int targetPoints) { if (targetPoints<2) mitkThrow() << "Minimum two points required for resampling!"; MITK_INFO << "Resampling fibers (number of points " << targetPoints << ")"; bool unequal_fibs = true; while (unequal_fibs) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); unequal_fibs = false; for (unsigned int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; double seg_len = 0; { weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); if (numPoints!=targetPoints) seg_len = static_cast(this->GetFiberLength(i)/(targetPoints-1)); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } vtkSmartPointer container = vtkSmartPointer::New(); vnl_vector_fixed< double, 3 > lastV = vertices.at(0); { vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block()); container->GetPointIds()->InsertNextId(id); } for (unsigned int j=1; j vec = vertices.at(j) - lastV; double new_dist = vec.magnitude(); if (new_dist >= seg_len && seg_len>0) { vnl_vector_fixed< double, 3 > newV = lastV; if ( new_dist-seg_len <= mitk::eps ) { vec.normalize(); newV += vec * seg_len; } else { // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p') vnl_vector_fixed< double, 3 > p = vertices.at(j-1); vnl_vector_fixed< double, 3 > d = vertices.at(j) - p; double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2])); double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - seg_len*seg_len; double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a); double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a); if (v1>0) newV = p + d * v1; else if (v2>0) newV = p + d * v2; else MITK_INFO << "ERROR1 - linear resampling"; j--; } //#pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block()); container->GetPointIds()->InsertNextId(id); } lastV = newV; } else if ( (j==vertices.size()-1 && new_dist>0.0001) || seg_len<=0.0000001) { //#pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } //#pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); vtkNewCells->InsertNextCell(container); if (container->GetNumberOfPoints()!=targetPoints) unequal_fibs = true; } } if (vtkNewCells->GetNumberOfCells()>0) { SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } } void mitk::FiberBundle::ResampleLinear(double pointDistance) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Resampling fibers (linear)"; boost::progress_display disp(static_cast(m_FiberPolyData->GetNumberOfCells())); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); std::vector< vtkSmartPointer > resampled_streamlines; resampled_streamlines.resize(static_cast(m_FiberPolyData->GetNumberOfCells())); #pragma omp parallel for for (int i=0; i(m_FiberPolyData->GetNumberOfCells()); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; #pragma omp critical { ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } vtkSmartPointer container = vtkSmartPointer::New(); vnl_vector_fixed< double, 3 > lastV = vertices.at(0); #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block()); container->GetPointIds()->InsertNextId(id); } for (unsigned int j=1; j vec = vertices.at(j) - lastV; double new_dist = vec.magnitude(); if (new_dist >= pointDistance) { vnl_vector_fixed< double, 3 > newV = lastV; if ( new_dist-pointDistance <= mitk::eps ) { vec.normalize(); newV += vec * pointDistance; } else { // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p') vnl_vector_fixed< double, 3 > p = vertices.at(j-1); vnl_vector_fixed< double, 3 > d = vertices.at(j) - p; double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2])); double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - pointDistance*pointDistance; double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a); double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a); if (v1>0) newV = p + d * v1; else if (v2>0) newV = p + d * v2; else MITK_INFO << "ERROR1 - linear resampling"; j--; } #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block()); container->GetPointIds()->InsertNextId(id); } lastV = newV; } else if (j==vertices.size()-1 && new_dist>0.0001) { #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } #pragma omp critical { resampled_streamlines[static_cast(i)] = container; } } for (auto container : resampled_streamlines) { vtkNewCells->InsertNextCell(container); } if (vtkNewCells->GetNumberOfCells()>0) { m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } // reapply selected colorcoding in case PolyData structure has changed bool mitk::FiberBundle::Equals(mitk::FiberBundle* fib, double eps) { if (fib==nullptr) { MITK_INFO << "Reference bundle is nullptr!"; return false; } if (m_NumFibers!=fib->GetNumFibers()) { MITK_INFO << "Unequal number of fibers!"; MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers(); return false; } for (unsigned int i=0; iGetCell(i); auto numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i); auto numPoints2 = cell2->GetNumberOfPoints(); vtkPoints* points2 = cell2->GetPoints(); if (numPoints2!=numPoints) { MITK_INFO << "Unequal number of points in fiber " << i << "!"; MITK_INFO << numPoints2 << " vs. " << numPoints; return false; } for (int j=0; jGetPoint(j); double* p2 = points2->GetPoint(j); if (fabs(p1[0]-p2[0])>eps || fabs(p1[1]-p2[1])>eps || fabs(p1[2]-p2[2])>eps) { MITK_INFO << "Unequal points in fiber " << i << " at position " << j << "!"; MITK_INFO << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2]; MITK_INFO << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2]; return false; } } } return true; } void mitk::FiberBundle::PrintSelf(std::ostream &os, itk::Indent indent) const { os << this->GetNameOfClass() << ":\n"; os << indent << "Number of fibers: " << this->GetNumFibers() << std::endl; os << indent << "Min. fiber length: " << this->GetMinFiberLength() << std::endl; os << indent << "Max. fiber length: " << this->GetMaxFiberLength() << std::endl; os << indent << "Mean fiber length: " << this->GetMeanFiberLength() << std::endl; os << indent << "Median fiber length: " << this->GetMedianFiberLength() << std::endl; os << indent << "STDEV fiber length: " << this->GetLengthStDev() << std::endl; os << indent << "Number of points: " << this->GetNumberOfPoints() << std::endl; os << indent << "Extent x: " << this->GetGeometry()->GetExtentInMM(0) << "mm" << std::endl; os << indent << "Extent y: " << this->GetGeometry()->GetExtentInMM(1) << "mm" << std::endl; os << indent << "Extent z: " << this->GetGeometry()->GetExtentInMM(2) << "mm" << std::endl; os << indent << "Diagonal: " << this->GetGeometry()->GetDiagonalLength() << "mm" << std::endl; os << "\nReference geometry:" << std::endl; - if (this->GetReferenceGeometry().IsNotNull()) - { - os << "Matrix:\n" << this->GetReferenceGeometry()->GetIndexToWorldTransform()->GetMatrix(); - os << "Origin: " << this->GetReferenceGeometry()->GetOrigin() << std::endl; - os << "Spacing: " << this->GetReferenceGeometry()->GetSpacing() << std::endl; - } - else - os << "NONE" << std::endl; + os << indent << "Size: [" << std::defaultfloat << m_TrackVisHeader.dim[0] << " " << m_TrackVisHeader.dim[1] << " " << m_TrackVisHeader.dim[2] << "]" << std::endl; + os << indent << "Voxel size: [" << m_TrackVisHeader.voxel_size[0] << " " << m_TrackVisHeader.voxel_size[1] << " " << m_TrackVisHeader.voxel_size[2] << "]" << std::endl; + os << indent << "Origin: [" << m_TrackVisHeader.origin[0] << " " << m_TrackVisHeader.origin[1] << " " << m_TrackVisHeader.origin[2] << "]" << std::endl; + os << indent << "Matrix: " << std::scientific << std::endl; + os << indent << "[[" << m_TrackVisHeader.vox_to_ras[0][0] << ", " << m_TrackVisHeader.vox_to_ras[0][1] << ", " << m_TrackVisHeader.vox_to_ras[0][2] << ", " << m_TrackVisHeader.vox_to_ras[0][3] << "]" << std::endl; + os << indent << " [" << m_TrackVisHeader.vox_to_ras[1][0] << ", " << m_TrackVisHeader.vox_to_ras[1][1] << ", " << m_TrackVisHeader.vox_to_ras[1][2] << ", " << m_TrackVisHeader.vox_to_ras[1][3] << "]" << std::endl; + os << indent << " [" << m_TrackVisHeader.vox_to_ras[2][0] << ", " << m_TrackVisHeader.vox_to_ras[2][1] << ", " << m_TrackVisHeader.vox_to_ras[2][2] << ", " << m_TrackVisHeader.vox_to_ras[2][3] << "]" << std::endl; + os << indent << " [" << m_TrackVisHeader.vox_to_ras[3][0] << ", " << m_TrackVisHeader.vox_to_ras[3][1] << ", " << m_TrackVisHeader.vox_to_ras[3][2] << ", " << m_TrackVisHeader.vox_to_ras[3][3] << "]]" << std::defaultfloat << std::endl; if (m_FiberWeights!=nullptr) { std::vector< float > weights; for (int i=0; iGetSize(); i++) weights.push_back(m_FiberWeights->GetValue(i)); std::sort(weights.begin(), weights.end()); os << "\nFiber weight statistics" << std::endl; os << indent << "Min: " << weights.front() << std::endl; os << indent << "1% quantile: " << weights.at(static_cast(weights.size()*0.01)) << std::endl; os << indent << "5% quantile: " << weights.at(static_cast(weights.size()*0.05)) << std::endl; os << indent << "25% quantile: " << weights.at(static_cast(weights.size()*0.25)) << std::endl; os << indent << "Median: " << weights.at(static_cast(weights.size()*0.5)) << std::endl; os << indent << "75% quantile: " << weights.at(static_cast(weights.size()*0.75)) << std::endl; os << indent << "95% quantile: " << weights.at(static_cast(weights.size()*0.95)) << std::endl; os << indent << "99% quantile: " << weights.at(static_cast(weights.size()*0.99)) << std::endl; os << indent << "Max: " << weights.back() << std::endl; } else os << indent << "\n\nNo fiber weight array found." << std::endl; Superclass::PrintSelf(os, indent); } -mitk::FiberBundle::TrackVis_header mitk::FiberBundle::GetTrackVisHeader() const +mitk::FiberBundle::TrackVis_header mitk::FiberBundle::GetTrackVisHeader() { + if (m_TrackVisHeader.hdr_size==0) + { + mitk::Geometry3D::Pointer geom = dynamic_cast(this->GetGeometry()); + SetTrackVisHeader(geom); + } return m_TrackVisHeader; } void mitk::FiberBundle::SetTrackVisHeader(const mitk::FiberBundle::TrackVis_header &TrackVisHeader) { m_TrackVisHeader = TrackVisHeader; } +void mitk::FiberBundle::SetTrackVisHeader(mitk::BaseGeometry* geometry) +{ + vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New(); + matrix->Identity(); + + if (geometry==nullptr) + return; + + for(int i=0; i<3 ;i++) + { + m_TrackVisHeader.dim[i] = geometry->GetExtent(i); + m_TrackVisHeader.voxel_size[i] = geometry->GetSpacing()[i]; + m_TrackVisHeader.origin[i] = geometry->GetOrigin()[i]; + matrix = geometry->GetVtkMatrix(); + } + + for (int i=0; i<4; ++i) + for (int j=0; j<4; ++j) + m_TrackVisHeader.vox_to_ras[i][j] = matrix->GetElement(i, j); + + m_TrackVisHeader.n_scalars = 0; + m_TrackVisHeader.n_properties = 0; + sprintf(m_TrackVisHeader.voxel_order,"LPS"); + m_TrackVisHeader.image_orientation_patient[0] = 1.0; + m_TrackVisHeader.image_orientation_patient[1] = 0.0; + m_TrackVisHeader.image_orientation_patient[2] = 0.0; + m_TrackVisHeader.image_orientation_patient[3] = 0.0; + m_TrackVisHeader.image_orientation_patient[4] = 1.0; + m_TrackVisHeader.image_orientation_patient[5] = 0.0; + m_TrackVisHeader.pad1[0] = 0; + m_TrackVisHeader.pad1[1] = 0; + m_TrackVisHeader.pad2[0] = 0; + m_TrackVisHeader.pad2[1] = 0; + m_TrackVisHeader.invert_x = 0; + m_TrackVisHeader.invert_y = 0; + m_TrackVisHeader.invert_z = 0; + m_TrackVisHeader.swap_xy = 0; + m_TrackVisHeader.swap_yz = 0; + m_TrackVisHeader.swap_zx = 0; + m_TrackVisHeader.n_count = 0; + m_TrackVisHeader.version = 2; + m_TrackVisHeader.hdr_size = 1000; + std::string id = "TRACK"; + strcpy(m_TrackVisHeader.id_string, id.c_str()); +} + /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundle::UpdateOutputInformation() { } void mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::FiberBundle::VerifyRequestedRegion() { return true; } void mitk::FiberBundle::SetRequestedRegion(const itk::DataObject* ) { } diff --git a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.h b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.h index 94ded14..8777861 100644 --- a/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.h +++ b/Modules/DiffusionCore/IODataStructures/mitkFiberBundle.h @@ -1,247 +1,244 @@ /*=================================================================== 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 _MITK_FiberBundle_H #define _MITK_FiberBundle_H //includes for MITK datastructure #include #include #include #include #include #include #include #include //includes storing fiberdata #include #include #include #include #include #include #include #include namespace mitk { /** * \brief Base Class for Fiber Bundles; */ class MITKDIFFUSIONCORE_EXPORT FiberBundle : public BaseData { public: typedef itk::Image ItkUcharImgType; // fiber colorcodings static const char* FIBER_ID_ARRAY; void UpdateOutputInformation() override; void SetRequestedRegionToLargestPossibleRegion() override; bool RequestedRegionIsOutsideOfTheBufferedRegion() override; bool VerifyRequestedRegion() override; void SetRequestedRegion(const itk::DataObject*) override; mitkClassMacro( FiberBundle, BaseData ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) mitkNewMacro1Param(Self, vtkSmartPointer) // custom constructor // colorcoding related methods void ColorFibersByFiberWeights(bool opacity, bool normalize); void ColorFibersByCurvature(bool opacity, bool normalize); void ColorFibersByLength(bool opacity, bool normalize); void ColorFibersByScalarMap(mitk::Image::Pointer, bool opacity, bool normalize); template void ColorFibersByScalarMap(const mitk::PixelType pixelType, mitk::Image::Pointer, bool opacity, bool normalize); void ColorFibersByOrientation(); void SetFiberOpacity(vtkDoubleArray *FAValArray); void ResetFiberOpacity(); void SetFiberColors(vtkSmartPointer fiberColors); void SetFiberColors(float r, float g, float b, float alpha=255); vtkSmartPointer GetFiberColors() const { return m_FiberColors; } // fiber compression void Compress(float error = 0.0); // fiber resampling void ResampleSpline(float pointDistance=1); void ResampleSpline(float pointDistance, double tension, double continuity, double bias ); void ResampleLinear(double pointDistance=1); void ResampleToNumPoints(unsigned int targetPoints); mitk::FiberBundle::Pointer FilterByWeights(float weight_thr, bool invert=false); bool RemoveShortFibers(float lengthInMM); bool RemoveLongFibers(float lengthInMM); bool ApplyCurvatureThreshold(float minRadius, bool deleteFibers); void MirrorFibers(unsigned int axis); void RotateAroundAxis(double x, double y, double z); void TranslateFibers(double x, double y, double z); void ScaleFibers(double x, double y, double z, bool subtractCenter=true); void TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz); void TransformFibers(itk::ScalableAffineTransform< mitk::ScalarType >::Pointer transform); void RemoveDir(vnl_vector_fixed dir, double threshold); template< class TType=float > void TransformPoint(itk::Point& point, itk::Matrix< TType, 3, 3>& rot, TType& tx, TType& ty, TType& tz) { mitk::Point3D center = this->GetGeometry()->GetCenter(); point[0] -= center[0]; point[1] -= center[1]; point[2] -= center[2]; point = rot*point; point[0] += center[0]+tx; point[1] += center[1]+ty; point[2] += center[2]+tz; } template< class TType=float > void TransformPoint(itk::Point& point, TType rx, TType ry, TType rz, TType tx, TType ty, TType tz) { auto rot = mitk::imv::GetRotationMatrixItk(rx, ry, rz); mitk::Point3D center = this->GetGeometry()->GetCenter(); point[0] -= center[0]; point[1] -= center[1]; point[2] -= center[2]; point = rot*point; point[0] += center[0]+tx; point[1] += center[1]+ty; point[2] += center[2]+tz; } itk::Matrix< double, 3, 3 > TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz); // add/subtract fibers FiberBundle::Pointer AddBundle(FiberBundle* fib); mitk::FiberBundle::Pointer AddBundles(std::vector< mitk::FiberBundle::Pointer > fibs); FiberBundle::Pointer SubtractBundle(FiberBundle* fib); // fiber subset extraction FiberBundle::Pointer ExtractFiberSubset(DataNode *roi, DataStorage* storage); std::vector ExtractFiberIdSubset(DataNode* roi, DataStorage* storage); FiberBundle::Pointer RemoveFibersOutside(ItkUcharImgType* mask, bool invert=false); float GetOverlap(ItkUcharImgType* mask); std::tuple GetDirectionalOverlap(ItkUcharImgType* mask, mitk::PeakImage::ItkPeakImageType* peak_image); float GetNumEpFractionInMask(ItkUcharImgType* mask, bool different_label); mitk::FiberBundle::Pointer SubsampleFibers(float factor, bool random_seed); // get/set data float GetFiberLength(unsigned int index) const { return m_FiberLengths.at(index); } vtkSmartPointer GetFiberWeights() const { return m_FiberWeights; } float GetFiberWeight(unsigned int fiber) const; void SetFiberWeights(float newWeight); void SetFiberWeight(unsigned int fiber, float weight); void SetFiberWeights(vtkSmartPointer weights); void SetFiberPolyData(vtkSmartPointer, bool updateGeometry = true); vtkSmartPointer GetFiberPolyData() const; itkGetConstMacro( NumFibers, unsigned int) //itkGetMacro( FiberSampling, int) itkGetConstMacro( MinFiberLength, float ) itkGetConstMacro( MaxFiberLength, float ) itkGetConstMacro( MeanFiberLength, float ) itkGetConstMacro( MedianFiberLength, float ) itkGetConstMacro( LengthStDev, float ) itkGetConstMacro( UpdateTime2D, itk::TimeStamp ) itkGetConstMacro( UpdateTime3D, itk::TimeStamp ) void RequestUpdate2D(){ m_UpdateTime2D.Modified(); } void RequestUpdate3D(){ m_UpdateTime3D.Modified(); } void RequestUpdate(){ m_UpdateTime2D.Modified(); m_UpdateTime3D.Modified(); } unsigned int GetNumberOfPoints() const; // copy fiber bundle mitk::FiberBundle::Pointer GetDeepCopy(); // compare fiber bundles bool Equals(FiberBundle* fib, double eps=0.01); - itkSetMacro( ReferenceGeometry, mitk::BaseGeometry::Pointer ) - itkGetConstMacro( ReferenceGeometry, mitk::BaseGeometry::Pointer ) - vtkSmartPointer GeneratePolyDataByIds(std::vector fiberIds, vtkSmartPointer weights); // Structure to hold metadata of a TrackVis file struct TrackVis_header { char id_string[6]; short int dim[3]; float voxel_size[3]; float origin[3]; short int n_scalars; char scalar_name[10][20]; short int n_properties; char property_name[10][20]; float vox_to_ras[4][4]; char reserved[444]; char voxel_order[4]; char pad2[4]; float image_orientation_patient[6]; char pad1[2]; unsigned char invert_x; unsigned char invert_y; unsigned char invert_z; unsigned char swap_xy; unsigned char swap_yz; unsigned char swap_zx; int n_count; int version; int hdr_size; }; - TrackVis_header GetTrackVisHeader() const; + TrackVis_header GetTrackVisHeader(); void SetTrackVisHeader(const TrackVis_header &TrackVisHeader); + void SetTrackVisHeader(BaseGeometry *geometry); protected: FiberBundle( vtkPolyData* fiberPolyData = nullptr ); ~FiberBundle() override; void GenerateFiberIds(); void UpdateFiberGeometry(); void PrintSelf(std::ostream &os, itk::Indent indent) const override; private: // actual fiber container vtkSmartPointer m_FiberPolyData; // contains fiber ids vtkSmartPointer m_FiberIdDataSet; unsigned int m_NumFibers; vtkSmartPointer m_FiberColors; vtkSmartPointer m_FiberWeights; std::vector< float > m_FiberLengths; float m_MinFiberLength; float m_MaxFiberLength; float m_MeanFiberLength; float m_MedianFiberLength; float m_LengthStDev; itk::TimeStamp m_UpdateTime2D; itk::TimeStamp m_UpdateTime3D; - mitk::BaseGeometry::Pointer m_ReferenceGeometry; TrackVis_header m_TrackVisHeader; }; } // namespace mitk #endif /* _MITK_FiberBundle_H */ diff --git a/Modules/DiffusionCore/IODataStructures/mitkTrackvis.cpp b/Modules/DiffusionCore/IODataStructures/mitkTrackvis.cpp index 92e5d2c..62bb244 100644 --- a/Modules/DiffusionCore/IODataStructures/mitkTrackvis.cpp +++ b/Modules/DiffusionCore/IODataStructures/mitkTrackvis.cpp @@ -1,304 +1,282 @@ #include #include TrackVisFiberReader::TrackVisFiberReader() { m_Filename = ""; m_FilePointer = nullptr; } TrackVisFiberReader::~TrackVisFiberReader() { if (m_FilePointer) fclose( m_FilePointer ); } - // Create a TrackVis file and store standard metadata. The file is ready to append fibers. // --------------------------------------------------------------------------------------- -short TrackVisFiberReader::create(std::string filename , const mitk::FiberBundle *fib, bool print_header) +short TrackVisFiberReader::create(std::string filename , mitk::FiberBundle *fib, bool print_header) { - if (fib->GetTrackVisHeader().hdr_size > 0) - { - m_Header = fib->GetTrackVisHeader(); - } - else - { - - vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New(); - matrix->Identity(); - - for(int i=0; i<3 ;i++) - { - if (fib->GetReferenceGeometry().IsNotNull()) - { - m_Header.dim[i] = fib->GetReferenceGeometry()->GetExtent(i); - m_Header.voxel_size[i] = fib->GetReferenceGeometry()->GetSpacing()[i]; - m_Header.origin[i] = fib->GetReferenceGeometry()->GetOrigin()[i]; - matrix = fib->GetReferenceGeometry()->GetVtkMatrix(); - } - else - { - m_Header.dim[i] = fib->GetGeometry()->GetExtent(i); - m_Header.voxel_size[i] = 1; - m_Header.origin[i] = 0; - } - } - - - for (int i=0; i<4; ++i) - for (int j=0; j<4; ++j) - m_Header.vox_to_ras[i][j] = matrix->GetElement(i, j); - - m_Header.n_scalars = 0; - m_Header.n_properties = 0; - sprintf(m_Header.voxel_order,"LPS"); - m_Header.image_orientation_patient[0] = 1.0; - m_Header.image_orientation_patient[1] = 0.0; - m_Header.image_orientation_patient[2] = 0.0; - m_Header.image_orientation_patient[3] = 0.0; - m_Header.image_orientation_patient[4] = 1.0; - m_Header.image_orientation_patient[5] = 0.0; - m_Header.pad1[0] = 0; - m_Header.pad1[1] = 0; - m_Header.pad2[0] = 0; - m_Header.pad2[1] = 0; - m_Header.invert_x = 0; - m_Header.invert_y = 0; - m_Header.invert_z = 0; - m_Header.swap_xy = 0; - m_Header.swap_yz = 0; - m_Header.swap_zx = 0; - m_Header.n_count = 0; - m_Header.version = 2; - m_Header.hdr_size = 1000; - std::string id = "TRACK"; - strcpy(m_Header.id_string, id.c_str()); - } + m_Header = fib->GetTrackVisHeader(); if (print_header) this->print_header(); // write the header to the file m_FilePointer = fopen(filename.c_str(),"w+b"); if (m_FilePointer == nullptr) { printf("[ERROR] Unable to create file '%s'\n",filename.c_str()); return 0; } if (fwrite((char*)&m_Header, 1, 1000, m_FilePointer) != 1000) MITK_ERROR << "TrackVis::create : Error occurding during writing fiber."; this->m_Filename = filename; return 1; } // Open an existing TrackVis file and read metadata information. // The file pointer is positiond at the beginning of fibers data // ------------------------------------------------------------- short TrackVisFiberReader::open( std::string filename ) { m_FilePointer = std::fopen(filename.c_str(), "rb"); if (m_FilePointer == nullptr) { printf("[ERROR] Unable to open file '%s'\n",filename.c_str()); return 0; } this->m_Filename = filename; return fread((char*)(&m_Header), 1, 1000, m_FilePointer); } short TrackVisFiberReader::write(const mitk::FiberBundle *fib) { vtkSmartPointer poly = fib->GetFiberPolyData(); - if (fib->GetReferenceGeometry().IsNotNull()) { + mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); + vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New(); + matrix->Identity(); + for (int i=0; i<4; ++i) + for (int j=0; j<4; ++j) + { + if (j<3) + matrix->SetElement(i, j, m_Header.vox_to_ras[i][j]/m_Header.voxel_size[j]); + else + matrix->SetElement(i, j, m_Header.vox_to_ras[i][j]); + } + + if (m_Header.voxel_order[0]=='R') + { + matrix->SetElement(0,0,-matrix->GetElement(0,0)); + matrix->SetElement(0,1,-matrix->GetElement(0,1)); + matrix->SetElement(0,2,-matrix->GetElement(0,2)); + matrix->SetElement(0,3,-matrix->GetElement(0,3)); + } + if (m_Header.voxel_order[1]=='A') + { + matrix->SetElement(1,0,-matrix->GetElement(1,0)); + matrix->SetElement(1,1,-matrix->GetElement(1,1)); + matrix->SetElement(1,2,-matrix->GetElement(1,2)); + matrix->SetElement(1,3,-matrix->GetElement(1,3)); + } + if (m_Header.voxel_order[2]=='I') + { + matrix->SetElement(2,0,-matrix->GetElement(2,0)); + matrix->SetElement(2,1,-matrix->GetElement(2,1)); + matrix->SetElement(2,2,-matrix->GetElement(2,2)); + matrix->SetElement(2,3,-matrix->GetElement(2,3)); + } + + geometry->SetIndexToWorldTransformByVtkMatrix(matrix); + vtkSmartPointer transformFilter = vtkSmartPointer::New(); transformFilter->SetInputData(poly); - transformFilter->SetTransform(fib->GetReferenceGeometry()->GetVtkTransform()->GetInverse()); + transformFilter->SetTransform(geometry->GetVtkTransform()->GetInverse()); transformFilter->Update(); poly = transformFilter->GetOutput(); } for (unsigned int i=0; iGetNumFibers(); i++) { vtkCell* cell = poly->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); unsigned int numSaved, pos = 0; //float* tmp = new float[3*maxSteps]; std::vector< float > tmp; tmp.reserve(3*numPoints); numSaved = numPoints; for(unsigned int i=0; iGetPoint(i); tmp[pos++] = p[0]; tmp[pos++] = p[1]; tmp[pos++] = p[2]; } // write the coordinates to the file if ( fwrite((char*)&numSaved, 1, 4, m_FilePointer) != 4 ) { printf( "[ERROR] Problems saving the fiber!\n" ); return 1; } if ( fwrite((char*)&(tmp.front()), 1, 4*pos, m_FilePointer) != 4*pos ) { printf( "[ERROR] Problems saving the fiber!\n" ); return 1; } } return 0; } void TrackVisFiberReader::print_header() { std::cout << "--------------------------------------------------------" << std::endl; std::cout << "see http://trackvis.org/docs/?subsect=fileformat" << std::endl; std::cout << "ONLY vox_to_ras AND voxel_order HEADER ENTRIES ARE USED FOR FIBER COORDINATE TRANSFORMATIONS!" << std::endl; std::cout << "\nid_string (should be \"TRACK\"): " << m_Header.id_string << std::endl; std::cout << "dim: [" << std::defaultfloat << m_Header.dim[0] << " " << m_Header.dim[1] << " " << m_Header.dim[2] << "]" << std::endl; std::cout << "voxel_size: [" << m_Header.voxel_size[0] << " " << m_Header.voxel_size[1] << " " << m_Header.voxel_size[2] << "]" << std::endl; std::cout << "origin: [" << m_Header.origin[0] << " " << m_Header.origin[1] << " " << m_Header.origin[2] << "]" << std::endl; std::cout << "vox_to_world: " << std::scientific << std::endl; std::cout << "[[" << m_Header.vox_to_ras[0][0] << ", " << m_Header.vox_to_ras[0][1] << ", " << m_Header.vox_to_ras[0][2] << ", " << m_Header.vox_to_ras[0][3] << "]" << std::endl; std::cout << " [" << m_Header.vox_to_ras[1][0] << ", " << m_Header.vox_to_ras[1][1] << ", " << m_Header.vox_to_ras[1][2] << ", " << m_Header.vox_to_ras[1][3] << "]" << std::endl; std::cout << " [" << m_Header.vox_to_ras[2][0] << ", " << m_Header.vox_to_ras[2][1] << ", " << m_Header.vox_to_ras[2][2] << ", " << m_Header.vox_to_ras[2][3] << "]" << std::endl; std::cout << " [" << m_Header.vox_to_ras[3][0] << ", " << m_Header.vox_to_ras[3][1] << ", " << m_Header.vox_to_ras[3][2] << ", " << m_Header.vox_to_ras[3][3] << "]]" << std::defaultfloat << std::endl; std::cout << "voxel_order: " << m_Header.voxel_order[0] << m_Header.voxel_order[1] << m_Header.voxel_order[2] << std::endl; std::cout << "pad1: " << m_Header.pad1[0] << m_Header.pad1[1] << std::endl; std::cout << "pad2: " << m_Header.pad2[0] << m_Header.pad2[1] << m_Header.pad2[2] << std::endl; std::cout << "image_orientation_patient: [" << m_Header.image_orientation_patient[0] << " " << m_Header.image_orientation_patient[1] << " " << m_Header.image_orientation_patient[2] << " " << m_Header.image_orientation_patient[3] << " " << m_Header.image_orientation_patient[4] << " " << m_Header.image_orientation_patient[5] << "]" << std::endl; std::cout << "invert_x: " << static_cast(m_Header.invert_x) << std::endl; std::cout << "invert_y: " << static_cast(m_Header.invert_y) << std::endl; std::cout << "invert_z: " << static_cast(m_Header.invert_z) << std::endl; std::cout << "swap_xy: " << static_cast(m_Header.swap_xy) << std::endl; std::cout << "swap_yz: " << static_cast(m_Header.swap_yz) << std::endl; std::cout << "swap_zx: " << static_cast(m_Header.swap_zx) << std::endl; std::cout << "n_count: " << m_Header.n_count << std::endl; std::cout << "version: " << m_Header.version << std::endl; std::cout << "hdr_size: " << m_Header.hdr_size << std::endl; std::cout << "\nNot printed: n_scalars, scalar_name, n_properties, property_name, reserved" << std::endl; std::cout << "--------------------------------------------------------" << std::endl; } short TrackVisFiberReader::read( mitk::FiberBundle* fib, bool use_matrix, bool print_header) { int numPoints; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); if (print_header) this->print_header(); while (fread((char*)&numPoints, 1, 4, m_FilePointer)==4) { if ( numPoints <= 0 ) { printf( "[ERROR] Trying to read a fiber with %d points!\n", numPoints ); return -1; } vtkSmartPointer container = vtkSmartPointer::New(); float tmp[3]; for(int i=0; iInsertNextPoint(tmp); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); fiberPolyData->SetPoints(vtkNewPoints); fiberPolyData->SetLines(vtkNewCells); mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New(); matrix->Identity(); if (use_matrix) for (int i=0; i<4; ++i) for (int j=0; j<4; ++j) { if (j<3) matrix->SetElement(i, j, m_Header.vox_to_ras[i][j]/m_Header.voxel_size[j]); else matrix->SetElement(i, j, m_Header.vox_to_ras[i][j]); } if (m_Header.voxel_order[0]=='R') { matrix->SetElement(0,0,-matrix->GetElement(0,0)); matrix->SetElement(0,1,-matrix->GetElement(0,1)); matrix->SetElement(0,2,-matrix->GetElement(0,2)); matrix->SetElement(0,3,-matrix->GetElement(0,3)); } if (m_Header.voxel_order[1]=='A') { matrix->SetElement(1,0,-matrix->GetElement(1,0)); matrix->SetElement(1,1,-matrix->GetElement(1,1)); matrix->SetElement(1,2,-matrix->GetElement(1,2)); matrix->SetElement(1,3,-matrix->GetElement(1,3)); } if (m_Header.voxel_order[2]=='I') { matrix->SetElement(2,0,-matrix->GetElement(2,0)); matrix->SetElement(2,1,-matrix->GetElement(2,1)); matrix->SetElement(2,2,-matrix->GetElement(2,2)); matrix->SetElement(2,3,-matrix->GetElement(2,3)); } geometry->SetIndexToWorldTransformByVtkMatrix(matrix); vtkSmartPointer transformFilter = vtkSmartPointer::New(); transformFilter->SetInputData(fiberPolyData); transformFilter->SetTransform(geometry->GetVtkTransform()); transformFilter->Update(); fib->SetFiberPolyData(transformFilter->GetOutput()); - fib->SetReferenceGeometry(dynamic_cast(geometry.GetPointer())); + fib->SetTrackVisHeader(geometry.GetPointer()); fib->SetTrackVisHeader(m_Header); return numPoints; } // Update the field in the header to the new FIBER TOTAL. // ------------------------------------------------------ void TrackVisFiberReader::updateTotal( int totFibers ) { fseek(m_FilePointer, 1000-12, SEEK_SET); if (fwrite((char*)&totFibers, 1, 4, m_FilePointer) != 4) MITK_ERROR << "[ERROR] Problems saving the fiber!"; } void TrackVisFiberReader::writeHdr() { fseek(m_FilePointer, 0, SEEK_SET); if (fwrite((char*)&m_Header, 1, 1000, m_FilePointer) != 1000) MITK_ERROR << "[ERROR] Problems saving the fiber!"; } // Close the TrackVis file, but keep the metadata in the header. // ------------------------------------------------------------- void TrackVisFiberReader::close() { fclose(m_FilePointer); m_FilePointer = nullptr; } bool TrackVisFiberReader::IsTransformValid() { if (fabs(m_Header.image_orientation_patient[0])<=0.001 || fabs(m_Header.image_orientation_patient[3])<=0.001 || fabs(m_Header.image_orientation_patient[5])<=0.001) return false; return true; } diff --git a/Modules/DiffusionCore/IODataStructures/mitkTrackvis.h b/Modules/DiffusionCore/IODataStructures/mitkTrackvis.h index f45d628..667f0af 100644 --- a/Modules/DiffusionCore/IODataStructures/mitkTrackvis.h +++ b/Modules/DiffusionCore/IODataStructures/mitkTrackvis.h @@ -1,53 +1,53 @@ /*=================================================================== 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 _TRACKVIS #define _TRACKVIS #include #include #include #include #include #include #include #include // Class to handle TrackVis files. // ------------------------------- class MITKDIFFUSIONCORE_EXPORT TrackVisFiberReader { private: std::string m_Filename; FILE* m_FilePointer; mitk::FiberBundle::TrackVis_header m_Header; public: - short create(std::string m_Filename, const mitk::FiberBundle* fib, bool print_header); + short create(std::string m_Filename, mitk::FiberBundle *fib, bool print_header); short open(std::string m_Filename ); short read(mitk::FiberBundle* fib , bool use_matrix, bool print_header); short write(const mitk::FiberBundle* fib ); void writeHdr(); void updateTotal( int totFibers ); void close(); bool IsTransformValid(); void print_header(); TrackVisFiberReader(); ~TrackVisFiberReader(); }; #endif diff --git a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisWriter.cpp b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisWriter.cpp index 7978591..8e2c5b6 100644 --- a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisWriter.cpp +++ b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisWriter.cpp @@ -1,108 +1,108 @@ /*=================================================================== 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 "mitkFiberBundleTrackVisWriter.h" #include #include #include #include #include #include #include #include #include #include #include "mitkDiffusionIOMimeTypes.h" mitk::FiberBundleTrackVisWriter::FiberBundleTrackVisWriter() : mitk::AbstractFileWriter(mitk::FiberBundle::GetStaticNameOfClass(), mitk::DiffusionIOMimeTypes::FIBERBUNDLE_TRK_MIMETYPE_NAME(), "TrackVis Fiber Bundle Reader") { Options defaultOptions; defaultOptions["Print header"] = false; this->SetDefaultOptions(defaultOptions); RegisterService(); } mitk::FiberBundleTrackVisWriter::FiberBundleTrackVisWriter(const mitk::FiberBundleTrackVisWriter & other) :mitk::AbstractFileWriter(other) {} mitk::FiberBundleTrackVisWriter::~FiberBundleTrackVisWriter() {} mitk::FiberBundleTrackVisWriter * mitk::FiberBundleTrackVisWriter::Clone() const { return new mitk::FiberBundleTrackVisWriter(*this); } void mitk::FiberBundleTrackVisWriter::Write() { std::ostream* out; std::ofstream outStream; if( this->GetOutputStream() ) { out = this->GetOutputStream(); }else{ outStream.open( this->GetOutputLocation().c_str() ); out = &outStream; } if ( !out->good() ) { mitkThrow() << "Stream not good."; } try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); setlocale(LC_ALL, locale.c_str()); std::locale I("C"); out->imbue(I); std::string filename = this->GetOutputLocation().c_str(); mitk::FiberBundle::ConstPointer input = dynamic_cast(this->GetInput()); std::string ext = itksys::SystemTools::GetFilenameLastExtension(this->GetOutputLocation().c_str()); // default extension is .fib if(ext == "") { ext = ".trk"; this->SetOutputLocation(this->GetOutputLocation() + ext); } Options options = this->GetOptions(); bool print_header = us::any_cast(options["Print header"]); MITK_INFO << "Writing fiber bundle as TRK"; TrackVisFiberReader trk; - trk.create(filename, input.GetPointer(), print_header); + trk.create(filename, const_cast(input.GetPointer()), print_header); trk.writeHdr(); trk.write(input.GetPointer()); setlocale(LC_ALL, currLocale.c_str()); MITK_INFO << "TrackVis Fiber bundle written to " << filename; } catch(...) { throw; } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp index 7ad7ab5..9d41f0c 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkFiberProcessingView.cpp @@ -1,1731 +1,1730 @@ /*=================================================================== 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 "QmitkFiberProcessingView.h" #include #include #include #include #include #include #include #include #include #include #include #include "usModuleRegistry.h" #include #include "mitkNodePredicateDataType.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include const std::string QmitkFiberProcessingView::VIEW_ID = "org.mitk.views.fiberprocessing"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace mitk; QmitkFiberProcessingView::QmitkFiberProcessingView() : QmitkAbstractView() , m_Controls( 0 ) , m_CircleCounter(0) , m_PolygonCounter(0) , m_UpsamplingFactor(1) { } // Destructor QmitkFiberProcessingView::~QmitkFiberProcessingView() { RemoveObservers(); } void QmitkFiberProcessingView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkFiberProcessingViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_CircleButton, SIGNAL( clicked() ), this, SLOT( OnDrawCircle() ) ); connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ), this, SLOT( OnDrawPolygon() ) ); connect(m_Controls->PFCompoANDButton, SIGNAL(clicked()), this, SLOT(GenerateAndComposite()) ); connect(m_Controls->PFCompoORButton, SIGNAL(clicked()), this, SLOT(GenerateOrComposite()) ); connect(m_Controls->PFCompoNOTButton, SIGNAL(clicked()), this, SLOT(GenerateNotComposite()) ); connect(m_Controls->m_GenerateRoiImage, SIGNAL(clicked()), this, SLOT(GenerateRoiImage()) ); connect(m_Controls->m_JoinBundles, SIGNAL(clicked()), this, SLOT(JoinBundles()) ); connect(m_Controls->m_SubstractBundles, SIGNAL(clicked()), this, SLOT(SubstractBundles()) ); connect(m_Controls->m_CopyBundle, SIGNAL(clicked()), this, SLOT(CopyBundles()) ); connect(m_Controls->m_ExtractFibersButton, SIGNAL(clicked()), this, SLOT(Extract())); connect(m_Controls->m_RemoveButton, SIGNAL(clicked()), this, SLOT(Remove())); connect(m_Controls->m_ModifyButton, SIGNAL(clicked()), this, SLOT(Modify())); connect(m_Controls->m_ExtractionMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_RemovalMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ModificationMethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); connect(m_Controls->m_ExtractionBoxMask, SIGNAL(currentIndexChanged(int)), this, SLOT(OnMaskExtractionChanged())); m_Controls->m_ColorMapBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateDataType::Pointer isDwi = mitk::NodePredicateDataType::New("DiffusionImage"); mitk::NodePredicateDataType::Pointer isDti = mitk::NodePredicateDataType::New("TensorImage"); mitk::NodePredicateDataType::Pointer isOdf = mitk::NodePredicateDataType::New("OdfImage"); mitk::NodePredicateOr::Pointer isDiffusionImage = mitk::NodePredicateOr::New(isDwi, isDti); isDiffusionImage = mitk::NodePredicateOr::New(isDiffusionImage, isOdf); mitk::NodePredicateNot::Pointer noDiffusionImage = mitk::NodePredicateNot::New(isDiffusionImage); mitk::NodePredicateAnd::Pointer finalPredicate = mitk::NodePredicateAnd::New(isMitkImage, noDiffusionImage); m_Controls->m_ColorMapBox->SetPredicate(finalPredicate); } UpdateGui(); OnMaskExtractionChanged(); } void QmitkFiberProcessingView::OnMaskExtractionChanged() { m_Controls->m_FiberExtractionFractionLabel->setVisible(false); m_Controls->m_FiberExtractionFractionBox->setVisible(false); m_Controls->m_FiberExtractionThresholdLabel->setVisible(false); m_Controls->m_FiberExtractionThresholdBox->setVisible(false); m_Controls->m_InterpolateRoiBox->setVisible(false); m_Controls->m_BothEnds->setVisible(false); m_Controls->m_LabelsBox->setVisible(false); m_Controls->m_LabelsLabel->setVisible(false); if (m_Controls->m_ExtractionBoxMask->currentIndex() == 2 || m_Controls->m_ExtractionBoxMask->currentIndex() == 3) { m_Controls->m_FiberExtractionFractionLabel->setVisible(true); m_Controls->m_FiberExtractionFractionBox->setVisible(true); m_Controls->m_FiberExtractionThresholdLabel->setVisible(true); m_Controls->m_FiberExtractionThresholdBox->setVisible(true); m_Controls->m_InterpolateRoiBox->setVisible(true); } else if (m_Controls->m_ExtractionBoxMask->currentIndex() == 0 || m_Controls->m_ExtractionBoxMask->currentIndex() == 1) { if (m_Controls->m_ExtractionBoxMask->currentIndex() != 3) m_Controls->m_BothEnds->setVisible(true); m_Controls->m_InterpolateRoiBox->setVisible(true); m_Controls->m_FiberExtractionThresholdLabel->setVisible(true); m_Controls->m_FiberExtractionThresholdBox->setVisible(true); } else if (m_Controls->m_ExtractionBoxMask->currentIndex() == 4) { m_Controls->m_BothEnds->setVisible(true); m_Controls->m_LabelsBox->setVisible(true); m_Controls->m_LabelsLabel->setVisible(true); } } void QmitkFiberProcessingView::SetFocus() { m_Controls->toolBoxx->setFocus(); } void QmitkFiberProcessingView::Modify() { switch (m_Controls->m_ModificationMethodBox->currentIndex()) { case 0: { ResampleSelectedBundlesSpline(); break; } case 1: { ResampleSelectedBundlesLinear(); break; } case 2: { CompressSelectedBundles(); break; } case 3: { DoImageColorCoding(); break; } case 4: { MirrorFibers(); break; } case 5: { WeightFibers(); break; } case 6: { DoCurvatureColorCoding(); break; } case 7: { DoWeightColorCoding(); break; } case 8: { DoLengthColorCoding(); break; } } } void QmitkFiberProcessingView::WeightFibers() { float weight = this->m_Controls->m_BundleWeightBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->SetFiberWeights(weight); } } void QmitkFiberProcessingView::Remove() { switch (m_Controls->m_RemovalMethodBox->currentIndex()) { case 0: { RemoveDir(); break; } case 1: { PruneBundle(); break; } case 2: { ApplyCurvatureThreshold(); break; } case 3: { RemoveWithMask(false); break; } case 4: { RemoveWithMask(true); break; } case 5: { ApplyWeightThreshold(); break; } case 6: { ApplyDensityThreshold(); break; } } } void QmitkFiberProcessingView::Extract() { switch (m_Controls->m_ExtractionMethodBox->currentIndex()) { case 0: { ExtractWithPlanarFigure(); break; } case 1: { switch (m_Controls->m_ExtractionBoxMask->currentIndex()) { { case 0: ExtractWithMask(true, false, false); break; } { case 1: ExtractWithMask(true, true, false); break; } { case 2: ExtractWithMask(false, false, false); break; } { case 3: ExtractWithMask(false, true, false); break; } { case 4: ExtractWithMask(true, false, true); break; } } break; } } } void QmitkFiberProcessingView::PruneBundle() { int minLength = this->m_Controls->m_PruneFibersMinBox->value(); int maxLength = this->m_Controls->m_PruneFibersMaxBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); if (!fib->RemoveShortFibers(minLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); else if (!fib->RemoveLongFibers(maxLength)) QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyWeightThreshold() { float thr = this->m_Controls->m_WeightThresholdBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->FilterByWeights(thr); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByFiberWeights(false, true); node->SetData(newFib); } else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyDensityThreshold() { float thr = this->m_Controls->m_DensityThresholdBox->value(); float ol = this->m_Controls->m_DensityOverlapBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); itk::TractDensityImageFilter< ItkFloatImageType >::Pointer generator = itk::TractDensityImageFilter< ItkFloatImageType >::New(); generator->SetFiberBundle(fib); generator->SetBinaryOutput(false); generator->SetOutputAbsoluteValues(false); generator->Update(); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetRoiImages({generator->GetOutput()}); extractor->SetInputFiberBundle(fib); extractor->SetOverlapFraction(ol); extractor->SetInterpolate(true); extractor->SetThreshold(thr); extractor->SetNoNegatives(true); extractor->Update(); if (extractor->GetPositives().empty()) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } mitk::FiberBundle::Pointer newFib = extractor->GetPositives().at(0); if (newFib->GetNumFibers()>0) node->SetData(newFib); else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ApplyCurvatureThreshold() { int angle = this->m_Controls->m_CurvSpinBox->value(); int dist = this->m_Controls->m_CurvDistanceSpinBox->value(); std::vector< DataNode::Pointer > nodes = m_SelectedFB; for (auto node : nodes) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); itk::FiberCurvatureFilter::Pointer filter = itk::FiberCurvatureFilter::New(); filter->SetInputFiberBundle(fib); filter->SetAngularDeviation(angle); filter->SetDistance(dist); filter->SetRemoveFibers(m_Controls->m_RemoveCurvedFibersBox->isChecked()); filter->Update(); mitk::FiberBundle::Pointer newFib = filter->GetOutputFiberBundle(); if (newFib->GetNumFibers()>0) { newFib->ColorFibersByOrientation(); node->SetData(newFib); } else QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveDir() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); vnl_vector_fixed dir; dir[0] = m_Controls->m_ExtractDirX->value(); dir[1] = m_Controls->m_ExtractDirY->value(); dir[2] = m_Controls->m_ExtractDirZ->value(); fib->RemoveDir(dir,cos((float)m_Controls->m_ExtractAngle->value()*itk::Math::pi/180)); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::RemoveWithMask(bool removeInside) { if (m_RoiImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_RoiImageNode->GetData()); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(mitkMask, mask); mitk::FiberBundle::Pointer newFib = fib->RemoveFibersOutside(mask, removeInside); if (newFib->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } node->SetData(newFib); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ExtractWithMask(bool onlyEnds, bool invert, bool labelmap) { if (m_RoiImageNode.IsNull()) return; mitk::Image::Pointer mitkMask = dynamic_cast(m_RoiImageNode->GetData()); for (auto node : m_SelectedFB) { std::string roi_name = m_RoiImageNode->GetName(); mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(mitkMask, mask); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetInputFiberBundle(fib); extractor->SetRoiImages({mask}); extractor->SetRoiImageNames({roi_name}); extractor->SetThreshold(m_Controls->m_FiberExtractionThresholdBox->value()); extractor->SetOverlapFraction(m_Controls->m_FiberExtractionFractionBox->value()); extractor->SetBothEnds(m_Controls->m_BothEnds->isChecked()); extractor->SetInterpolate(m_Controls->m_InterpolateRoiBox->isChecked()); extractor->SetMinFibersPerTract(m_Controls->m_MinExtractedFibersBox->value()); extractor->SetSplitByRoi(true); if (invert) extractor->SetNoPositives(true); else extractor->SetNoNegatives(true); if (labelmap) { std::string labels_string = m_Controls->m_LabelsBox->text().toStdString(); if (labels_string!="ALL") { std::vector strs; boost::split(strs,labels_string,boost::is_any_of(" ,;\t")); std::vector< unsigned short > labels_vector; for (auto v : strs) { try{ unsigned short l = boost::lexical_cast(v); labels_vector.push_back(l); } catch(...) { } } extractor->SetLabels(labels_vector); } extractor->SetInterpolate(false); extractor->SetInputType(itk::FiberExtractionFilter::INPUT::LABEL_MAP); extractor->SetSplitLabels(true); onlyEnds = true; } if (onlyEnds) extractor->SetMode(itk::FiberExtractionFilter::MODE::ENDPOINTS); extractor->Update(); std::vector< mitk::FiberBundle::Pointer > newFibs; if (invert) newFibs = extractor->GetNegatives(); else newFibs = extractor->GetPositives(); if (newFibs.empty()) { QMessageBox::information(nullptr, "No output generated:", "No fibers could be extracted."); continue; } auto labels = extractor->GetPositiveLabels(); for (unsigned int i=0; iSetData(fib); std::string name = roi_name; if (iSetName(name); GetDataStorage()->Add(newNode, node); } node->SetVisibility(false); } } void QmitkFiberProcessingView::GenerateRoiImage() { if (m_SelectedPF.empty()) return; mitk::BaseGeometry::Pointer geometry; if (!m_SelectedFB.empty()) { mitk::FiberBundle::Pointer fib = dynamic_cast(m_SelectedFB.front()->GetData()); geometry = fib->GetGeometry(); } else if (m_SelectedImage) geometry = m_SelectedImage->GetGeometry(); else return; itk::Vector spacing = geometry->GetSpacing(); spacing /= m_UpsamplingFactor; mitk::Point3D newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); newOrigin[0] += bounds.GetElement(0); newOrigin[1] += bounds.GetElement(2); newOrigin[2] += bounds.GetElement(4); itk::Matrix direction; itk::ImageRegion<3> imageRegion; for (int i=0; i<3; i++) for (int j=0; j<3; j++) direction[j][i] = geometry->GetMatrixColumn(i)[j]/spacing[j]; imageRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); imageRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); imageRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); m_PlanarFigureImage = ItkUCharImageType::New(); m_PlanarFigureImage->SetSpacing( spacing ); // Set the image spacing m_PlanarFigureImage->SetOrigin( newOrigin ); // Set the image origin m_PlanarFigureImage->SetDirection( direction ); // Set the image direction m_PlanarFigureImage->SetRegions( imageRegion ); m_PlanarFigureImage->Allocate(); m_PlanarFigureImage->FillBuffer( 0 ); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); std::string name = m_SelectedPF.at(0)->GetName(); WritePfToImage(m_SelectedPF.at(0), tmpImage); for (unsigned int i=1; iGetName(); WritePfToImage(m_SelectedPF.at(i), tmpImage); } DataNode::Pointer node = DataNode::New(); tmpImage = Image::New(); tmpImage->InitializeByItk(m_PlanarFigureImage.GetPointer()); tmpImage->SetVolume(m_PlanarFigureImage->GetBufferPointer()); node->SetData(tmpImage); node->SetName(name); this->GetDataStorage()->Add(node); } void QmitkFiberProcessingView::WritePfToImage(mitk::DataNode::Pointer node, mitk::Image* image) { if (dynamic_cast(node->GetData())) { m_PlanarFigure = dynamic_cast(node->GetData()); AccessFixedDimensionByItk_2( image, InternalReorientImagePlane, 3, m_PlanarFigure->GetGeometry(), -1); AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 3, 2, node->GetName() ); } else if (dynamic_cast(node->GetData())) { DataStorage::SetOfObjects::ConstPointer children = GetDataStorage()->GetDerivations(node); for (unsigned int i=0; iSize(); i++) { WritePfToImage(children->at(i), image); } } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalReorientImagePlane( const itk::Image< TPixel, VImageDimension > *image, mitk::BaseGeometry* planegeo3D, int additionalIndex ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< float, VImageDimension > FloatImageType; typedef itk::ResampleImageFilter ResamplerType; typename ResamplerType::Pointer resampler = ResamplerType::New(); mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); float upsamp = m_UpsamplingFactor; float gausssigma = 0.5; // Spacing typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); spacing[0] = image->GetSpacing()[0] / upsamp; spacing[1] = image->GetSpacing()[1] / upsamp; spacing[2] = image->GetSpacing()[2]; resampler->SetOutputSpacing( spacing ); // Size typename ResamplerType::SizeType size; size[0] = planegeo->GetExtentInMM(0) / spacing[0]; size[1] = planegeo->GetExtentInMM(1) / spacing[1]; size[2] = 1; resampler->SetSize( size ); // Origin typename mitk::Point3D orig = planegeo->GetOrigin(); typename mitk::Point3D corrorig; planegeo3D->WorldToIndex(orig,corrorig); corrorig[0] += 0.5/upsamp; corrorig[1] += 0.5/upsamp; corrorig[2] += 0; planegeo3D->IndexToWorld(corrorig,corrorig); resampler->SetOutputOrigin(corrorig ); // Direction typename ResamplerType::DirectionType direction; typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); for(unsigned int c=0; cSetOutputDirection( direction ); // Gaussian interpolation if(gausssigma != 0) { double sigma[3]; for( unsigned int d = 0; d < 3; d++ ) sigma[d] = gausssigma * image->GetSpacing()[d]; double alpha = 2.0; typedef itk::GaussianInterpolateImageFunction GaussianInterpolatorType; typename GaussianInterpolatorType::Pointer interpolator = GaussianInterpolatorType::New(); interpolator->SetInputImage( image ); interpolator->SetParameters( sigma, alpha ); resampler->SetInterpolator( interpolator ); } else { typedef typename itk::LinearInterpolateImageFunction InterpolatorType; typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); interpolator->SetInputImage( image ); resampler->SetInterpolator( interpolator ); } resampler->SetInput( image ); resampler->SetDefaultPixelValue(0); resampler->Update(); if(additionalIndex < 0) { this->m_InternalImage = mitk::Image::New(); this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); } } template < typename TPixel, unsigned int VImageDimension > void QmitkFiberProcessingView::InternalCalculateMaskFromPlanarFigure( itk::Image< TPixel, VImageDimension > *image, unsigned int axis, std::string ) { typedef itk::Image< TPixel, VImageDimension > ImageType; // Generate mask image as new image with same header as input image and // initialize with "1". ItkUCharImageType::Pointer newMaskImage = ItkUCharImageType::New(); newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); newMaskImage->Allocate(); newMaskImage->FillBuffer( 1 ); // Generate VTK polygon from (closed) PlanarFigure polyline // (The polyline points are shifted by -0.5 in z-direction to make sure // that the extrusion filter, which afterwards elevates all points by +0.5 // in z-direction, creates a 3D object which is cut by the the plane z=0) const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const BaseGeometry *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); vtkPolyData *polyline = vtkPolyData::New(); polyline->Allocate( 1, 1 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // Create VTK polydata object of polyline contour vtkPoints *points = vtkPoints::New(); PlanarFigure::PolyLineType::const_iterator it; unsigned int numberOfPoints = 0; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected image Point2D point2D = *it; planarFigurePlaneGeometry->WorldToIndex(point2D, point2D); point2D[0] -= 0.5/m_UpsamplingFactor; point2D[1] -= 0.5/m_UpsamplingFactor; planarFigurePlaneGeometry->IndexToWorld(point2D, point2D); planarFigurePlaneGeometry->Map( point2D, point3D ); // Polygons (partially) outside of the image bounds can not be processed further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { float bounds[2] = {0,0}; bounds[0] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i0); bounds[1] = this->m_InternalImage->GetLargestPossibleRegion().GetSize().GetElement(i1); imageGeometry3D->WorldToIndex( point3D, point3D ); if (point3D[i0]<0) point3D[i0] = 0.0; else if (point3D[i0]>bounds[0]) point3D[i0] = bounds[0]-0.001; if (point3D[i1]<0) point3D[i1] = 0.0; else if (point3D[i1]>bounds[1]) point3D[i1] = bounds[1]-0.001; points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } else { imageGeometry3D->WorldToIndex( point3D, point3D ); // Add point to polyline array points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); numberOfPoints++; } } polyline->SetPoints( points ); points->Delete(); vtkIdType *ptIds = new vtkIdType[numberOfPoints]; for ( vtkIdType i = 0; i < numberOfPoints; ++i ) ptIds[i] = i; polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); // Extrude the generated contour polygon vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); extrudeFilter->SetInputData( polyline ); extrudeFilter->SetScaleFactor( 1 ); extrudeFilter->SetExtrusionTypeToNormalExtrusion(); extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); // Make a stencil from the extruded polygon vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); polyDataToImageStencil->SetInputConnection( extrudeFilter->GetOutputPort() ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< ItkUCharImageType > ImageImportType; typedef itk::VTKImageExport< ItkUCharImageType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( newMaskImage ); vtkImageImport *vtkImporter = vtkImageImport::New(); this->ConnectPipelines( itkExporter, vtkImporter ); vtkImporter->Update(); // Apply the generated image stencil to the input image vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(polyDataToImageStencil->GetOutputPort() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkImageExport *vtkExporter = vtkImageExport::New(); vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // calculate cropping bounding box m_InternalImageMask3D = itkImporter->GetOutput(); m_InternalImageMask3D->SetDirection(image->GetDirection()); itk::ImageRegionConstIterator itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); itk::ImageRegionIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); typename ImageType::SizeType lowersize = {{itk::NumericTraits::max(),itk::NumericTraits::max(),itk::NumericTraits::max()}}; typename ImageType::SizeType uppersize = {{0,0,0}}; while( !itmask.IsAtEnd() ) { if(itmask.Get() == 0) itimage.Set(0); else { typename ImageType::IndexType index = itimage.GetIndex(); typename ImageType::SizeType signedindex; signedindex[0] = index[0]; signedindex[1] = index[1]; signedindex[2] = index[2]; lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; } ++itmask; ++itimage; } typename ImageType::IndexType index; index[0] = lowersize[0]; index[1] = lowersize[1]; index[2] = lowersize[2]; typename ImageType::SizeType size; size[0] = uppersize[0] - lowersize[0] + 1; size[1] = uppersize[1] - lowersize[1] + 1; size[2] = uppersize[2] - lowersize[2] + 1; itk::ImageRegion<3> cropRegion = itk::ImageRegion<3>(index, size); // crop internal mask typedef itk::RegionOfInterestImageFilter< ItkUCharImageType, ItkUCharImageType > ROIMaskFilterType; typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); roi2->SetRegionOfInterest(cropRegion); roi2->SetInput(m_InternalImageMask3D); roi2->Update(); m_InternalImageMask3D = roi2->GetOutput(); Image::Pointer tmpImage = Image::New(); tmpImage->InitializeByItk(m_InternalImageMask3D.GetPointer()); tmpImage->SetVolume(m_InternalImageMask3D->GetBufferPointer()); Image::Pointer tmpImage2 = Image::New(); tmpImage2->InitializeByItk(m_PlanarFigureImage.GetPointer()); const BaseGeometry *pfImageGeometry3D = tmpImage2->GetGeometry( 0 ); const BaseGeometry *intImageGeometry3D = tmpImage->GetGeometry( 0 ); typedef itk::ImageRegionIteratorWithIndex IteratorType; IteratorType imageIterator (m_InternalImageMask3D, m_InternalImageMask3D->GetRequestedRegion()); imageIterator.GoToBegin(); while ( !imageIterator.IsAtEnd() ) { unsigned char val = imageIterator.Value(); if (val>0) { itk::Index<3> index = imageIterator.GetIndex(); Point3D point; point[0] = index[0]; point[1] = index[1]; point[2] = index[2]; intImageGeometry3D->IndexToWorld(point, point); pfImageGeometry3D->WorldToIndex(point, point); point[i0] += 0.5; point[i1] += 0.5; index[0] = point[0]; index[1] = point[1]; index[2] = point[2]; if (pfImageGeometry3D->IsIndexInside(index)) m_PlanarFigureImage->SetPixel(index, 1); } ++imageIterator; } // Clean up VTK objects polyline->Delete(); extrudeFilter->Delete(); polyDataToImageStencil->Delete(); vtkImporter->Delete(); imageStencilFilter->Delete(); //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? delete[] ptIds; } void QmitkFiberProcessingView::UpdateGui() { m_Controls->m_FibLabel->setText("mandatory"); m_Controls->m_PfLabel->setText("needed for extraction"); m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->m_RemoveButton->setEnabled(false); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); m_Controls->PFCompoANDButton->setEnabled(false); m_Controls->PFCompoORButton->setEnabled(false); m_Controls->PFCompoNOTButton->setEnabled(false); m_Controls->m_GenerateRoiImage->setEnabled(false); m_Controls->m_ExtractFibersButton->setEnabled(false); m_Controls->m_ModifyButton->setEnabled(false); m_Controls->m_CopyBundle->setEnabled(false); m_Controls->m_JoinBundles->setEnabled(false); m_Controls->m_SubstractBundles->setEnabled(false); // disable alle frames m_Controls->m_BundleWeightFrame->setVisible(false); m_Controls->m_ExtactionFramePF->setVisible(false); m_Controls->m_RemoveDirectionFrame->setVisible(false); m_Controls->m_RemoveLengthFrame->setVisible(false); m_Controls->m_RemoveCurvatureFrame->setVisible(false); m_Controls->m_RemoveByWeightFrame->setVisible(false); m_Controls->m_RemoveByDensityFrame->setVisible(false); m_Controls->m_SmoothFibersFrame->setVisible(false); m_Controls->m_CompressFibersFrame->setVisible(false); m_Controls->m_ColorFibersFrame->setVisible(false); m_Controls->m_MirrorFibersFrame->setVisible(false); m_Controls->m_MaskExtractionFrame->setVisible(false); m_Controls->m_ColorMapBox->setVisible(false); bool pfSelected = !m_SelectedPF.empty(); bool fibSelected = !m_SelectedFB.empty(); bool multipleFibsSelected = (m_SelectedFB.size()>1); bool maskSelected = m_RoiImageNode.IsNotNull(); bool imageSelected = m_SelectedImage.IsNotNull(); // toggle visibility of elements according to selected method switch ( m_Controls->m_ExtractionMethodBox->currentIndex() ) { case 0: m_Controls->m_ExtactionFramePF->setVisible(true); break; case 1: m_Controls->m_MaskExtractionFrame->setVisible(true); break; } switch ( m_Controls->m_RemovalMethodBox->currentIndex() ) { case 0: m_Controls->m_RemoveDirectionFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 1: m_Controls->m_RemoveLengthFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 2: m_Controls->m_RemoveCurvatureFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 3: break; case 4: break; case 5: m_Controls->m_RemoveByWeightFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; case 6: m_Controls->m_RemoveByDensityFrame->setVisible(true); if ( fibSelected ) m_Controls->m_RemoveButton->setEnabled(true); break; } switch ( m_Controls->m_ModificationMethodBox->currentIndex() ) { case 0: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 1: m_Controls->m_SmoothFibersFrame->setVisible(true); break; case 2: m_Controls->m_CompressFibersFrame->setVisible(true); break; case 3: m_Controls->m_ColorFibersFrame->setVisible(true); m_Controls->m_ColorMapBox->setVisible(true); break; case 4: m_Controls->m_MirrorFibersFrame->setVisible(true); if (m_SelectedSurfaces.size()>0) m_Controls->m_ModifyButton->setEnabled(true); break; case 5: m_Controls->m_BundleWeightFrame->setVisible(true); break; case 6: m_Controls->m_ColorFibersFrame->setVisible(true); break; case 7: m_Controls->m_ColorFibersFrame->setVisible(true); break; case 8: m_Controls->m_ColorFibersFrame->setVisible(true); break; } // are fiber bundles selected? if ( fibSelected ) { m_Controls->m_CopyBundle->setEnabled(true); m_Controls->m_ModifyButton->setEnabled(true); m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); m_Controls->m_FibLabel->setText(QString(m_SelectedFB.at(0)->GetName().c_str())); // one bundle and one planar figure needed to extract fibers if (pfSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==0) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_SelectedPF.at(0)->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } // more than two bundles needed to join/subtract if (multipleFibsSelected) { m_Controls->m_FibLabel->setText("multiple bundles selected"); m_Controls->m_JoinBundles->setEnabled(true); m_Controls->m_SubstractBundles->setEnabled(true); } if (maskSelected && m_Controls->m_ExtractionMethodBox->currentIndex()==1) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_RoiImageNode->GetName().c_str())); m_Controls->m_ExtractFibersButton->setEnabled(true); } if (maskSelected && (m_Controls->m_RemovalMethodBox->currentIndex()==3 || m_Controls->m_RemovalMethodBox->currentIndex()==4) ) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_PfLabel->setText(QString(m_RoiImageNode->GetName().c_str())); m_Controls->m_RemoveButton->setEnabled(true); } } // are planar figures selected? if (pfSelected) { if ( fibSelected || m_SelectedImage.IsNotNull()) m_Controls->m_GenerateRoiImage->setEnabled(true); if (m_SelectedPF.size() > 1) { m_Controls->PFCompoANDButton->setEnabled(true); m_Controls->PFCompoORButton->setEnabled(true); } else m_Controls->PFCompoNOTButton->setEnabled(true); } // is image selected if (imageSelected || maskSelected) { m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); } } void QmitkFiberProcessingView::NodeRemoved(const mitk::DataNode* node ) { for (auto fnode: m_SelectedFB) if (node == fnode) { m_SelectedFB.clear(); break; } berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } void QmitkFiberProcessingView::NodeAdded(const mitk::DataNode* ) { if (!m_Controls->m_InteractiveBox->isChecked()) { berry::IWorkbenchPart::Pointer nullPart; QList nodes; OnSelectionChanged(nullPart, nodes); } } void QmitkFiberProcessingView::OnEndInteraction() { if (m_Controls->m_InteractiveBox->isChecked()) ExtractWithPlanarFigure(true); } void QmitkFiberProcessingView::AddObservers() { typedef itk::SimpleMemberCommand< QmitkFiberProcessingView > SimpleCommandType; for (auto node : m_SelectedPF) { mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); if (figure!=nullptr) { figure->RemoveAllObservers(); // add observer for event when interaction with figure starts SimpleCommandType::Pointer endInteractionCommand = SimpleCommandType::New(); endInteractionCommand->SetCallbackFunction( this, &QmitkFiberProcessingView::OnEndInteraction); m_EndInteractionObserverTag = figure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), endInteractionCommand ); } } } void QmitkFiberProcessingView::RemoveObservers() { for (auto node : m_SelectedPF) { mitk::PlanarFigure* figure = dynamic_cast(node->GetData()); if (figure!=nullptr) figure->RemoveAllObservers(); } } void QmitkFiberProcessingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { RemoveObservers(); //reset existing Vectors containing FiberBundles and PlanarFigures from a previous selection std::vector lastSelectedFB = m_SelectedFB; m_SelectedFB.clear(); m_SelectedPF.clear(); m_SelectedSurfaces.clear(); m_SelectedImage = nullptr; m_RoiImageNode = nullptr; for (auto node: nodes) { if ( dynamic_cast(node->GetData()) ) m_SelectedFB.push_back(node); else if (dynamic_cast(node->GetData()) || dynamic_cast(node->GetData()) || dynamic_cast(node->GetData())) m_SelectedPF.push_back(node); else if (dynamic_cast(node->GetData())) { m_SelectedImage = dynamic_cast(node->GetData()); if (m_SelectedImage->GetDimension()==3) m_RoiImageNode = node; } else if (dynamic_cast(node->GetData())) m_SelectedSurfaces.push_back(dynamic_cast(node->GetData())); } // if we perform interactive fiber extraction, we want to avoid auto-selection of the extracted bundle if (m_SelectedFB.empty() && m_Controls->m_InteractiveBox->isChecked()) m_SelectedFB = lastSelectedFB; // if no fibers or surfaces are selected, select topmost if (m_SelectedFB.empty() && m_SelectedSurfaces.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedFB.clear(); m_SelectedFB.push_back(nodes->at(i)); } } } // if no plar figure is selected, select topmost if (m_SelectedPF.empty()) { int maxLayer = 0; itk::VectorContainer::ConstPointer nodes = this->GetDataStorage()->GetAll(); for (unsigned int i=0; iSize(); i++) if (dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData()) || dynamic_cast(nodes->at(i)->GetData())) { mitk::DataStorage::SetOfObjects::ConstPointer sources = GetDataStorage()->GetSources(nodes->at(i)); if (sources->Size()>0) continue; int layer = 0; nodes->at(i)->GetPropertyValue("layer", layer); if (layer>=maxLayer) { maxLayer = layer; m_SelectedPF.clear(); m_SelectedPF.push_back(nodes->at(i)); } } } AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::OnDrawPolygon() { mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); figure->ClosedOn(); this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); } void QmitkFiberProcessingView::OnDrawCircle() { mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_CircleCounter)); } void QmitkFiberProcessingView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *, mitk::BaseProperty* ) { // initialize figure's geometry with empty geometry mitk::PlaneGeometry::Pointer emptygeometry = mitk::PlaneGeometry::New(); figure->SetPlaneGeometry( emptygeometry ); //set desired data to DataNode where Planarfigure is stored mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetName(name.toStdString()); newNode->SetData(figure); newNode->SetBoolProperty("planarfigure.3drendering", true); newNode->SetBoolProperty("planarfigure.3drendering.fill", true); mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(newNode->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode(newNode); } // figure drawn on the topmost layer / image GetDataStorage()->Add(newNode ); RemoveObservers(); for(unsigned int i = 0; i < m_SelectedPF.size(); i++) m_SelectedPF[i]->SetSelected(false); newNode->SetSelected(true); m_SelectedPF.clear(); m_SelectedPF.push_back(newNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::ExtractWithPlanarFigure(bool interactive) { if ( m_SelectedFB.empty() || m_SelectedPF.empty() ){ QMessageBox::information( nullptr, "Warning", "No fibe bundle selected!"); return; } try { std::vector fiberBundles = m_SelectedFB; mitk::DataNode::Pointer planarFigure = m_SelectedPF.at(0); for (unsigned int i=0; i(fiberBundles.at(i)->GetData()); mitk::FiberBundle::Pointer extFB = fib->ExtractFiberSubset(planarFigure, GetDataStorage()); if (interactive && 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); } float op = 5.0/sqrt(fib->GetNumFibers()); float currentOp = 0; fiberBundles.at(i)->GetFloatProperty("opacity", currentOp); if (currentOp!=op) { fib->SetFiberColors(255, 255, 255); fiberBundles.at(i)->SetFloatProperty("opacity", op); fiberBundles.at(i)->SetBoolProperty("Fiber2DfadeEFX", false); } m_InteractiveNode->SetData(extFB); } else { if (extFB->GetNumFibers()<=0) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers."); continue; } mitk::DataNode::Pointer node; node = mitk::DataNode::New(); node->SetData(extFB); QString name(fiberBundles.at(i)->GetName().c_str()); name += "*"; node->SetName(name.toStdString()); fiberBundles.at(i)->SetVisibility(false); GetDataStorage()->Add(node); } } } catch(const std::out_of_range& ) { QMessageBox::warning( nullptr, "Fiber extraction failed", "Did you only create the planar figure, using the circle or polygon button, but forgot to actually place it in the image afterwards? \nAfter creating a planar figure, simply left-click at the desired position in the image or on the tractogram to place it."); } } void QmitkFiberProcessingView::GenerateAndComposite() { mitk::PlanarFigureComposite::Pointer PFCAnd = mitk::PlanarFigureComposite::New(); PFCAnd->setOperationType(mitk::PlanarFigureComposite::AND); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("AND"); newPFCNode->SetData(PFCAnd); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); RemoveObservers(); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::GenerateOrComposite() { mitk::PlanarFigureComposite::Pointer PFCOr = mitk::PlanarFigureComposite::New(); PFCOr->setOperationType(mitk::PlanarFigureComposite::OR); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("OR"); newPFCNode->SetData(PFCOr); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); UpdateGui(); } void QmitkFiberProcessingView::GenerateNotComposite() { mitk::PlanarFigureComposite::Pointer PFCNot = mitk::PlanarFigureComposite::New(); PFCNot->setOperationType(mitk::PlanarFigureComposite::NOT); mitk::DataNode::Pointer newPFCNode; newPFCNode = mitk::DataNode::New(); newPFCNode->SetName("NOT"); newPFCNode->SetData(PFCNot); RemoveObservers(); AddCompositeToDatastorage(newPFCNode, m_SelectedPF); m_SelectedPF.clear(); m_SelectedPF.push_back(newPFCNode); AddObservers(); UpdateGui(); } void QmitkFiberProcessingView::AddCompositeToDatastorage(mitk::DataNode::Pointer pfc, std::vector children, mitk::DataNode::Pointer parentNode ) { pfc->SetSelected(true); if (parentNode.IsNotNull()) GetDataStorage()->Add(pfc, parentNode); else GetDataStorage()->Add(pfc); for (auto child : children) { if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); newChild->SetBoolProperty("planarfigure.3drendering", true); newChild->SetBoolProperty("planarfigure.3drendering.fill", true); GetDataStorage()->Add(newChild, pfc); GetDataStorage()->Remove(child); } else if (dynamic_cast(child->GetData())) { mitk::DataNode::Pointer newChild; newChild = mitk::DataNode::New(); newChild->SetData(dynamic_cast(child->GetData())); newChild->SetName( child->GetName() ); std::vector< mitk::DataNode::Pointer > grandChildVector; mitk::DataStorage::SetOfObjects::ConstPointer grandchildren = GetDataStorage()->GetDerivations(child); for( mitk::DataStorage::SetOfObjects::const_iterator it = grandchildren->begin(); it != grandchildren->end(); ++it ) grandChildVector.push_back(*it); AddCompositeToDatastorage(newChild, grandChildVector, pfc); GetDataStorage()->Remove(child); } } UpdateGui(); } void QmitkFiberProcessingView::CopyBundles() { if ( m_SelectedFB.empty() ){ QMessageBox::information( nullptr, "Warning", "Select at least one fiber bundle!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least one fiber bundle!"; return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); mitk::FiberBundle::Pointer newFib = fib->GetDeepCopy(); node->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); name += "_copy"; mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newFib); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); } UpdateGui(); } void QmitkFiberProcessingView::JoinBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } m_SelectedFB.at(0)->SetVisibility(false); mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); std::vector< mitk::FiberBundle::Pointer > tractograms; for (unsigned int i=1; iSetVisibility(false); tractograms.push_back(dynamic_cast(m_SelectedFB.at(i)->GetData())); } newBundle = newBundle->AddBundles(tractograms); mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName("Joined_Tractograms"); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } void QmitkFiberProcessingView::SubstractBundles() { if ( m_SelectedFB.size()<2 ){ QMessageBox::information( nullptr, "Warning", "Select at least two fiber bundles!"); MITK_WARN("QmitkFiberProcessingView") << "Select at least two fiber bundles!"; return; } mitk::FiberBundle::Pointer newBundle = dynamic_cast(m_SelectedFB.at(0)->GetData()); m_SelectedFB.at(0)->SetVisibility(false); QString name(""); name += QString(m_SelectedFB.at(0)->GetName().c_str()); for (unsigned int i=1; iSubtractBundle(dynamic_cast(m_SelectedFB.at(i)->GetData())); if (newBundle.IsNull()) break; name += "-"+QString(m_SelectedFB.at(i)->GetName().c_str()); m_SelectedFB.at(i)->SetVisibility(false); } if (newBundle.IsNull()) { QMessageBox::information(nullptr, "No output generated:", "The resulting fiber bundle contains no fibers. Did you select the fiber bundles in the correct order? X-Y is not equal to Y-X!"); return; } mitk::DataNode::Pointer fbNode = mitk::DataNode::New(); fbNode->SetData(newBundle); fbNode->SetName(name.toStdString()); fbNode->SetVisibility(true); GetDataStorage()->Add(fbNode); UpdateGui(); } void QmitkFiberProcessingView::ResampleSelectedBundlesSpline() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleSpline(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::ResampleSelectedBundlesLinear() { double factor = this->m_Controls->m_SmoothFibersBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ResampleLinear(factor); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::CompressSelectedBundles() { double factor = this->m_Controls->m_ErrorThresholdBox->value(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->Compress(factor); fib->ColorFibersByOrientation(); } RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkFiberProcessingView::DoImageColorCoding() { if (m_Controls->m_ColorMapBox->GetSelectedNode().IsNull()) { QMessageBox::information(nullptr, "Bundle coloring aborted:", "No image providing the scalar values for coloring the selected bundle available."); return; } for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByScalarMap(dynamic_cast(m_Controls->m_ColorMapBox->GetSelectedNode()->GetData()), m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoCurvatureColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByCurvature(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoLengthColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByLength(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::DoWeightColorCoding() { for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); fib->ColorFibersByFiberWeights(m_Controls->m_FiberOpacityBox->isChecked(), m_Controls->m_NormalizeColorValues->isChecked()); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } void QmitkFiberProcessingView::MirrorFibers() { unsigned int axis = this->m_Controls->m_MirrorFibersBox->currentIndex(); for (auto node : m_SelectedFB) { mitk::FiberBundle::Pointer fib = dynamic_cast(node->GetData()); if (m_SelectedImage.IsNotNull()) - fib->SetReferenceGeometry(m_SelectedImage->GetGeometry()); + fib->SetTrackVisHeader(m_SelectedImage->GetGeometry()); fib->MirrorFibers(axis); } - for (auto surf : m_SelectedSurfaces) { vtkSmartPointer poly = surf->GetVtkPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); for (int i=0; iGetNumberOfPoints(); i++) { double* point = poly->GetPoint(i); point[axis] *= -1; vtkNewPoints->InsertNextPoint(point); } poly->SetPoints(vtkNewPoints); surf->CalculateBoundingBox(); } if (auto renderWindowPart = this->GetRenderWindowPart()) { renderWindowPart->RequestUpdate(); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp index e78d3cf..baa5202 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp @@ -1,608 +1,608 @@ /*=================================================================== 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 // Qmitk #include "QmitkGibbsTrackingView.h" // Qt #include #include #include #include // MITK #include #include #include #include #include #include #include #include #include #include // ITK #include #include #include // MISC #include QmitkTrackingWorker::QmitkTrackingWorker(QmitkGibbsTrackingView* view) : m_View(view) { } void QmitkTrackingWorker::run() { m_View->m_GlobalTracker = QmitkGibbsTrackingView::GibbsTrackingFilterType::New(); m_View->m_GlobalTracker->SetOdfImage(m_View->m_ItkOdfImage); m_View->m_GlobalTracker->SetTensorImage(m_View->m_ItkTensorImage); m_View->m_GlobalTracker->SetMaskImage(m_View->m_MaskImage); m_View->m_GlobalTracker->SetStartTemperature((float)m_View->m_Controls->m_StartTempSlider->value()/100); m_View->m_GlobalTracker->SetEndTemperature((float)m_View->m_Controls->m_EndTempSlider->value()/10000); m_View->m_GlobalTracker->SetIterations(m_View->m_Controls->m_IterationsBox->text().toDouble()); m_View->m_GlobalTracker->SetParticleWeight((float)m_View->m_Controls->m_ParticleWeightSlider->value()/10000); m_View->m_GlobalTracker->SetParticleWidth((float)(m_View->m_Controls->m_ParticleWidthSlider->value())/10); m_View->m_GlobalTracker->SetParticleLength((float)(m_View->m_Controls->m_ParticleLengthSlider->value())/10); m_View->m_GlobalTracker->SetInexBalance((float)m_View->m_Controls->m_InExBalanceSlider->value()/10); m_View->m_GlobalTracker->SetMinFiberLength(m_View->m_Controls->m_FiberLengthSlider->value()); m_View->m_GlobalTracker->SetCurvatureThreshold(cos((float)m_View->m_Controls->m_CurvatureThresholdSlider->value()*itk::Math::pi/180)); m_View->m_GlobalTracker->SetRandomSeed(m_View->m_Controls->m_RandomSeedSlider->value()); try{ m_View->m_GlobalTracker->Update(); } catch( mitk::Exception& e ) { MITK_ERROR << "Internal error occured: " << e.what() << "\nAborting"; } m_View->m_TrackingThread.quit(); } const std::string QmitkGibbsTrackingView::VIEW_ID = "org.mitk.views.gibbstracking"; QmitkGibbsTrackingView::QmitkGibbsTrackingView() : QmitkAbstractView() , m_Controls( 0 ) , m_TrackingNode(nullptr) , m_FiberBundle(nullptr) , m_MaskImage(nullptr) , m_TensorImage(nullptr) , m_OdfImage(nullptr) , m_ItkOdfImage(nullptr) , m_ItkTensorImage(nullptr) , m_ImageNode(nullptr) , m_MaskImageNode(nullptr) , m_FiberBundleNode(nullptr) , m_ThreadIsRunning(false) , m_ElapsedTime(0) , m_GlobalTracker(nullptr) , m_TrackingWorker(this) { 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); } QmitkGibbsTrackingView::~QmitkGibbsTrackingView() { if (m_GlobalTracker.IsNull()) return; m_GlobalTracker->SetAbortTracking(true); m_TrackingThread.wait(); } // update tracking status and generate fiber bundle void QmitkGibbsTrackingView::TimerUpdate() { UpdateTrackingStatus(); GenerateFiberBundle(); } // tell global tractography filter to stop after current step void QmitkGibbsTrackingView::StopGibbsTracking() { if (m_GlobalTracker.IsNull()) return; m_GlobalTracker->SetAbortTracking(true); m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStop->setText("Stopping Tractography ..."); m_TrackingNode = nullptr; } // update gui elements and generate fiber bundle after tracking is finished void QmitkGibbsTrackingView::AfterThread() { m_ThreadIsRunning = false; m_TrackingTimer->stop(); UpdateGUI(); if( !m_GlobalTracker->GetIsInValidState() ) { QMessageBox::critical( nullptr, "Gibbs Tracking", "An internal error occured. Tracking aborted.\n Please check the log for details." ); m_FiberBundleNode = nullptr; return; } UpdateTrackingStatus(); if(m_Controls->m_ParticleWeightSlider->value()==0) { m_Controls->m_ParticleWeightLabel->setText(QString::number(m_GlobalTracker->GetParticleWeight())); m_Controls->m_ParticleWeightSlider->setValue(m_GlobalTracker->GetParticleWeight()*10000); } if(m_Controls->m_ParticleWidthSlider->value()==0) { m_Controls->m_ParticleWidthLabel->setText(QString::number(m_GlobalTracker->GetParticleWidth())); m_Controls->m_ParticleWidthSlider->setValue(m_GlobalTracker->GetParticleWidth()*10); } if(m_Controls->m_ParticleLengthSlider->value()==0) { m_Controls->m_ParticleLengthLabel->setText(QString::number(m_GlobalTracker->GetParticleLength())); m_Controls->m_ParticleLengthSlider->setValue(m_GlobalTracker->GetParticleLength()*10); } GenerateFiberBundle(); m_FiberBundleNode = 0; m_GlobalTracker = 0; // images not needed anymore ( relevant only for computation ) // we need to release them to remove the memory access block created through CastToItk<> calls this->m_ItkOdfImage = 0; this->m_ItkTensorImage = 0; } // start tracking timer and update gui elements before tracking is started void QmitkGibbsTrackingView::BeforeThread() { m_ThreadIsRunning = true; m_TrackingTime = QTime::currentTime(); m_ElapsedTime = 0; m_TrackingTimer->start(1000); UpdateGUI(); } // setup gui elements and signal/slot connections void QmitkGibbsTrackingView::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::QmitkGibbsTrackingViewControls; m_Controls->setupUi( parent ); AdvancedSettings(); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); connect( m_Controls->m_TrackingStop, SIGNAL(clicked()), this, SLOT(StopGibbsTracking()) ); connect( m_Controls->m_TrackingStart, SIGNAL(clicked()), this, SLOT(StartGibbsTracking()) ); connect( m_Controls->m_AdvancedSettingsCheckbox, SIGNAL(clicked()), this, SLOT(AdvancedSettings()) ); connect( m_Controls->m_SaveTrackingParameters, SIGNAL(clicked()), this, SLOT(SaveTrackingParameters()) ); connect( m_Controls->m_LoadTrackingParameters, SIGNAL(clicked()), this, SLOT(LoadTrackingParameters()) ); connect( m_Controls->m_ParticleWidthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWidth(int)) ); connect( m_Controls->m_ParticleLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleLength(int)) ); connect( m_Controls->m_InExBalanceSlider, SIGNAL(valueChanged(int)), this, SLOT(SetInExBalance(int)) ); connect( m_Controls->m_FiberLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetFiberLength(int)) ); connect( m_Controls->m_ParticleWeightSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWeight(int)) ); connect( m_Controls->m_StartTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetStartTemp(int)) ); connect( m_Controls->m_EndTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetEndTemp(int)) ); connect( m_Controls->m_CurvatureThresholdSlider, SIGNAL(valueChanged(int)), this, SLOT(SetCurvatureThreshold(int)) ); connect( m_Controls->m_RandomSeedSlider, SIGNAL(valueChanged(int)), this, SLOT(SetRandomSeed(int)) ); connect( m_Controls->m_OutputFileButton, SIGNAL(clicked()), this, SLOT(SetOutputFile()) ); m_Controls->m_InputImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isOdfImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTensorImagePredicate = mitk::TNodePredicateDataType::New(); m_Controls->m_InputImageBox->SetPredicate( mitk::NodePredicateOr::New(isOdfImagePredicate, isTensorImagePredicate) ); m_Controls->m_MaskImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskImageBox->SetZeroEntryText("--"); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateDimension::Pointer is3D = mitk::NodePredicateDimension::New(3); m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, mitk::NodePredicateAnd::New(isImagePredicate, is3D)) ); connect( (QObject*)(m_Controls->m_MaskImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGUI())); connect( (QObject*)(m_Controls->m_InputImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGUI())); } } void QmitkGibbsTrackingView::SetFocus() { m_Controls->m_TrackingStart->setFocus(); } void QmitkGibbsTrackingView::SetInExBalance(int value) { m_Controls->m_InExBalanceLabel->setText(QString::number((float)value/10)); } void QmitkGibbsTrackingView::SetFiberLength(int value) { m_Controls->m_FiberLengthLabel->setText(QString::number(value)+"mm"); } void QmitkGibbsTrackingView::SetRandomSeed(int value) { if (value>=0) m_Controls->m_RandomSeedLabel->setText(QString::number(value)); else m_Controls->m_RandomSeedLabel->setText("auto"); } void QmitkGibbsTrackingView::SetParticleWeight(int value) { if (value>0) m_Controls->m_ParticleWeightLabel->setText(QString::number((float)value/10000)); else m_Controls->m_ParticleWeightLabel->setText("auto"); } void QmitkGibbsTrackingView::SetStartTemp(int value) { m_Controls->m_StartTempLabel->setText(QString::number((float)value/100)); } void QmitkGibbsTrackingView::SetEndTemp(int value) { m_Controls->m_EndTempLabel->setText(QString::number((float)value/10000)); } void QmitkGibbsTrackingView::SetParticleWidth(int value) { if (value>0) m_Controls->m_ParticleWidthLabel->setText(QString::number((float)value/10)+" mm"); else m_Controls->m_ParticleWidthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetParticleLength(int value) { if (value>0) m_Controls->m_ParticleLengthLabel->setText(QString::number((float)value/10)+" mm"); else m_Controls->m_ParticleLengthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetCurvatureThreshold(int value) { m_Controls->m_CurvatureThresholdLabel->setText(QString::number(value)+"°"); } // called if datamanager selection changes void QmitkGibbsTrackingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { UpdateGUI(); } void QmitkGibbsTrackingView::NodeRemoved(const mitk::DataNode * node) { if (m_ThreadIsRunning) { if (node==m_TrackingNode.GetPointer()) { StopGibbsTracking(); } } } // update gui elements displaying trackings status void QmitkGibbsTrackingView::UpdateTrackingStatus() { if (m_GlobalTracker.IsNull()) return; m_ElapsedTime += m_TrackingTime.elapsed()/1000; m_TrackingTime.restart(); unsigned long hours = m_ElapsedTime/3600; unsigned long minutes = (m_ElapsedTime%3600)/60; unsigned long seconds = m_ElapsedTime%60; m_Controls->m_ProposalAcceptance->setText(QString::number(m_GlobalTracker->GetProposalAcceptance()*100)+"%"); m_Controls->m_TrackingTimeLabel->setText( QString::number(hours)+QString("h ")+QString::number(minutes)+QString("m ")+QString::number(seconds)+QString("s") ); m_Controls->m_NumConnectionsLabel->setText( QString::number(m_GlobalTracker->GetNumConnections()) ); m_Controls->m_NumParticlesLabel->setText( QString::number(m_GlobalTracker->GetNumParticles()) ); m_Controls->m_CurrentStepLabel->setText( QString::number(100*m_GlobalTracker->GetCurrentIteration()/m_GlobalTracker->GetIterations())+"%" ); m_Controls->m_AcceptedFibersLabel->setText( QString::number(m_GlobalTracker->GetNumAcceptedFibers()) ); } // update gui elements (enable/disable elements and set tooltips) void QmitkGibbsTrackingView::UpdateGUI() { if (!m_ThreadIsRunning && m_Controls->m_InputImageBox->GetSelectedNode().IsNotNull()) { m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStart->setEnabled(true); m_Controls->m_LoadTrackingParameters->setEnabled(true); m_Controls->m_IterationsBox->setEnabled(true); m_Controls->m_AdvancedFrame->setEnabled(true); m_Controls->m_TrackingStop->setText("Stop Tractography"); m_Controls->m_TrackingStart->setToolTip("Start tractography. No further change of parameters possible."); m_Controls->m_TrackingStop->setToolTip(""); m_Controls->m_MaskImageBox->setEnabled(true); m_Controls->m_InputImageBox->setEnabled(true); } else if (!m_ThreadIsRunning) { m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStart->setEnabled(false); m_Controls->m_LoadTrackingParameters->setEnabled(true); m_Controls->m_IterationsBox->setEnabled(true); m_Controls->m_AdvancedFrame->setEnabled(true); m_Controls->m_TrackingStop->setText("Stop Tractography"); m_Controls->m_TrackingStart->setToolTip("No ODF image selected."); m_Controls->m_TrackingStop->setToolTip(""); m_Controls->m_MaskImageBox->setEnabled(true); m_Controls->m_InputImageBox->setEnabled(true); } else { m_Controls->m_TrackingStop->setEnabled(true); m_Controls->m_TrackingStart->setEnabled(false); m_Controls->m_LoadTrackingParameters->setEnabled(false); m_Controls->m_IterationsBox->setEnabled(false); m_Controls->m_AdvancedFrame->setEnabled(false); m_Controls->m_AdvancedFrame->setVisible(false); m_Controls->m_AdvancedSettingsCheckbox->setChecked(false); m_Controls->m_TrackingStart->setToolTip("Tracking in progress."); m_Controls->m_TrackingStop->setToolTip("Stop tracking and display results."); m_Controls->m_MaskImageBox->setEnabled(false); m_Controls->m_InputImageBox->setEnabled(false); } } // show/hide advanced settings frame void QmitkGibbsTrackingView::AdvancedSettings() { m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); } // check for mask and odf and start tracking thread void QmitkGibbsTrackingView::StartGibbsTracking() { if(m_ThreadIsRunning) { MITK_WARN("QmitkGibbsTrackingView")<<"Thread already running!"; return; } m_GlobalTracker = nullptr; if (m_Controls->m_InputImageBox->GetSelectedNode().IsNull()) { QMessageBox::information( nullptr, "Warning", "Please load and select a Odf image before starting image processing."); return; } m_ImageNode = m_Controls->m_InputImageBox->GetSelectedNode(); if (dynamic_cast(m_ImageNode->GetData())) m_OdfImage = dynamic_cast(m_ImageNode->GetData()); else if (dynamic_cast(m_ImageNode->GetData())) m_TensorImage = dynamic_cast(m_ImageNode->GetData()); if (m_OdfImage.IsNull() && m_TensorImage.IsNull()) return; // cast odf to itk m_TrackingNode = m_ImageNode; m_ItkTensorImage = nullptr; m_ItkOdfImage = nullptr; m_MaskImage = nullptr; if (m_OdfImage.IsNotNull()) { m_ItkOdfImage = ItkOdfImgType::New(); mitk::CastToItkImage(m_OdfImage, m_ItkOdfImage); } else { m_ItkTensorImage = ItkTensorImage::New(); mitk::CastToItkImage(m_TensorImage, m_ItkTensorImage); } // mask image found? // catch exceptions thrown by the itkAccess macros try{ if(m_Controls->m_MaskImageBox->GetSelectedNode().IsNotNull()) { m_MaskImageNode = m_Controls->m_MaskImageBox->GetSelectedNode(); if (dynamic_cast(m_MaskImageNode->GetData())) mitk::CastToItkImage(dynamic_cast(m_MaskImageNode->GetData()), m_MaskImage); } } catch(...){}; // start worker thread m_TrackingThread.start(QThread::LowestPriority); } // generate mitkFiberBundle from tracking filter output void QmitkGibbsTrackingView::GenerateFiberBundle() { if (m_GlobalTracker.IsNull() || (!(m_Controls->m_VisualizationCheckbox->isChecked() || m_Controls->m_VisualizeOnceButton->isChecked()) && m_ThreadIsRunning)) return; if (m_Controls->m_VisualizeOnceButton->isChecked()) m_Controls->m_VisualizeOnceButton->setChecked(false); vtkSmartPointer fiberBundle = m_GlobalTracker->GetFiberBundle(); if ( m_GlobalTracker->GetNumAcceptedFibers()==0 ) return; m_FiberBundle = mitk::FiberBundle::New(fiberBundle); - m_FiberBundle->SetReferenceGeometry(dynamic_cast(m_ImageNode->GetData())->GetGeometry()); + m_FiberBundle->SetTrackVisHeader(dynamic_cast(m_ImageNode->GetData())->GetGeometry()); if (m_FiberBundleNode.IsNotNull()){ GetDataStorage()->Remove(m_FiberBundleNode); m_FiberBundleNode = 0; } m_GlobalTracker->SetDicomProperties(m_FiberBundle); m_FiberBundleNode = mitk::DataNode::New(); m_FiberBundleNode->SetData(m_FiberBundle); QString name("FiberBundle_"); name += m_ImageNode->GetName().c_str(); name += "_Gibbs"; m_FiberBundleNode->SetName(name.toStdString()); m_FiberBundleNode->SetVisibility(true); if (!m_OutputFileName.isEmpty() && !m_ThreadIsRunning) { try { mitk::IOUtil::Save(m_FiberBundle.GetPointer(),m_OutputFileName.toStdString()); QMessageBox::information(nullptr, "Fiber bundle saved to", m_OutputFileName); } catch (itk::ExceptionObject &ex) { QMessageBox::information(nullptr, "Fiber bundle could not be saved", QString("%1\n%2\n%3\n%4\n%5\n%6").arg(ex.GetNameOfClass()).arg(ex.GetFile()).arg(ex.GetLine()).arg(ex.GetLocation()).arg(ex.what()).arg(ex.GetDescription())); } } if(m_ImageNode.IsNull()) GetDataStorage()->Add(m_FiberBundleNode); else GetDataStorage()->Add(m_FiberBundleNode, m_ImageNode); } void QmitkGibbsTrackingView::SetOutputFile() { // SELECT FOLDER DIALOG m_OutputFileName = QFileDialog::getSaveFileName(0, tr("Set file name"), QDir::currentPath()+"/FiberBundle.fib", tr("Fiber Bundle (*.fib)") ); if (m_OutputFileName.isEmpty()) m_Controls->m_OutputFileLabel->setText("N/A"); else m_Controls->m_OutputFileLabel->setText(m_OutputFileName); } // save current tracking paramters as xml file (.gtp) void QmitkGibbsTrackingView::SaveTrackingParameters() { TiXmlDocument documentXML; TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); documentXML.LinkEndChild( declXML ); TiXmlElement* mainXML = new TiXmlElement("global_tracking_parameter_file"); mainXML->SetAttribute("file_version", "0.1"); documentXML.LinkEndChild(mainXML); TiXmlElement* paramXML = new TiXmlElement("parameter_set"); paramXML->SetAttribute("iterations", m_Controls->m_IterationsBox->text().toStdString()); paramXML->SetAttribute("particle_length", QString::number((float)m_Controls->m_ParticleLengthSlider->value()/10).toStdString()); paramXML->SetAttribute("particle_width", QString::number((float)m_Controls->m_ParticleWidthSlider->value()/10).toStdString()); paramXML->SetAttribute("particle_weight", QString::number((float)m_Controls->m_ParticleWeightSlider->value()/10000).toStdString()); paramXML->SetAttribute("temp_start", QString::number((float)m_Controls->m_StartTempSlider->value()/100).toStdString()); paramXML->SetAttribute("temp_end", QString::number((float)m_Controls->m_EndTempSlider->value()/10000).toStdString()); paramXML->SetAttribute("inexbalance", QString::number((float)m_Controls->m_InExBalanceSlider->value()/10).toStdString()); paramXML->SetAttribute("fiber_length", QString::number(m_Controls->m_FiberLengthSlider->value()).toStdString()); paramXML->SetAttribute("curvature_threshold", QString::number(m_Controls->m_CurvatureThresholdSlider->value()).toStdString()); mainXML->LinkEndChild(paramXML); QString filename = QFileDialog::getSaveFileName( 0, tr("Save Parameters"), QDir::currentPath()+"/param.gtp", tr("Global Tracking Parameters (*.gtp)") ); if(filename.isEmpty() || filename.isNull()) return; if(!filename.endsWith(".gtp")) filename += ".gtp"; documentXML.SaveFile( filename.toStdString() ); } // load current tracking paramters from xml file (.gtp) void QmitkGibbsTrackingView::LoadTrackingParameters() { QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QDir::currentPath(), tr("Global Tracking Parameters (*.gtp)") ); if(filename.isEmpty() || filename.isNull()) return; TiXmlDocument doc( filename.toStdString() ); doc.LoadFile(); TiXmlHandle hDoc(&doc); TiXmlElement* pElem; TiXmlHandle hRoot(0); pElem = hDoc.FirstChildElement().Element(); hRoot = TiXmlHandle(pElem); pElem = hRoot.FirstChildElement("parameter_set").Element(); QString iterations(pElem->Attribute("iterations")); m_Controls->m_IterationsBox->setText(iterations); QString particleLength(pElem->Attribute("particle_length")); float pLength = particleLength.toFloat(); QString particleWidth(pElem->Attribute("particle_width")); float pWidth = particleWidth.toFloat(); if (pLength==0) m_Controls->m_ParticleLengthLabel->setText("auto"); else m_Controls->m_ParticleLengthLabel->setText(particleLength+" mm"); if (pWidth==0) m_Controls->m_ParticleWidthLabel->setText("auto"); else m_Controls->m_ParticleWidthLabel->setText(particleWidth+" mm"); m_Controls->m_ParticleWidthSlider->setValue(pWidth*10); m_Controls->m_ParticleLengthSlider->setValue(pLength*10); QString partWeight(pElem->Attribute("particle_weight")); m_Controls->m_ParticleWeightSlider->setValue(partWeight.toFloat()*10000); m_Controls->m_ParticleWeightLabel->setText(partWeight); QString startTemp(pElem->Attribute("temp_start")); m_Controls->m_StartTempSlider->setValue(startTemp.toFloat()*100); m_Controls->m_StartTempLabel->setText(startTemp); QString endTemp(pElem->Attribute("temp_end")); m_Controls->m_EndTempSlider->setValue(endTemp.toFloat()*10000); m_Controls->m_EndTempLabel->setText(endTemp); QString inExBalance(pElem->Attribute("inexbalance")); m_Controls->m_InExBalanceSlider->setValue(inExBalance.toFloat()*10); m_Controls->m_InExBalanceLabel->setText(inExBalance); QString fiberLength(pElem->Attribute("fiber_length")); m_Controls->m_FiberLengthSlider->setValue(fiberLength.toInt()); m_Controls->m_FiberLengthLabel->setText(fiberLength+"mm"); QString curvThres(pElem->Attribute("curvature_threshold")); m_Controls->m_CurvatureThresholdSlider->setValue(curvThres.toInt()); m_Controls->m_CurvatureThresholdLabel->setText(curvThres+"°"); } 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 3bf1d89..78db5ee 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,1159 +1,1159 @@ /*=================================================================== 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->m_LoadParametersButton, SIGNAL(clicked()), this, SLOT(LoadParameters()) ); 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_MaxTractLengthBox, 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_PriorFlipXBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_PriorFlipYBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_PriorFlipZBox, 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()) ); connect( m_Controls->m_PeakJitterBox, SIGNAL(editingFinished()), 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_MaxTractLengthBox->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(); m_Controls->m_PeakJitterBox->editingFinished(); StartStopTrackingGui(false); } m_ParameterFile = QDir::currentPath()+"/param.stp"; UpdateGui(); } void QmitkStreamlineTrackingView::ParametersToGui(mitk::StreamlineTractographyParameters& params) { m_Controls->m_SeedRadiusBox->setValue(params.m_InteractiveRadiusMm); m_Controls->m_NumSeedsBox->setValue(params.m_NumInteractiveSeeds); m_Controls->m_InteractiveBox->setChecked(params.m_EnableInteractive); m_Controls->m_FiberErrorBox->setValue(params.m_Compression); m_Controls->m_ResampleFibersBox->setChecked(params.m_CompressFibers); m_Controls->m_SeedRadiusBox->setValue(params.m_InteractiveRadiusMm); m_Controls->m_NumFibersBox->setValue(params.m_MaxNumFibers); m_Controls->m_ScalarThresholdBox->setValue(params.m_Cutoff); m_Controls->m_fBox->setValue(params.m_F); m_Controls->m_gBox->setValue(params.m_G); m_Controls->m_OdfCutoffBox->setValue(params.m_OdfCutoff); m_Controls->m_SharpenOdfsBox->setChecked(params.m_SharpenOdfs); m_Controls->m_PriorWeightBox->setValue(params.m_Weight); m_Controls->m_PriorAsMaskBox->setChecked(params.m_RestrictToPrior); m_Controls->m_NewDirectionsFromPriorBox->setChecked(params.m_NewDirectionsFromPrior); m_Controls->m_PriorFlipXBox->setChecked(params.m_PriorFlipX); m_Controls->m_PriorFlipYBox->setChecked(params.m_PriorFlipY); m_Controls->m_PriorFlipZBox->setChecked(params.m_PriorFlipZ); m_Controls->m_FlipXBox->setChecked(params.m_FlipX); m_Controls->m_FlipYBox->setChecked(params.m_FlipY); m_Controls->m_FlipZBox->setChecked(params.m_FlipZ); m_Controls->m_InterpolationBox->setChecked(params.m_InterpolateTractographyData); m_Controls->m_MaskInterpolationBox->setChecked(params.m_InterpolateRoiImages); m_Controls->m_SeedsPerVoxelBox->setValue(params.m_SeedsPerVoxel); m_Controls->m_StepSizeBox->setValue(params.GetStepSizeVox()); m_Controls->m_SamplingDistanceBox->setValue(params.GetSamplingDistanceVox()); m_Controls->m_StopVotesBox->setChecked(params.m_StopVotes); m_Controls->m_FrontalSamplesBox->setChecked(params.m_OnlyForwardSamples); m_Controls->m_TrialsPerSeedBox->setValue(params.m_TrialsPerSeed); m_Controls->m_NumSamplesBox->setValue(params.m_NumSamples); m_Controls->m_LoopCheckBox->setValue(params.GetLoopCheckDeg()); m_Controls->m_AngularThresholdBox->setValue(params.GetAngularThresholdDeg()); m_Controls->m_MinTractLengthBox->setValue(params.m_MinTractLengthMm); m_Controls->m_MaxTractLengthBox->setValue(params.m_MaxTractLengthMm); m_Controls->m_OutputProbMap->setChecked(params.m_OutputProbMap); m_Controls->m_FixSeedBox->setChecked(params.m_FixRandomSeed); m_Controls->m_PeakJitterBox->setValue(params.m_PeakJitter); switch (params.m_Mode) { case mitk::TrackingDataHandler::MODE::DETERMINISTIC: m_Controls->m_ModeBox->setCurrentIndex(0); break; case mitk::TrackingDataHandler::MODE::PROBABILISTIC: m_Controls->m_ModeBox->setCurrentIndex(1); break; } switch (params.m_EpConstraints) { case itk::StreamlineTrackingFilter::EndpointConstraints::NONE: m_Controls->m_EpConstraintsBox->setCurrentIndex(0); break; case itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET: m_Controls->m_EpConstraintsBox->setCurrentIndex(1); break; case itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_TARGET_LABELDIFF: m_Controls->m_EpConstraintsBox->setCurrentIndex(2); break; case itk::StreamlineTrackingFilter::EndpointConstraints::EPS_IN_SEED_AND_TARGET: m_Controls->m_EpConstraintsBox->setCurrentIndex(3); break; case itk::StreamlineTrackingFilter::EndpointConstraints::MIN_ONE_EP_IN_TARGET: m_Controls->m_EpConstraintsBox->setCurrentIndex(4); break; case itk::StreamlineTrackingFilter::EndpointConstraints::ONE_EP_IN_TARGET: m_Controls->m_EpConstraintsBox->setCurrentIndex(5); break; case itk::StreamlineTrackingFilter::EndpointConstraints::NO_EP_IN_TARGET: m_Controls->m_EpConstraintsBox->setCurrentIndex(6); break; } } std::shared_ptr QmitkStreamlineTrackingView::GetParametersFromGui() { // NOT IN GUI // unsigned int m_NumPreviousDirections = 1; // bool m_AvoidStop = true; // bool m_RandomSampling = false; // float m_DeflectionMod = 1.0; // bool m_ApplyDirectionMatrix = false; std::shared_ptr params = std::make_shared(); params->m_InteractiveRadiusMm = m_Controls->m_SeedRadiusBox->value(); params->m_NumInteractiveSeeds = m_Controls->m_NumSeedsBox->value(); params->m_EnableInteractive = m_Controls->m_InteractiveBox->isChecked(); params->m_Compression = m_Controls->m_FiberErrorBox->value(); params->m_CompressFibers = m_Controls->m_ResampleFibersBox->isChecked(); params->m_InteractiveRadiusMm = 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_PriorFlipX = m_Controls->m_PriorFlipXBox->isChecked(); params->m_PriorFlipY = m_Controls->m_PriorFlipYBox->isChecked(); params->m_PriorFlipZ = m_Controls->m_PriorFlipZBox->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->SetSamplingDistanceVox(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_MinTractLengthMm = m_Controls->m_MinTractLengthBox->value(); params->m_MaxTractLengthMm = m_Controls->m_MaxTractLengthBox->value(); params->m_OutputProbMap = m_Controls->m_OutputProbMap->isChecked(); params->m_FixRandomSeed = m_Controls->m_FixSeedBox->isChecked(); params->m_PeakJitter = static_cast(m_Controls->m_PeakJitterBox->value()); 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; 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 Tractography Parameters"), m_ParameterFile, tr("Streamline Tractography Parameters (*.stp)") ); if(filename.isEmpty() || filename.isNull()) return; m_ParameterFile = filename; auto params = GetParametersFromGui(); params->SaveParameters(m_ParameterFile.toStdString()); } void QmitkStreamlineTrackingView::LoadParameters() { QString filename = QFileDialog::getOpenFileName( 0, tr("Load Tractography Parameters"), m_ParameterFile, tr("Streamline Tractography Parameters (*.stp)") ); if(filename.isEmpty() || filename.isNull()) return; m_ParameterFile = filename; mitk::StreamlineTractographyParameters params; params.LoadParameters(m_ParameterFile.toStdString()); ParametersToGui(params); } 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 (!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()); + fib->SetTrackVisHeader(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", params->GetMinVoxelSizeMm()/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", params->GetMinVoxelSizeMm()/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", params->GetMinVoxelSizeMm()/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); m_Controls->m_PeakJitterBox->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); } // stuff that is only important for probabilistic tractography if (m_Controls->m_ModeBox->currentIndex()==1) { m_Controls->m_TrialsPerSeedBox->setEnabled(true); m_Controls->m_TrialsPerSeedLabel->setEnabled(true); m_Controls->m_PeakJitterBox->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 (params->m_EpConstraints == itk::StreamlineTrackingFilter::EndpointConstraints::NONE) m_Tracker->SetTargetRegions(nullptr); 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)->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); } } } } 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); } } } 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) { 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); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } else { 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); 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_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingHandler)->SetPeakImage(mitk::convert::GetItkPeakFromPeakImage(dynamic_cast(m_InputImageNodes.at(0)->GetData()))); } } 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); m_LastPrior = m_Controls->m_PriorImageSelectionWidget->GetSelectedNode(); } prior_params->m_FlipX = m_Controls->m_PriorFlipXBox->isChecked(); prior_params->m_FlipY = m_Controls->m_PriorFlipYBox->isChecked(); prior_params->m_FlipZ = m_Controls->m_PriorFlipZBox->isChecked(); m_TrackingPriorHandler->SetParameters(prior_params); m_Tracker->SetTrackingPriorHandler(m_TrackingPriorHandler); } 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); } 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 (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->SetParameters(params); m_Tracker->SetTrackingHandler(m_TrackingHandler); m_Tracker->SetVerbose(!m_Controls->m_InteractiveBox->isChecked()); m_ParentNode = m_InputImageNodes.at(0); m_TrackingThread.start(QThread::LowestPriority); }